iaf_psc_alpha_multisynapse

class brainpy.state.iaf_psc_alpha_multisynapse(in_size, E_L=Quantity(-70., "mV"), C_m=Quantity(250., "pF"), tau_m=Quantity(10., "ms"), t_ref=Quantity(2., "ms"), V_th=Quantity(-55., "mV"), V_reset=Quantity(-70., "mV"), tau_syn=Quantity([2.], "ms"), I_e=Quantity(0., "pA"), V_min=None, 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 iaf_psc_alpha_multisynapse neuron model.

Current-based leaky integrate-and-fire neuron with an arbitrary number of receptor-indexed alpha-shaped synaptic current channels.

Description

iaf_psc_alpha_multisynapse mirrors NEST models/iaf_psc_alpha_multisynapse.{h,cpp} and generalizes iaf_psc_alpha from two fixed excitatory/inhibitory channels to n_receptors independently parameterized current ports.

Each receptor k (1-based, NEST convention) carries its own alpha time constant tau_syn[k-1]. Synaptic weights are signed currents in pA; positive values are depolarizing and negative values are hyperpolarizing.

1. Continuous-Time Dynamics and Receptor States

Membrane dynamics are

\[\frac{dV_m}{dt} = -\frac{V_m - E_L}{\tau_m} + \frac{\sum_k I_k + I_e + I_0}{C_m},\]

where \(I_0\) is the one-step delayed continuous current buffer (NEST ring-buffer semantics) and \(I_k\) are the per-receptor alpha currents.

For each receptor \(k\), the alpha current kernel is represented by a two-state linear system (y1[k], y2[k]):

\[\frac{d\,y1_k}{dt} = -\frac{y1_k}{\tau_{\mathrm{syn},k}}, \qquad \frac{d\,y2_k}{dt} = y1_k - \frac{y2_k}{\tau_{\mathrm{syn},k}}.\]

The effective synaptic current for receptor \(k\) is \(I_k = y2_k\). An incoming spike with weight \(w_k\) (pA) is injected into y1[k] with the NEST alpha normalization factor:

\[y1_k \leftarrow y1_k + \frac{e}{\tau_{\mathrm{syn},k}} w_k.\]

This normalization ensures that a single spike with weight \(w_k\) produces a current kernel that peaks exactly at \(w_k\) when \(t = \tau_{\mathrm{syn},k}\):

\[I_k(t) = w_k \frac{t}{\tau_{\mathrm{syn},k}} \exp\!\left(1 - \frac{t}{\tau_{\mathrm{syn},k}}\right), \quad t \ge 0.\]

2. Exact Discrete Propagator, Derivation Constraints, and Stability

With fixed step \(h = dt\), exact matrix propagation of the linear subsystem is used. For each receptor \(k\):

\[y1_{k,n+1} = P_{11,k}\,y1_{k,n} + \frac{e}{\tau_{\mathrm{syn},k}} w_{k,n},\]

\[y2_{k,n+1} = P_{21,k}\,y1_{k,n} + P_{22,k}\,y2_{k,n},\]

where \(P_{11,k} = P_{22,k} = e^{-h/\tau_{\mathrm{syn},k}}\) and \(P_{21,k} = h\,e^{-h/\tau_{\mathrm{syn},k}}\).

Membrane relative voltage \(y_3 = V_m - E_L\) is updated as

\[y_{3,n+1} = P_{33}\,y_{3,n} + P_{30}(I_{0,n} + I_e) + \sum_k \left(P_{31,k}\,y1_{k,n} + P_{32,k}\,y2_{k,n}\right),\]

with \(P_{33} = e^{-h/\tau_m}\) and \(P_{30} = \tau_m(1 - e^{-h/\tau_m})/C_m\). Coefficients \(P_{31,k}\), \(P_{32,k}\) are computed via iaf_psc_alpha._alpha_propagator_p31_p32(), which applies the stable near-singular limit for \(\tau_m \approx \tau_{\mathrm{syn},k}\):

\[P_{32}^{\mathrm{sing}} = \frac{h}{C_m} e^{-h/\tau_m}, \qquad P_{31}^{\mathrm{sing}} = \frac{h^2}{2C_m} e^{-h/\tau_m},\]

preventing catastrophic cancellation when \(\tau_m = \tau_{\mathrm{syn},k}\).

3. Update Order per Simulation Step (NEST Semantics)

Per-step execution order:

1. Integrate membrane with exact propagator for neurons not refractory (\(r = 0\)).

2. Decrement refractory counters for neurons currently refractory (\(r > 0\)).

3. Propagate all receptor alpha states y1, y2 forward by one step.

4. Inject receptor-specific spike weights into y1, including default delta input mapped to receptor 1 when n_receptors > 0.

5. Apply threshold test, hard reset, refractory assignment, and spike emission.

6. Store buffered continuous current for the next step.

4. Assumptions, Constraints, and Computational Implications

• C_m > 0, tau_m > 0, all tau_syn > 0, t_ref >= 0, and V_reset < V_th are enforced at construction.

• update(x=...) uses one-step delayed current buffering: current provided at step n contributes through i_const at step n+1, matching NEST ring-buffer event semantics.

• The update path is fully vectorized over self.varshape and scales as \(O(\prod \mathrm{varshape} \times n\_receptors)\) per call.

• Internal propagator arithmetic is performed in NumPy float64 before writing back to BrainUnit-typed states.

• When n_receptors == 0, all spike event inputs are silently ignored.

Parameters:

• in_size (Size) – Population shape specification. Per-neuron parameters and state variables are broadcast/initialized over self.varshape derived from in_size.

• E_L (ArrayLike, optional) – Resting potential \(E_L\) in mV; scalar or array broadcastable to self.varshape. Default is -70. * u.mV.

• C_m (ArrayLike, optional) – Membrane capacitance \(C_m\) in pF; broadcastable to self.varshape and strictly positive. Default is 250. * u.pF.

• tau_m (ArrayLike, optional) – Membrane time constant \(\tau_m\) in ms; broadcastable and strictly positive. Default is 10. * u.ms.

• t_ref (ArrayLike, optional) – Absolute refractory period \(t_\mathrm{ref}\) in ms; broadcastable and nonnegative. Converted to integer grid steps by ceil(t_ref / dt). Default is 2. * u.ms.

• V_th (ArrayLike, optional) – Spike threshold \(V_\mathrm{th}\) in mV; broadcastable to self.varshape. Default is -55. * u.mV.

• V_reset (ArrayLike, optional) – Post-spike reset potential \(V_\mathrm{reset}\) in mV; broadcastable and constrained by V_reset < V_th elementwise. Default is -70. * u.mV.

• tau_syn (ArrayLike, optional) – Receptor alpha time constants in ms. Values are converted to a 1-D float64 array with shape (n_receptors,); every entry must be strictly positive. The number of entries defines n_receptors. Default is (2.0,) * u.ms (one receptor).

• I_e (ArrayLike, optional) – Constant injected current \(I_e\) in pA; scalar or array broadcastable to self.varshape. Default is 0. * u.pA.

• V_min (ArrayLike or None, optional) – Optional lower clamp \(V_\mathrm{min}\) in mV applied to the membrane candidate update before thresholding. None disables clamping. Default is None.

• V_initializer (Callable, optional) – Initializer for membrane state V used by init_state(). Default is braintools.init.Constant(-70. * u.mV).

• spk_fun (Callable, optional) – Surrogate spike function used by get_spike() and returned by update(). Default is braintools.surrogate.ReluGrad().

• spk_reset (str, optional) – Reset policy inherited from Neuron. 'hard' reproduces NEST hard reset behavior. Default is 'hard'.

• ref_var (bool, optional) – If True, allocates optional boolean state self.refractory for external refractory inspection. Default is False.

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

Parameter Mapping

Table 4 Parameter mapping to model symbols

Parameter	Type / shape / unit	Default	Math symbol	Semantics
in_size	Size; scalar or tuple	required	–	Defines population/state shape self.varshape.
E_L	ArrayLike, broadcastable to self.varshape (mV)	-70. * u.mV	\(E_L\)	Leak reversal (resting) potential.
C_m	ArrayLike, broadcastable (pF), > 0	250. * u.pF	\(C_m\)	Membrane capacitance in subthreshold integration.
tau_m	ArrayLike, broadcastable (ms), > 0	10. * u.ms	\(\tau_m\)	Membrane leak time constant.
t_ref	ArrayLike, broadcastable (ms), >= 0	2. * u.ms	\(t_\mathrm{ref}\)	Absolute refractory duration in physical time.
V_th and V_reset	ArrayLike, broadcastable (mV), with V_reset < V_th	-55. * u.mV, -70. * u.mV	\(V_\mathrm{th}\), \(V_\mathrm{reset}\)	Threshold and post-spike reset levels.
tau_syn	ArrayLike, flattened to (n_receptors,) (ms), each > 0	(2.0,) * u.ms	\(\tau_{\mathrm{syn},k}\)	Receptor-specific alpha time constants; length defines n_receptors.
I_e	ArrayLike, broadcastable (pA)	0. * u.pA	\(I_e\)	Constant current added each update step.
V_min	ArrayLike broadcastable (mV) or None	None	\(V_\mathrm{min}\)	Optional lower clamp on candidate membrane voltage.
V_initializer	Callable	Constant(-70. * u.mV)	–	Initializer for membrane state V.
spk_fun	Callable	ReluGrad()	–	Surrogate nonlinearity used for spike output.
spk_reset	str	'hard'	–	Reset mode from Neuron.
ref_var	bool	False	–	If True, exposes boolean state self.refractory.
name	str | None	None	–	Optional node name.

Raises:

• ValueError – Raised at initialization or update time if any of the following holds: - C_m <= 0, tau_m <= 0, any tau_syn <= 0, t_ref < 0, or V_reset >= V_th. - A spike event receptor index is outside [1, n_receptors].

• TypeError – If parameters or inputs are not unit-compatible with the expected conversions (mV, ms, pF, pA).

• KeyError – If simulation context entries (for example t or dt) are missing when update() is called.

• AttributeError – If update() is called before init_state() creates required state holders.

Notes

• State variables are V, y1_syn, y2_syn, i_const, refractory_step_count, and last_spike_time; refractory is added only when ref_var=True.

• Spike weights from spike_events and sum_delta_inputs are signed currents in pA: positive for depolarizing, negative for hyperpolarizing receptors. This differs from conductance-based multisynapse models where weights must be non-negative.

• update(x=...) stores x into i_const for use on the next step, matching NEST current-event buffering semantics.

• If n_receptors == 0, all spike event inputs are silently ignored and sum_delta_inputs is discarded.

• Default delta input from sum_delta_inputs is routed to receptor 1 when n_receptors > 0, replicating NEST default port behavior.

Examples

>>> import brainstate
>>> import saiunit as u
>>> from brainpy_state._nest.iaf_psc_alpha_multisynapse import (
...     iaf_psc_alpha_multisynapse,
... )
>>> with brainstate.environ.context(dt=0.1 * u.ms):
...     neu = iaf_psc_alpha_multisynapse(
...         in_size=3,
...         tau_syn=(1.5, 3.0) * u.ms,
...         I_e=180.0 * u.pA,
...     )
...     neu.init_state()
...     with brainstate.environ.context(t=0.0 * u.ms):
...         spk = neu.update(
...             spike_events=[{'receptor_type': 2, 'weight': 40.0 * u.pA}]
...         )
...     _ = spk.shape

>>> import brainstate
>>> import saiunit as u
>>> from brainpy_state._nest.iaf_psc_alpha_multisynapse import (
...     iaf_psc_alpha_multisynapse,
... )
>>> with brainstate.environ.context(dt=0.1 * u.ms):
...     neu = iaf_psc_alpha_multisynapse(in_size=1, tau_syn=(2.0,) * u.ms)
...     neu.init_state()
...     with brainstate.environ.context(t=0.0 * u.ms):
...         _ = neu.update(x=250.0 * u.pA)
...     with brainstate.environ.context(t=0.1 * u.ms):
...         spk_next = neu.update()
...     _ = spk_next

References

[1]

NEST source: models/iaf_psc_alpha_multisynapse.h and models/iaf_psc_alpha_multisynapse.cpp.

[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

[3]

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

[4]

Morrison A, Straube S, Plesser HE, Diesmann M (2007). Exact subthreshold integration with continuous spike times in discrete time neural network simulations. Neural Computation 19(1):47-79. DOI: https://doi.org/10.1162/neco.2007.19.1.47

get_spike(V=None)

Evaluate surrogate spike output for a voltage tensor.

Parameters:

V (ArrayLike or None, optional) – Voltage input in mV, broadcast-compatible with self.varshape. If None, uses current membrane state self.V.value.

Returns:

out – Surrogate spike output from self.spk_fun with the same shape as V (or self.V.value when V is None). The input to spk_fun is scaled as (V - V_th) / (V_th - V_reset) so the surrogate activates positively for suprathreshold voltages.

Return type:

dict

init_state(**kwargs)

Initialize runtime states for membrane, synaptic, and refractory variables.

Parameters:

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

Raises:

• ValueError – If initializers cannot broadcast to self.varshape.

• TypeError – If initializer outputs are incompatible with expected unit/array conversions for voltage, current, or integer refractory states.

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

Advance the neuron by one simulation step.

Parameters:

• x (ArrayLike, optional) – Continuous current input in pA for this step. x is accumulated through sum_current_inputs() and stored in i_const for use on the next call (one-step delayed buffering matching NEST ring-buffer semantics). Default is 0. * u.pA.

• spike_events (iterable or None, optional) –

  Receptor-indexed spike weight events to inject this step. Each entry must be either:

  • A (receptor_type, weight) tuple where receptor_type is a 1-based integer in [1, n_receptors] and weight is a scalar or array in pA (broadcastable to self.varshape).

  • A dict with keys 'receptor_type' (or 'receptor') and 'weight'.

  Multiple events for the same receptor are accumulated additively. None injects no receptor spike events. Default is None. Ignored when w_by_rec is provided.

• w_by_rec (array-like or None, optional) – Pre-computed per-receptor spike weights in pA (dimensionless), shape broadcastable to self.varshape + (n_receptors,). When provided, bypasses spike_events parsing and sum_delta_inputs, making the update JIT-compatible for use inside brainstate.transform.for_loop. Default is None.

Returns:

out – Spike output tensor from get_spike(), shape self.V.value.shape. On threshold crossings, the voltage presented to spk_fun is nudged above threshold by 1e-12 mV-equivalent to preserve positive surrogate activation.

Return type:

jax.Array

Raises:

• ValueError – If any receptor index in spike_events is outside [1, n_receptors].

• KeyError – If simulation context does not provide t or dt.

• AttributeError – If required states are missing because init_state() was not called.

• TypeError – If x or stored states are not unit-compatible with expected pA / mV conversions.