mat2_psc_exp

class brainpy.state.mat2_psc_exp(in_size, E_L=Quantity(-70., "mV"), C_m=Quantity(100., "pF"), tau_m=Quantity(5., "ms"), t_ref=Quantity(2., "ms"), tau_syn_ex=Quantity(1., "ms"), tau_syn_in=Quantity(3., "ms"), I_e=Quantity(0., "pA"), tau_1=Quantity(10., "ms"), tau_2=Quantity(200., "ms"), alpha_1=Quantity(37., "mV"), alpha_2=Quantity(2., "mV"), omega=Quantity(-51., "mV"), V_initializer=Constant(value=-70. mV), spk_fun=ReluGrad(alpha=0.3, width=1.0), spk_reset='hard', ref_var=False, name=None)

NEST-compatible mat2_psc_exp neuron model.

Non-resetting leaky integrate-and-fire neuron model with exponential postsynaptic currents and a two-timescale adaptive threshold.

1. Model Overview

mat2_psc_exp implements a leaky integrate-and-fire model with exponential shaped postsynaptic currents (PSCs) and a multi-timescale adaptive threshold (MAT) [3]. Key features:

• No voltage reset: The membrane potential continues to integrate through spikes, providing biological realism for high-firing-rate regimes.

• Two-timescale threshold adaptation: Separate fast (τ₁) and slow (τ₂) threshold components capture short-term spike frequency adaptation and long-term gain control.

• Absolute refractory period: Prevents multiple spikes within t_ref without clamping the membrane potential.

• Exact integration: Subthreshold dynamics use the exponential Euler propagator [1] for numerical stability.

2. Mathematical Formulation

2.1 Subthreshold Membrane Dynamics

The membrane potential evolves according to:

\[\frac{dV_m}{dt} = -\frac{V_m - E_L}{\tau_m} + \frac{I_{\mathrm{syn,ex}} + I_{\mathrm{syn,in}} + I_e + I_0}{C_m}\]

where:

• \(V_m\) – membrane potential (absolute voltage)

• \(E_L\) – resting potential

• \(\tau_m = C_m / g_L\) – membrane time constant

• \(I_{\mathrm{syn,ex}}, I_{\mathrm{syn,in}}\) – synaptic currents

• \(I_e\) – constant external current

• \(I_0\) – buffered step current input (updated each time step)

2.2 Synaptic Currents

Exponentially decaying currents with infinitely fast rise:

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

Incoming spike weights are added instantaneously: \(I_{\mathrm{syn}} \leftarrow I_{\mathrm{syn}} + w\).

2.3 Adaptive Threshold

The effective spike threshold is the sum of a baseline and two decaying components:

\[V_{th}(t) = \omega + V_{th,1}(t) + V_{th,2}(t)\]

where:

\[\frac{dV_{th,1}}{dt} = -\frac{V_{th,1}}{\tau_1} \qquad \frac{dV_{th,2}}{dt} = -\frac{V_{th,2}}{\tau_2}\]

On each spike at time \(t_{\text{spike}}\):

\[V_{th,1}(t_{\text{spike}}^+) = V_{th,1}(t_{\text{spike}}^-) + \alpha_1 \qquad V_{th,2}(t_{\text{spike}}^+) = V_{th,2}(t_{\text{spike}}^-) + \alpha_2\]

2.4 Spike Emission

A spike is emitted when:

\[V_m \geq \omega + V_{th,1} + V_{th,2} \quad \text{and} \quad t - t_{\text{last\_spike}} \geq t_{\text{ref}}\]

After spiking:

• Threshold components jump by \(\alpha_1, \alpha_2\)

• Refractory counter is set to \(\lceil t_{\text{ref}} / \Delta t \rceil\)

• Membrane potential is NOT reset (continues integrating)

3. Numerical Integration

The model uses exact integration for the linear subthreshold system [1]. For a time step \(h = \Delta t\):

3.1 Membrane Potential Propagators

\[\begin{split}V_m(t+h) &= V_m(t) e^{-h/\tau_m} + E_L (1 - e^{-h/\tau_m}) \\ &\quad + P_{21}^{\text{ex}} I_{\text{syn,ex}}(t) + P_{21}^{\text{in}} I_{\text{syn,in}}(t) + P_{20} (I_e + I_0)\end{split}\]

where:

\[\begin{split}P_{21}^{\text{ex}} &= -\frac{\tau_m}{C_m (1 - \tau_m/\tau_{\text{syn,ex}})} e^{-h/\tau_{\text{syn,ex}}} (e^{h(1/\tau_{\text{syn,ex}} - 1/\tau_m)} - 1) \\ P_{21}^{\text{in}} &= -\frac{\tau_m}{C_m (1 - \tau_m/\tau_{\text{syn,in}})} e^{-h/\tau_{\text{syn,in}}} (e^{h(1/\tau_{\text{syn,in}} - 1/\tau_m)} - 1) \\ P_{20} &= -\frac{\tau_m}{C_m} (e^{-h/\tau_m} - 1)\end{split}\]

3.2 Synaptic and Threshold Propagators

\[\begin{split}I_{\text{syn}}(t+h) &= I_{\text{syn}}(t) e^{-h/\tau_{\text{syn}}} + w_{\text{spike}} \\ V_{th,1}(t+h) &= V_{th,1}(t) e^{-h/\tau_1} \\ V_{th,2}(t+h) &= V_{th,2}(t) e^{-h/\tau_2}\end{split}\]

3.3 Numerical Stability Constraint

The propagators become ill-conditioned when \(\tau_m \approx \tau_{\text{syn,ex}}\) or \(\tau_m \approx \tau_{\text{syn,in}}\) due to division by \((1 - \tau_m/\tau_{\text{syn}})\). The constructor enforces strict inequality.

4. Update Order (NEST-Compatible)

For each time step (matching NEST’s mat2_psc_exp.cpp):

1. Integrate membrane potential using exact propagators

2. Decay adaptive threshold components (\(V_{th,1}\), \(V_{th,2}\))

3. Decay synaptic currents and add incoming spike weights

4. Detect spikes: if not refractory and \(V_m \geq V_{th}\), emit spike, jump threshold components, reset refractory counter

5. Update refractory state: decrement counter if active

6. Buffer current inputs for next step (\(I_0\))

5. Surrogate Gradient Handling

For differentiable training, the output spike signal passes through a surrogate gradient function (spk_fun). The voltage is scaled relative to the adaptive threshold:

\[v_{\text{scaled}} = \frac{V_m - V_{th}}{|\omega - E_L|}\]

where the denominator provides a characteristic voltage scale (~19 mV with defaults).

Parameters:

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

• E_L (Quantity, ArrayLike, optional) – Resting membrane potential (default: -70 mV). Broadcastable to in_size.

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

• tau_m (Quantity, ArrayLike, optional) – Membrane time constant (default: 5 ms). Must be strictly positive and differ from tau_syn_ex and tau_syn_in to avoid numerical degeneracy.

• t_ref (Quantity, ArrayLike, optional) – Duration of absolute refractory period (default: 2 ms). Must be strictly positive.

• tau_syn_ex (Quantity, ArrayLike, optional) – Time constant of excitatory postsynaptic current (default: 1 ms). Must be strictly positive and differ from tau_m.

• tau_syn_in (Quantity, ArrayLike, optional) – Time constant of inhibitory postsynaptic current (default: 3 ms). Must be strictly positive and differ from tau_m.

• I_e (Quantity, ArrayLike, optional) – Constant external input current (default: 0 pA). Broadcastable to in_size.

• tau_1 (Quantity, ArrayLike, optional) – Short time constant of adaptive threshold (default: 10 ms). Must be strictly positive.

• tau_2 (Quantity, ArrayLike, optional) – Long time constant of adaptive threshold (default: 200 ms). Must be strictly positive.

• alpha_1 (Quantity, ArrayLike, optional) – Amplitude of short-timescale threshold jump on spike (default: 37 mV).

• alpha_2 (Quantity, ArrayLike, optional) – Amplitude of long-timescale threshold jump on spike (default: 2 mV).

• omega (Quantity, ArrayLike, optional) – Resting spike threshold (default: -51 mV). This is an absolute voltage, not relative to E_L. With defaults, the threshold is 19 mV above resting.

• V_initializer (Callable, optional) – Initializer for membrane potential (default: Constant(-70 mV)). Called as V_initializer(shape, batch_size) to produce initial voltages.

• spk_fun (Callable, optional) – Surrogate gradient function for spike generation (default: ReluGrad()). Maps scaled voltage to differentiable spike signal.

• spk_reset (str, optional) – Reset mode for spike output (default: 'hard'). Options: 'hard' (stop gradient) or 'soft' (preserve gradient). Does NOT affect membrane voltage reset (which never occurs in this model).

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

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

Parameter Mapping

Correspondence between constructor parameters and mathematical symbols:

Parameter	Default	Math Symbol	Description
in_size	(required)	—	Population shape
E_L	-70 mV	\(E_L\)	Resting membrane potential
C_m	100 pF	\(C_m\)	Membrane capacitance
tau_m	5 ms	\(\tau_m\)	Membrane time constant
t_ref	2 ms	\(t_{\text{ref}}\)	Duration of absolute refractory period
tau_syn_ex	1 ms	\(\tau_{\text{syn,ex}}\)	Time constant of excitatory PSC
tau_syn_in	3 ms	\(\tau_{\text{syn,in}}\)	Time constant of inhibitory PSC
I_e	0 pA	\(I_e\)	Constant external input current
tau_1	10 ms	\(\tau_1\)	Short time constant of adaptive threshold
tau_2	200 ms	\(\tau_2\)	Long time constant of adaptive threshold
alpha_1	37 mV	\(\alpha_1\)	Amplitude of short-timescale threshold jump
alpha_2	2 mV	\(\alpha_2\)	Amplitude of long-timescale threshold jump
omega	-51 mV	\(\omega\)	Resting spike threshold (absolute voltage)
V_initializer	Constant(-70 mV)	—	Membrane potential initializer
spk_fun	ReluGrad()	—	Surrogate spike function
spk_reset	'hard'	—	Reset mode (for gradient handling)
ref_var	False	—	If True, expose refractory boolean state

State Variables

After init_state(), the following state variables are available:

Variable	Type	Description
V	HiddenState (mV)	Membrane potential (absolute voltage)
V_th_1	ShortTermState	Short-timescale adaptive threshold component (mV)
V_th_2	ShortTermState	Long-timescale adaptive threshold component (mV)
i_syn_ex	ShortTermState	Excitatory postsynaptic current (pA)
i_syn_in	ShortTermState	Inhibitory postsynaptic current (pA)
i_0	ShortTermState	Buffered DC input current (pA, one-step delayed)
refractory_step_count	ShortTermState	Refractory countdown (integer steps remaining)
last_spike_time	ShortTermState	Time of last spike (ms)
refractory	ShortTermState	Boolean refractory flag (only if ref_var=True)

Raises:

• ValueError – If C_m <= 0, tau_m <= 0, tau_syn_ex <= 0, tau_syn_in <= 0, t_ref <= 0, tau_1 <= 0, or tau_2 <= 0.

• ValueError – If tau_m == tau_syn_ex or tau_m == tau_syn_in (numerical degeneracy).

See also

amat2_psc_exp

Variant with exponential synaptic currents and additional state variables.

iaf_psc_exp

Standard LIF with voltage reset (no adaptive threshold).

aeif_psc_exp

Adaptive exponential IF with spike-triggered adaptation current.

Notes

Biological Interpretation: The MAT model captures spike frequency adaptation without explicit adaptation currents. The fast threshold component (τ₁ ~ 10 ms) models sodium channel inactivation, while the slow component (τ₂ ~ 200 ms) models calcium-dependent potassium currents.

Comparison to NEST: This implementation matches NEST’s mat2_psc_exp update order (see NEST 3.7+ mat2_psc_exp.cpp). Key differences:

• Surrogate gradients: brainpy.state adds differentiable spike signals via spk_fun for gradient-based learning; NEST uses exact spike times.

• Batch dimension: brainpy.state supports batch processing for parallel simulations; NEST operates on single neuron instances.

• Precision: brainpy.state uses float32 (JAX default); NEST uses float64. Minor numerical differences may occur for long simulations.

Performance Notes:

• Propagator computation (exponentials, expm1) dominates runtime for small populations.

• For large populations (>10k neurons), vectorized operations amortize this cost.

• Use jax.jit compilation for optimal performance.

References

[1] (1,2)

Rotter S and Diesmann M (1999). Exact simulation of time-invariant linear systems with applications to neuronal modeling. Biological Cybernetics 81:381-402. DOI: https://doi.org/10.1007/s004220050570

[2]

Diesmann M, Gewaltig M-O, Rotter S, Aertsen A (2001). State space analysis of synchronous spiking in cortical neural networks. Neurocomputing 38-40:565-571. DOI: https://doi.org/10.1016/S0925-2312(01)00409-X

[3]

Kobayashi R, Tsubo Y and Shinomoto S (2009). Made-to-order spiking neuron model equipped with a multi-timescale adaptive threshold. Frontiers in Computational Neuroscience 3:9. DOI: https://doi.org/10.3389/neuro.10.009.2009

Examples

Basic Usage:

>>> import brainpy.state as bp
>>> import saiunit as u
>>> # Create a population of 100 MAT neurons
>>> neurons = bp.mat2_psc_exp(100, tau_1=10*u.ms, tau_2=200*u.ms)
>>> neurons.init_all_states()
>>> # Inject step current and simulate
>>> with brainstate.environ.context(dt=0.1*u.ms):
...     spikes = neurons.update(500*u.pA)  # 500 pA step current

Demonstrating Adaptive Threshold:

>>> # Single neuron with strong adaptation
>>> neuron = bp.mat2_psc_exp(1, alpha_1=50*u.mV, alpha_2=5*u.mV)
>>> neuron.init_all_states()
>>> with brainstate.environ.context(dt=0.1*u.ms):
...     V_trace = []
...     for _ in range(1000):  # 100 ms simulation
...         spk = neuron.update(800*u.pA)
...         V_trace.append(neuron.V.value)
>>> # Plot V_trace to observe spike frequency adaptation

Network with Excitatory/Inhibitory Synapses:

>>> exc = bp.mat2_psc_exp(800)
>>> inh = bp.mat2_psc_exp(200, tau_syn_ex=0.5*u.ms)
>>> exc.init_all_states()
>>> inh.init_all_states()
>>> # Connect populations via projections (see brainpy.state.Projection)

get_spike(V=None, V_th=None)

Compute surrogate gradient spike signal.

Parameters:

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

• V_th (Quantity, ArrayLike, optional) – Effective threshold (mV). If None, computes as omega + V_th_1 + V_th_2.

Returns:

spike – Differentiable spike signal from surrogate function. Shape matches V.

Return type:

ArrayLike

Notes

The voltage is scaled relative to the adaptive threshold before passing through the surrogate function, providing a normalized input that improves gradient stability.

init_state(**kwargs)

State initialization function.

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

Advance the neuron state by one time step.

Implements the NEST-compatible update order for the MAT2 model with exact integration of subthreshold dynamics.

Parameters:

• x (Quantity, ArrayLike, optional) – External input current (pA) for this time step (default: 0 pA). Broadcastable to population shape. This current is buffered and applied in the next time step (one-step delay).

• spike_delta (Quantity, optional) – Instantaneous spike-weight input (pA) to add to synaptic currents. When provided, bypasses sum_delta_inputs() — useful for JIT-compiled brainstate.transform.for_loop simulations where delta inputs are pre-computed as a JAX array indexed by step. Positive values go to i_syn_ex; negative values go to i_syn_in.

Returns:

spike – Differentiable spike signal for this time step. Shape matches population size.

Return type:

ArrayLike

Notes

Update sequence (NEST-compatible):

1. Integrate membrane potential using exact propagators

2. Decay adaptive threshold components (V_th_1, V_th_2)

3. Decay synaptic currents and add incoming spike weights

4. Detect spikes: if not refractory and V_m >= V_th, emit spike

5. On spike: jump threshold components, reset refractory counter

6. Buffer external current for next step

Key implementation details:

• Membrane potential is never reset (non-resetting LIF)

• Spike detection compares V_m against the adaptive threshold V_th = ω + V_th_1 + V_th_2

• Refractory period is implemented as an integer countdown; no voltage clamping

• External current x is stored in i_0 and applied in the next time step

Numerical stability:

The exact integration scheme requires tau_m ≠ tau_syn_ex and tau_m ≠ tau_syn_in. Violations of this constraint are caught during initialization.