iaf_tum_2000

iaf_tum_2000#

class brainpy.state.iaf_tum_2000(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_ex=Quantity(2., "ms"), tau_syn_in=Quantity(2., "ms"), I_e=Quantity(0., "pA"), rho=Quantity(0.01, "Hz"), delta=Quantity(0., "mV"), tau_fac=Quantity(1000., "ms"), tau_psc=Quantity(2., "ms"), tau_rec=Quantity(400., "ms"), U=0.5, x=0.0, y=0.0, u=0.0, 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_tum_2000 neuron model.

Description

iaf_tum_2000 is a leaky integrate-and-fire neuron with exponential postsynaptic currents and integrated Tsodyks-Markram short-term synaptic plasticity. The model extends iaf_psc_exp by maintaining presynaptic resource states x (readily-releasable pool), y (cleft/active fraction), and u (release probability), and emitting a per-spike spike_offset signal that encodes the jump in y at each spike event. This signal is used for receptor-1 coupling between iaf_tum_2000 neurons, enabling dynamic synaptic efficacy.

The implementation follows NEST models/iaf_tum_2000.{h,cpp} update ordering and event semantics exactly, including NEST-style buffered input handling and receptor-type routing.

1. Membrane and synaptic dynamics

Subthreshold voltage evolution follows the same equation as iaf_psc_exp:

\[\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 I_0 is the buffered current from the previous time step. Synaptic currents decay exponentially:

\[\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}}}.\]

Receptor-1 current input I_1 is filtered through the excitatory kernel:

\[I_{\mathrm{syn,ex}} \leftarrow I_{\mathrm{syn,ex}} + (1 - e^{-h/\tau_{\mathrm{syn,ex}}}) I_1,\]

where h = dt is the simulation time step.

2. Tsodyks-Markram short-term plasticity on spike

When a neuron emits a spike at time t_spike, the Tsodyks states are updated. Let t_last be the previous spike time (with NEST-compatible first-spike convention: t_last = 0 if the internal last-spike time is negative, indicating no prior spike), and h_ts = t_spike - t_last.

Define propagators:

\[\begin{split}P_{uu} = \begin{cases} 0, & \tau_{\mathrm{fac}}=0 \\ e^{-h_{ts}/\tau_{\mathrm{fac}}}, & \text{otherwise} \end{cases}, \quad P_{yy} = e^{-h_{ts}/\tau_{\mathrm{psc}}},\end{split}\]
\[P_{zz} = \mathrm{expm1}(-h_{ts}/\tau_{\mathrm{rec}}) = e^{-h_{ts}/\tau_{\mathrm{rec}}} - 1,\]
\[P_{xy} = \frac{P_{zz}\tau_{\mathrm{rec}} - (P_{yy}-1)\tau_{\mathrm{psc}}}{\tau_{\mathrm{psc}}-\tau_{\mathrm{rec}}}.\]

With \(z = 1 - x - y\) (inactive/recovered fraction), NEST performs state propagation in this exact order:

\[\begin{split}u &\leftarrow u P_{uu}, \\ x &\leftarrow x + P_{xy}y - P_{zz}z, \\ y &\leftarrow y P_{yy},\end{split}\]

followed by utilization jump and resource transfer:

\[\begin{split}u &\leftarrow u + U(1-u), \\ \Delta y &= u x, \\ x &\leftarrow x - \Delta y, \\ y &\leftarrow y + \Delta y.\end{split}\]

spike_offset is set to \(\Delta y\) on spike steps, zero otherwise.

3. NEST update ordering

Per time step, the model follows this precise sequence:

  1. Update membrane potential if not refractory (exact exponential propagator).

  2. Decay synaptic currents \(I_{\mathrm{syn,ex}}\) and \(I_{\mathrm{syn,in}}\).

  3. Add filtered receptor-1 current to \(I_{\mathrm{syn,ex}}\).

  4. Add arriving spike inputs (positive weights to excitatory, non-positive to inhibitory).

  5. Perform threshold test (deterministic or escape-noise), assign refractory and reset.

  6. On emitted spike, update Tsodyks states (using the order above) and set spike_offset.

  7. Buffer current inputs i_0 and i_1 for the next step.

4. Escape-noise threshold dynamics

Spike generation uses deterministic thresholding when \(\delta < 10^{-10}\): \(V_{\mathrm{rel}} \ge \theta\), where \(\theta = V_{th} - E_L\).

For \(\delta > 0\), the model uses exponential hazard:

\[\phi(V) = \rho \exp\left(\frac{V_{\mathrm{rel}} - \theta}{\delta}\right),\]

with step spike probability \(p=\phi(V)\,h\,10^{-3}\) (h in ms, \(\phi\) in 1/s). Stochastic decisions use numpy.random.random.

5. Event semantics and receptor routing

The update() method accepts spike_events as an iterable of event descriptors in one of these formats:

  • (receptor_type, weight)

  • (receptor_type, weight, offset)

  • (receptor_type, weight, offset, multiplicity)

  • (receptor_type, weight, offset, multiplicity, sender_model)

  • dict with keys receptor_type/receptor, weight, offset, multiplicity, sender_model

Receptors:

  • Receptor 0 (DEFAULT): regular spike input, effective weight is weight * multiplicity.

  • Receptor 1 (TSODYKS): Tsodyks-coupled input, effective weight is weight * multiplicity * offset, where offset is typically the spike_offset from the presynaptic iaf_tum_2000 neuron.

For receptor 1, the sender_model field must be "iaf_tum_2000" (default assumption if not provided); otherwise a ValueError is raised, mirroring NEST’s connection constraints.

6. Stability constraints and computational implications

  • Construction validates: V_reset < V_th, C_m > 0, tau_m > 0, tau_syn_ex > 0, tau_syn_in > 0, tau_psc > 0, tau_rec > 0, tau_fac >= 0, t_ref >= 0, rho >= 0, delta >= 0, x + y <= 1, u ∈ [0,1].

  • Tsodyks state propagation uses the same singularity-free logic as NEST to handle tau_psc == tau_rec or tau_fac == 0 cases gracefully.

  • Per-call cost is \(O(\prod \mathrm{varshape})\) with vectorized NumPy operations in float64 for coefficient evaluation.

  • Buffered current semantics match NEST ring-buffer timing: x and x_filtered supplied at step n are stored and consumed at step n+1.

Parameters:

  • in_size (Size) – Population shape specification. All per-neuron parameters are broadcast to self.varshape.

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

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

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

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

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

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

  • tau_syn_ex (ArrayLike, optional) – Excitatory synaptic decay constant \(\tau_{\mathrm{syn,ex}}\) in ms; broadcastable and strictly positive. Default is 2. * bu.ms.

  • tau_syn_in (ArrayLike, optional) – Inhibitory synaptic decay constant \(\tau_{\mathrm{syn,in}}\) in ms; broadcastable and strictly positive. Default is 2. * bu.ms.

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

  • rho (ArrayLike, optional) – Escape-noise base firing intensity \(\rho\) in 1/s; broadcastable and nonnegative. Used only in stochastic mode (delta > 0). Default is 0.01 / bu.second.

  • delta (ArrayLike, optional) – Escape-noise soft-threshold width \(\delta\) in mV; broadcastable and nonnegative. delta == 0 reproduces deterministic thresholding. Default is 0. * bu.mV.

  • tau_fac (ArrayLike, optional) – Facilitation time constant \(\tau_{\mathrm{fac}}\) in ms; broadcastable and nonnegative. tau_fac == 0 disables facilitation (\(P_{uu}=0\)). Default is 1000. * bu.ms.

  • tau_psc (ArrayLike, optional) – Tsodyks postsynaptic current time constant \(\tau_{\mathrm{psc}}\) in ms; broadcastable and strictly positive. Used in state propagators. Default is 2. * bu.ms.

  • tau_rec (ArrayLike, optional) – Resource recovery time constant \(\tau_{\mathrm{rec}}\) in ms; broadcastable and strictly positive. Default is 400. * bu.ms.

  • U (ArrayLike, optional) – Utilization increment factor \(U\) (dimensionless); broadcastable and must lie in [0, 1]. Represents the per-spike increase in release probability. Default is 0.5.

  • x (ArrayLike, optional) – Initial readily-releasable resource fraction (dimensionless); broadcastable. Must satisfy x + y <= 1 and x >= 0. Default is 0.0.

  • y (ArrayLike, optional) – Initial cleft/active fraction (dimensionless); broadcastable. Must satisfy x + y <= 1 and y >= 0. Default is 0.0.

  • u (ArrayLike, optional) – Initial release probability (dimensionless); broadcastable and must lie in [0, 1]. Default is 0.0.

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

  • spk_fun (Callable, optional) – Surrogate spike nonlinearity used by get_spike(). Default is braintools.surrogate.ReluGrad().

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

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

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

Parameter Mapping

Parameter

Type / shape / unit

Default

Math symbol

Semantics

in_size

Size; scalar/tuple

required

–

Defines neuron population shape self.varshape.

E_L

ArrayLike, broadcastable (mV)

-70. * bu.mV

\(E_L\)

Resting membrane potential.

C_m

ArrayLike, broadcastable (pF), > 0

250. * bu.pF

\(C_m\)

Membrane capacitance in voltage integration.

tau_m

ArrayLike, broadcastable (ms), > 0

10. * bu.ms

\(\tau_m\)

Membrane leak time constant.

t_ref

ArrayLike, broadcastable (ms), >= 0

2. * bu.ms

\(t_{\mathrm{ref}}\)

Absolute refractory duration.

V_th and V_reset

ArrayLike, broadcastable (mV), with V_reset < V_th

-55. * bu.mV, -70. * bu.mV

\(V_{th}\), \(V_{\mathrm{reset}}\)

Threshold and post-spike reset voltages.

tau_syn_ex and tau_syn_in

ArrayLike, broadcastable (ms), each > 0

2. * bu.ms

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

Exponential PSC decay constants.

I_e

ArrayLike, broadcastable (pA)

0. * bu.pA

\(I_e\)

Constant current injected every step.

rho and delta

ArrayLike, broadcastable; rho in 1/s and delta in mV, both >= 0

0.01 / bu.second, 0. * bu.mV

\(\rho\), \(\delta\)

Escape-noise hazard parameters.

tau_fac

ArrayLike, broadcastable (ms), >= 0

1000. * bu.ms

\(\tau_{\mathrm{fac}}\)

Facilitation decay time constant; 0 disables.

tau_psc

ArrayLike, broadcastable (ms), > 0

2. * bu.ms

\(\tau_{\mathrm{psc}}\)

Tsodyks PSC time constant.

tau_rec

ArrayLike, broadcastable (ms), > 0

400. * bu.ms

\(\tau_{\mathrm{rec}}\)

Resource recovery time constant.

U

ArrayLike, broadcastable (dimensionless), ∈ [0,1]

0.5

\(U\)

Utilization increment per spike.

x

ArrayLike, broadcastable (dimensionless), x+y <= 1

0.0

\(x\)

Initial readily-releasable fraction.

y

ArrayLike, broadcastable (dimensionless), x+y <= 1

0.0

\(y\)

Initial cleft/active fraction.

u

ArrayLike, broadcastable (dimensionless), ∈ [0,1]

0.0

\(u\)

Initial release probability.

V_initializer

Callable

Constant(-70. * bu.mV)

–

Initializer for membrane state V.

spk_fun

Callable

ReluGrad()

–

Surrogate function for output spikes.

spk_reset

str (typically 'hard')

'hard'

–

Reset behavior selection in base class.

ref_var

bool

False

–

Enables explicit boolean refractory state variable.

name

str or None

None

–

Optional instance name.
Raises:

ValueError – Raised at construction when any validated constraint is violated: V_reset >= V_th, nonpositive C_m/tau_m/synaptic time constants/Tsodyks time constants, negative tau_fac/t_ref/rho/delta, U not in [0,1], u not in [0,1], or x + y > 1.

Examples

>>> import brainstate
>>> import saiunit as bu
>>> from brainpy_state._nest.iaf_tum_2000 import iaf_tum_2000
>>> brainstate.environ.set(dt=0.1 * bu.ms, t=0.0 * bu.ms)
>>> neu = iaf_tum_2000(
...     in_size=(2,),
...     I_e=250. * bu.pA,
...     tau_fac=500. * bu.ms,
...     tau_rec=400. * bu.ms,
...     U=0.3
... )
>>> neu.init_state()
>>> out = neu.update(x=0. * bu.pA, x_filtered=0. * bu.pA)
>>> out.shape
(2,)

Notes

  • Shares the exact exponential propagator implementation with iaf_psc_exp (via propagator_exp() from _utils).

  • The Tsodyks update order matches NEST iaf_tum_2000.cpp exactly to ensure identical dynamics in network simulations.

  • Receptor-1 connections require both presynaptic and postsynaptic neurons to be iaf_tum_2000 models, enforced via runtime validation.

  • spike_offset can be recorded and monitored for debugging or analysis of dynamic synaptic efficacy.

  • The model is grid-based with one-step input buffering matching NEST’s ring-buffer semantics.

get_spike(V=None)[source]#

Compute surrogate spike output given membrane potential.

Applies the surrogate spike function (self.spk_fun) to a normalized voltage that ranges from 0 (at reset) to 1 (at threshold). This enables differentiable spike computation for gradient-based learning.

Parameters:

V (ArrayLike or None, optional) – Membrane potential in mV; broadcastable to self.varshape. If None, uses self.V.value. Default is None.

Returns:

out – Surrogate spike activation, shape matching the input V (or self.V.value). The output is typically in [0, 1] for sub-threshold voltages and close to 1 for supra-threshold voltages, depending on the surrogate function used.

Return type:

dict

Notes

Voltage normalization:

\[v_{\mathrm{scaled}} = \frac{V - V_{th}}{V_{th} - V_{\mathrm{reset}}}.\]

The surrogate function self.spk_fun (default braintools.surrogate.ReluGrad()) is then applied to v_scaled, providing a differentiable approximation of the Heaviside step function.

init_state(**kwargs)[source]#

Initialize all neuron state variables.

Creates and allocates state variables for membrane potential, synaptic currents, refractory counter, Tsodyks-Markram plasticity states, and buffered inputs. All states are allocated as brainstate.HiddenState or brainstate.ShortTermState with shape self.varshape.

Parameters:

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

Notes

State variables created:

  • VHiddenState (mV)

    Membrane potential, initialized via self.V_initializer.

  • i_syn_exShortTermState (pA)

    Excitatory synaptic current, initialized to zero.

  • i_syn_inShortTermState (pA)

    Inhibitory synaptic current, initialized to zero.

  • i_0ShortTermState (pA)

    Buffered receptor-0 current input, initialized to zero.

  • i_1ShortTermState (pA)

    Buffered receptor-1 current input, initialized to zero.

  • refractory_step_countShortTermState (int32)

    Remaining refractory steps, initialized to zero.

  • last_spike_timeShortTermState (ms)

    Time of last emitted spike, initialized to -1e7 * bu.ms (no prior spike).

  • xShortTermState (dimensionless)

    Readily-releasable resource fraction, initialized to self.x_init.

  • yShortTermState (dimensionless)

    Cleft/active fraction, initialized to self.y_init.

  • uShortTermState (dimensionless)

    Release probability, initialized to self.u_init.

  • spike_offsetShortTermState (dimensionless)

    Per-spike \(\Delta y\) signal for receptor-1 coupling, initialized to zero.

  • refractoryShortTermState (bool), optional

    Boolean refractory flag, allocated only if ref_var=True.

property receptor_types#

Return a dictionary of available receptor type labels.

Returns:

Mapping from receptor name (str) to receptor ID (int): {'DEFAULT': 0, 'TSODYKS': 1}.

Return type:

dict

property recordables#

Return a list of state variable names available for recording.

Returns:

State variable names: ['V_m', 'I_syn_ex', 'I_syn_in', 'x', 'y', 'u', 'spike_offset']. Note that the membrane potential is exposed as 'V_m' (matching NEST convention) but stored internally as self.V.

Return type:

list of str

update(x=Quantity(0., 'pA'), x_filtered=Quantity(0., 'pA'), spike_events=None, _w_ex_jnp=None, _w_in_jnp=None)[source]#

Advance the neuron state by one simulation step.

Performs a complete integration step following NEST iaf_tum_2000 update order: membrane propagation (if not refractory), synaptic current decay, filtered-current injection, spike input addition, threshold test, Tsodyks state update on spike emission, and input buffering for the next step.

Parameters:

  • x (ArrayLike, optional) – Current input in pA for receptor-0 (standard current port). Scalar or array broadcastable to self.varshape. The value is buffered and applied in the next step (NEST ring-buffer semantics). Default is 0. * bu.pA.

  • x_filtered (ArrayLike, optional) – Current input in pA for receptor-1. It is buffered to self.i_1 and injected through excitatory exponential filtering at the next update step via (1 - P11_ex) * i_1. Scalar or array broadcastable to self.varshape. Default is 0. * bu.pA.

  • spike_events (Iterable or None, optional) –

    Collection of spike event descriptors for direct spike input. Each event can be:

    • (receptor_type, weight)

    • (receptor_type, weight, offset)

    • (receptor_type, weight, offset, multiplicity)

    • (receptor_type, weight, offset, multiplicity, sender_model)

    • dict with keys receptor_type/receptor (int or str), weight (ArrayLike in pA), offset (float, default 1.0), multiplicity (float, default 1.0), sender_model (str, default "iaf_tum_2000").

    Receptor types:

    • 0 (or 'DEFAULT'): regular spike input, effective weight is weight * multiplicity.

    • 1 (or 'TSODYKS'): Tsodyks-coupled input, effective weight is weight * multiplicity * offset, where offset is typically the spike_offset from the presynaptic neuron.

    For receptor 1, sender_model must be "iaf_tum_2000"; otherwise a ValueError is raised.

    Positive effective weights route to excitatory channel, non-positive to inhibitory channel. Default is None (no events).

Returns:

out – Surrogate spike output returned by get_spike(). The output is elementwise over the neuron state shape (and batch axis, if initialized). For emitted spikes, the voltage argument to get_spike() is nudged above threshold by 1e-12 mV to preserve positive spike activation under hard reset.

Return type:

jax.Array

Raises:

ValueError – If provided inputs cannot be broadcast to the internal state shape, or if receptor-1 events have sender_model != "iaf_tum_2000".

Notes

Update order (following NEST ``iaf_tum_2000.cpp``):

  1. Membrane propagation: If not refractory, update \(V_{\mathrm{rel}}\) using exact exponential propagators (same as iaf_psc_exp).

  2. Synaptic decay: Exponentially decay i_syn_ex and i_syn_in.

  3. Filtered current injection: Add (1 - exp(-h/tau_syn_ex)) * i_1 to i_syn_ex.

  4. Spike input addition: Add arriving spike inputs (from spike_events and registered delta inputs) to i_syn_ex and i_syn_in by sign.

  5. Threshold test: Determine spike condition (deterministic or escape-noise), assign refractory counter, and reset voltage.

  6. Tsodyks update: On emitted spike, update (u, x, y) states in NEST order, compute \(\Delta y\), and set spike_offset.

  7. Buffer inputs: Store x and x_filtered for next step.

Tsodyks state update on spike:

When a spike is emitted, inter-spike interval h_ts = t_spike - t_last is computed (with t_last = 0 if last_spike_time < 0). Propagators are:

\[\begin{split}P_{uu} = \begin{cases} 0, & \tau_{\mathrm{fac}}=0 \\ e^{-h_{ts}/\tau_{\mathrm{fac}}}, & \text{otherwise} \end{cases}, \quad P_{yy} = e^{-h_{ts}/\tau_{\mathrm{psc}}}, \quad P_{zz} = e^{-h_{ts}/\tau_{\mathrm{rec}}} - 1,\end{split}\]
\[P_{xy} = \frac{P_{zz}\tau_{\mathrm{rec}} - (P_{yy}-1)\tau_{\mathrm{psc}}}{\tau_{\mathrm{psc}}-\tau_{\mathrm{rec}}}.\]

Then states update as:

\[\begin{split}u \leftarrow u P_{uu}, \quad x \leftarrow x + P_{xy}y - P_{zz}(1-x-y), \quad y \leftarrow y P_{yy}, \\ u \leftarrow u + U(1-u), \quad \Delta y = u x, \quad x \leftarrow x - \Delta y, \quad y \leftarrow y + \Delta y.\end{split}\]

spike_offset is set to \(\Delta y\) on spike, zero otherwise.

Performance: Per-step computational cost is \(O(\prod \mathrm{varshape})\) with vectorized NumPy operations in float64 for coefficient computation and state updates.