iaf_chs_2007

class brainpy.state.iaf_chs_2007(in_size, tau_epsp=Quantity(8.5, 'ms'), tau_reset=Quantity(15.4, 'ms'), V_epsp=0.77, V_reset=2.31, V_noise=0.0, noise=None, gsl_error_tol=1e-06, V_initializer=Constant(value=0.0), spk_fun=ReluGrad(alpha=0.3, width=1.0), spk_reset='hard', name=None)

NEST-compatible iaf_chs_2007 spike-response neuron model.

Description

iaf_chs_2007 is a discrete-time linear spike-response neuron model developed for analyzing thalamic filtering of retinal spike trains (Carandini, Horton, and Sincich, 2007). The normalized membrane potential is the sum of three components:

• An alpha-shaped postsynaptic potential (PSP) waveform V_syn driven by excitatory spike inputs,

• A post-spike reset/after-hyperpolarization (AHP) component V_spike that decays exponentially after each spike emission,

• An optional externally prepared noise trace scaled by V_noise.

This implementation mirrors NEST’s C++ model models/iaf_chs_2007.{h,cpp} to ensure semantic equivalence: exact discrete-time update order, normalized voltage conventions (\(U_{th} = 1\), \(E_L = 0\)), positive-weight-only PSP accumulation, and external noise injection without refractory state.

1. Model equations and exact integration

Let \(h = dt\) (in ms) denote the global integration step, and \(k\) index discrete time steps. NEST precomputes four propagators at initialization:

\[P_{11} = e^{-h / \tau_{\mathrm{epsp}}}, \qquad P_{22} = P_{11}, \qquad P_{30} = e^{-h / \tau_{\mathrm{reset}}},\]

\[P_{21} = U_{\mathrm{epsp}} \, e \, P_{11} \, \frac{h}{\tau_{\mathrm{epsp}}}.\]

Here \(U_{\mathrm{epsp}}\) sets the EPSP peak amplitude and \(P_{21}\) is derived from the alpha-function kernel integral.

Each simulation step updates state in the following order:

\[V_{\mathrm{syn}}^{k+1} = P_{22} \, V_{\mathrm{syn}}^k + P_{21} \, i_{\mathrm{syn}}^k,\]

\[i_{\mathrm{syn}}^{k+1} = P_{11} \, i_{\mathrm{syn}}^k + \max(w_k, 0),\]

\[V_{\mathrm{spike}}^{k+1} = P_{30} \, V_{\mathrm{spike}}^k,\]

\[V_m^{k+1} = V_{\mathrm{syn}}^{k+1} + V_{\mathrm{spike}}^{k+1} + U_{\mathrm{noise}} \, \eta_k,\]

where \(w_k = \sum_{j} w_j \delta(t - t_j^{spike})\) collects incoming excitatory spike weights delivered at step \(k\), and \(\eta_k\) is the externally provided noise sample (if configured).

Spike emission uses the hard threshold \(U_{\mathrm{th}} = 1\):

\[V_m^{k+1} \ge 1 \quad \Longrightarrow \quad V_{\mathrm{spike}}^{k+1} \leftarrow V_{\mathrm{spike}}^{k+1} - U_{\mathrm{reset}}, \quad V_m^{k+1} \leftarrow V_m^{k+1} - U_{\mathrm{reset}}.\]

Both the reset/AHP component and the total membrane potential are decremented by \(U_{\mathrm{reset}}\) upon spike emission. No refractory clamping occurs; the neuron can spike again immediately if \(V_m\) remains above threshold.

2. Update ordering and NEST semantics

The per-step operation sequence is identical to NEST’s C++ update() routine:

1. Update V_syn from the previous step’s i_syn_ex.

2. Decay i_syn_ex by P11.

3. Add arriving excitatory spike weights (non-negative values only; negative inputs are clamped to zero).

4. Decay V_spike by P30.

5. Sample and add noise term if V_noise > 0 and noise buffer is non-empty.

6. Compute total V_m and apply threshold/reset/spike emission logic.

A critical consequence of this ordering is that a spike arriving in the current step immediately increments i_syn_ex, but the resulting V_syn contribution appears in V_m only from the next step onward (one-step synaptic delay).

3. Noise semantics

NEST expects noise to be externally prepared (e.g., from a Gaussian distribution with specified variance) and supplied as a pre-generated sequence. This implementation follows the same convention:

• Noise is active only if V_noise > 0. and the noise buffer is non-empty.

• One sample per neuron per step is consumed from the flat noise array using a position index.

• If the noise buffer is exhausted before the end of the simulation, an IndexError is raised. Users must provide a noise array of length at least equal to the total number of simulation steps.

4. Assumptions, constraints, and computational complexity

• All model parameters are scalar or broadcastable to self.varshape.

• Construction-time constraints enforce V_epsp >= 0, V_reset >= 0, tau_epsp > 0, tau_reset > 0 elementwise.

• The model operates in normalized voltage units where \(E_L = 0\) (rest), \(U_{th} = 1\) (threshold).

• Negative input weights are silently clamped to zero (matching NEST’s positive-weight-only convention for this model).

• Unlike standard LIF models, there is no refractory state or explicit continuous current input handler; the update(x=...) argument is unused by design.

• Per-step complexity is \(O(|\mathrm{state}|)\) for state propagation, plus \(O(K)\) for collecting K delta inputs.

Parameters:

• in_size (Size) – Population shape specification. Model parameters and states are broadcast to self.varshape derived from in_size.

• tau_epsp (ArrayLike, optional) – EPSP time constant \(\tau_{\mathrm{epsp}}\) in ms, broadcastable to self.varshape. Must be strictly positive elementwise. Controls the rise/decay timescale of the alpha-shaped PSP kernel. Default is 8.5 * u.ms.

• tau_reset (ArrayLike, optional) – Post-spike reset/AHP time constant \(\tau_{\mathrm{reset}}\) in ms, broadcastable to self.varshape. Must be strictly positive elementwise. Governs exponential decay of the V_spike component after each spike. Default is 15.4 * u.ms.

• V_epsp (ArrayLike, optional) – Normalized maximal EPSP amplitude \(U_{\mathrm{epsp}}\) (dimensionless), broadcastable to self.varshape. Must be non-negative elementwise. Sets the peak height of the PSP waveform per unit weight. Default is 0.77.

• V_reset (ArrayLike, optional) – Normalized reset/AHP magnitude \(U_{\mathrm{reset}}\) (dimensionless), broadcastable to self.varshape. Must be non-negative elementwise. Both V_spike and V_m are decremented by this value upon threshold crossing. Default is 2.31.

• V_noise (ArrayLike, optional) – Noise scale factor \(U_{\mathrm{noise}}\) (dimensionless), broadcastable to self.varshape. Multiplies externally provided noise samples. No noise is added if V_noise == 0 or if the noise buffer is empty. Default is 0.0.

• noise (Sequence[float] | np.ndarray | None, optional) – Externally prepared noise samples (dimensionless). If provided, must be a flat 1D sequence of at least num_steps values, where num_steps is the total simulation duration in steps. One sample per neuron per step is consumed sequentially. If None, no noise is applied. Default is None.

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

• V_initializer (Callable, optional) – Initializer used by init_state() for membrane potential V. Must return dimensionless values with shape compatible with self.varshape (and optional batch prefix). Default is braintools.init.Constant(0.0).

• spk_fun (Callable, optional) – Surrogate spike function used by get_spike() and update(). Receives normalized threshold distance tensor. Default is braintools.surrogate.ReluGrad().

• spk_reset (str, optional) – Reset policy forwarded to Neuron. 'hard' matches NEST’s hard subtraction reset. Default is 'hard'.

• name (str or None, optional) – Optional node name.

Parameter Mapping

Table 11 Parameter mapping to model symbols

Parameter	Type / shape / unit	Default	Math symbol	Semantics
in_size	Size; scalar/tuple	required	–	Defines self.varshape for parameter/state broadcasting.
tau_epsp	ArrayLike, broadcastable (ms), > 0	8.5 * u.ms	\(\tau_{\mathrm{epsp}}\)	Alpha-kernel time constant for EPSP waveform.
tau_reset	ArrayLike, broadcastable (ms), > 0	15.4 * u.ms	\(\tau_{\mathrm{reset}}\)	Exponential decay time constant for post-spike AHP.
V_epsp	ArrayLike, broadcastable (dimensionless), >= 0	0.77	\(U_{\mathrm{epsp}}\)	Normalized peak EPSP amplitude per unit weight.
V_reset	ArrayLike, broadcastable (dimensionless), >= 0	2.31	\(U_{\mathrm{reset}}\)	Normalized reset/AHP decrement applied on spike.
V_noise	ArrayLike, broadcastable (dimensionless)	0.0	\(U_{\mathrm{noise}}\)	Noise scale factor; zero disables noise injection.
noise	Sequence[float] | np.ndarray | None (dimensionless)	None	\(\eta_k\)	Externally prepared noise samples; must have length >= num_steps.
gsl_error_tol	ArrayLike, broadcastable, unitless, > 0	1e-6	–	Local absolute tolerance for the embedded RKF45 error estimate.
V_initializer	Callable returning (dimensionless)	Constant(0.0)	\(V_m(0)\)	Initial membrane potential at t=0.
spk_fun	Callable	ReluGrad()	–	Differentiable surrogate for spike emission.
spk_reset	str ('hard' or 'soft')	'hard'	–	Reset mode; 'hard' matches NEST semantics.

Raises:

• ValueError –

  • If V_epsp < 0 or V_reset < 0 elementwise. - If tau_epsp <= 0 or tau_reset <= 0 elementwise.

• IndexError – If the noise buffer is exhausted before the end of the simulation. Ensure the noise array has length at least equal to the total number of simulation steps.

Notes

Key differences from standard LIF models:

• Operates in normalized voltage units where \(E_L = 0\) and \(U_{th} = 1\). There is no explicit leak conductance or external current integration beyond the PSP and AHP components.

• No refractory state: The neuron can spike immediately again if \(V_m\) remains above threshold after reset.

• Positive-weight-only synapses: Negative incoming spike weights are clamped to zero, matching NEST’s convention for this model.

• External noise injection: Noise is not generated internally but must be pre-prepared and supplied as a flat array. The model consumes one sample per neuron per step from the noise buffer using a sequential index.

• No continuous current input: Unlike iaf_psc_exp and similar models, iaf_chs_2007 in NEST has no CurrentEvent handler. The update(x=...) argument is present for API compatibility but is intentionally unused.

Spike-response interpretation:

The model is a discrete-time linear spike-response model (SRM) where each incoming spike triggers an alpha-shaped PSP and each output spike triggers an exponential AHP. The total membrane potential is the linear sum of these components plus optional noise. This differs from integrate-and-fire models that compute a continuous leak current; here, the “leak” is implicit in the exponential decay of the PSP and AHP kernels.

Computational considerations

• State propagation uses exact exponential integration with precomputed propagators P11, P22, P30, and P21, ensuring machine-precision accuracy regardless of dt (within numerical precision of exp and floating-point arithmetic).

• The model performs all computations in NumPy on the host before transferring results to JAX, which may limit GPU acceleration efficiency. This design choice ensures exact NEST compatibility, including identical floating-point rounding behavior.

• For large-scale simulations, consider using the standard LIF or ExpIF models, which are fully JIT-compatible and GPU-accelerated.

References

[1]

Carandini M, Horton JC, Sincich LC (2007). Thalamic filtering of retinal spike trains by postsynaptic summation. Journal of Vision 7(14):20, 1-11. DOI: https://doi.org/10.1167/7.14.20

[2]

Rotter S, 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

Examples

Create a population of 100 neurons and drive them with Poisson spike input:

>>> import brainpy.state as bp
>>> import saiunit as u
>>> import brainstate as bs
>>> import brainevent as be
>>> import numpy as np
>>>
>>> # Prepare external noise (Gaussian)
>>> num_steps = 1000
>>> noise = np.random.randn(num_steps)
>>>
>>> # Create model
>>> model = bp.iaf_chs_2007(
...     in_size=100,
...     tau_epsp=8.5 * u.ms,
...     tau_reset=15.4 * u.ms,
...     V_epsp=0.77,
...     V_reset=2.31,
...     V_noise=0.1,
...     noise=noise
... )
>>>
>>> # Poisson spike source
>>> poisson = be.nn.PoissonEncoder(in_size=100)
>>>
>>> # Projection
>>> proj = bp.AlignPostProj(
...     comm=be.nn.AllToAll(pre_size=100, post_size=100, w_init=0.05),
...     out=bp.CUBA(),
...     post=model
... )
>>>
>>> # Simulate
>>> with bs.environ.context(dt=0.1 * u.ms):
...     model.init_state()
...     for _ in range(num_steps):
...         inp = poisson.generate()
...         proj(inp)
...         spk = model.update()
...

get_spike(V=None)

Compute spike output using the surrogate spike function.

Applies the surrogate gradient function spk_fun to the normalized threshold distance \(V_m - U_{th}\), where \(U_{th} = 1\). This produces a differentiable spike output suitable for gradient-based training.

Parameters:

V (ArrayLike or None, optional) – Membrane potential tensor (dimensionless), with shape matching self.varshape (plus optional batch prefix). If None, uses self.V.value (the current membrane state). Default is None.

Returns:

Spike output tensor with the same shape as V. Values are typically in [0, 1] or continuous approximations depending on the chosen spk_fun.

Return type:

ArrayLike

Notes

Unlike typical LIF models where \(U_{reset} < U_{th}\), this model has \(U_{reset} > U_{th}\) by default (2.31 vs 1.0). The spike function operates on \(V_m - U_{th}\) to detect threshold crossings; the surrogate gradient enables backpropagation through spike events.

init_state(**kwargs)

Initialize all state variables for the neuron population.

Creates and registers the following states with brainstate:

• i_syn_ex: Excitatory synaptic current state (ShortTermState, initialized to zeros).

• V_syn: EPSP waveform state (ShortTermState, initialized to zeros).

• V_spike: Post-spike reset/AHP state (ShortTermState, initialized to zeros).

• V: Normalized membrane potential (HiddenState, initialized via V_initializer).

• position: Current index into the noise buffer (ShortTermState, initialized to zeros as int64).

• last_spike_time: Last spike time in ms (ShortTermState, initialized to -1e7 * u.ms to indicate no previous spike).

Parameters:

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

Notes

All states except V are initialized to zero. The membrane potential V is initialized using self.V_initializer, which defaults to Constant(0.0) (normalized resting potential).

reset_state(batch_size=None, **kwargs)

Reset all state variables to their initial values.

Resets all state variables to the same values as init_state(), but operates on existing state instances rather than creating new ones. This is useful for re-initializing a model between simulation runs without destroying the computational graph.

Parameters:

• batch_size (int or None, optional) – Optional batch dimension size. If provided, states will have shape (batch_size, *self.varshape). If None, states have shape self.varshape.

• **kwargs – Additional keyword arguments (currently unused).

Notes

The position index into the noise buffer is reset to zero. If you want to continue consuming the noise sequence from where it left off, manually preserve and restore self.position.value before/after calling this method.

update(x=0.0)

Advance the neuron state by one simulation time step.

Implements the discrete-time update rule following NEST’s exact sequence using precomputed exact propagators:

1. Update V_syn from previous i_syn_ex using propagator P21 (one-step synaptic delay).

2. Decay i_syn_ex by P11.

3. Decay V_spike by P30.

4. Sample noise term if V_noise > 0 and noise buffer is non-empty.

5. Add arriving excitatory delta inputs (non-negative weights only) to i_syn_ex.

6. Compute total membrane potential V_m = V_syn + V_spike + noise and apply threshold/reset logic.

7. Return spike output.

Parameters:

x (ArrayLike, optional) – Continuous input current (unused by design). Default is 0.0.

Returns:

Spike output tensor with shape matching self.V.value.

Return type:

ArrayLike

Raises:

IndexError – If the noise buffer is exhausted before the end of the simulation.