pp_psc_delta

class brainpy.state.pp_psc_delta(in_size, tau_m=Quantity(10., "ms"), C_m=Quantity(250., "pF"), dead_time=1.0, dead_time_random=False, dead_time_shape=1, with_reset=True, tau_sfa=(), q_sfa=(), c_1=0.0, c_2=1.238, c_3=0.25, I_e=Quantity(0., "pA"), t_ref_remaining=0.0, rng_key=None, V_initializer=Constant(value=0. mV), spk_fun=ReluGrad(alpha=0.3, width=1.0), spk_reset='hard', name=None)

Point process neuron with leaky integration of delta-shaped PSCs.

pp_psc_delta is an implementation of a leaky integrator where the potential jumps on each spike arrival. It produces spikes stochastically according to a transfer function operating on the membrane potential, and supports spike-frequency adaptation with multiple exponential kernels.

This is a brainpy.state re-implementation of the NEST simulator model of the same name, using NEST-standard parameterization and exact integration.

Parameters:

• in_size (int, tuple of int) – Population shape. Defines the number of neurons in the population.

• tau_m (Quantity, optional) – Membrane time constant. Must be a positive quantity with time units. Default: 10.0 ms.

• C_m (Quantity, optional) – Membrane capacitance. Must be a positive quantity with capacitance units. Default: 250.0 pF.

• dead_time (float, optional) – Duration of the dead time (absolute refractory period) in milliseconds. If set to 0, the model operates in Poisson mode with potentially multiple spikes per time step. If dead_time is nonzero but smaller than the simulation resolution, it is clamped to the resolution. Must be non-negative. Default: 1.0 ms.

• dead_time_random (bool, optional) – Whether to draw random dead time after each spike from a gamma distribution. If True, dead_time becomes the mean of the gamma distribution with shape parameter dead_time_shape. Default: False.

• dead_time_shape (int, optional) – Shape parameter of the gamma distribution for random dead times. Must be at least 1. Default: 1.

• with_reset (bool, optional) – Whether to reset the membrane potential to 0 after each spike. Default: True.

• tau_sfa (tuple of float, optional) – Adaptive threshold time constants in milliseconds. Each element defines the decay time constant of one adaptation kernel. Must be a sequence of positive values with the same length as q_sfa. Default: () (no adaptation).

• q_sfa (tuple of float, optional) – Adaptive threshold jump sizes in millivolts. Each element defines the increment added to the corresponding adaptation kernel on each spike. Must be a sequence with the same length as tau_sfa. Default: () (no adaptation).

• c_1 (float, optional) – Slope of the linear part of the transfer function in Hz/mV. Default: 0.0.

• c_2 (float, optional) – Prefactor of the exponential part of the transfer function in Hz. Can be used as an offset spike rate when c_3 = 0. Default: 1.238 Hz.

• c_3 (float, optional) – Coefficient of exponential nonlinearity in 1/mV. Must be non-negative. Set to 0 for purely linear transfer function. Default: 0.25 1/mV.

• I_e (Quantity, optional) – Constant external input current. Must be a quantity with current units. Default: 0.0 pA.

• t_ref_remaining (float, optional) – Remaining dead time at simulation start in milliseconds. Must be non-negative. Default: 0.0 ms.

• rng_key (jax.Array, optional) – JAX PRNG key for stochastic spike generation. If None, a default key is used. For reproducible results, provide an explicit key. Default: None.

• V_initializer (Callable, optional) – Initializer for the membrane potential (relative to resting potential). Default: Constant(0.0 * u.mV).

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

• spk_reset (str, optional) – Reset mode. Options: 'hard' (stop gradient), 'soft' (V -= V_th). Default: 'hard' (matches NEST behavior).

• name (str, optional) – Name of the neuron population. Default: None.

Raises:

• ValueError – If C_m <= 0 (capacitance must be strictly positive).

• ValueError – If tau_m <= 0 (membrane time constant must be strictly positive).

• ValueError – If dead_time < 0 (dead time must be non-negative).

• ValueError – If dead_time_shape < 1 (gamma shape parameter must be at least 1).

• ValueError – If t_ref_remaining < 0 (remaining refractory time must be non-negative).

• ValueError – If c_3 < 0 (exponential coefficient must be non-negative).

• ValueError – If any element of tau_sfa <= 0 (adaptation time constants must be positive).

• ValueError – If len(tau_sfa) != len(q_sfa) (adaptation parameter lists must match).

See also

iaf_psc_delta

Integrate-and-fire neuron with delta PSCs

gif_psc_exp

Generalized integrate-and-fire with exponential PSCs

Parameter Mapping

NEST Parameter	brainpy.state	Notes
tau_m	tau_m	Membrane time constant
C_m	C_m	Membrane capacitance
dead_time	dead_time	Refractory period duration
dead_time_random	dead_time_random	Enable random dead time
dead_time_shape	dead_time_shape	Gamma distribution shape parameter
with_reset	with_reset	Reset V_m after spike
tau_sfa	tau_sfa	Adaptation time constants (list)
q_sfa	q_sfa	Adaptation jump sizes (list)
c_1	c_1	Linear transfer function coefficient
c_2	c_2	Exponential transfer function prefactor
c_3	c_3	Exponential transfer function exponent
I_e	I_e	External input current
t_ref_remaining	t_ref_remaining	Initial refractory time
V_m	V.value	Membrane potential (relative to rest)
E_sfa	_q_val	Sum of all adaptation elements

1. Mathematical Model

1.1. Membrane Dynamics

The membrane potential \(V_\mathrm{m}\) (relative to resting potential) evolves according to a leaky integrator:

\[C_\mathrm{m} \frac{dV_\mathrm{m}}{dt} = -\frac{V_\mathrm{m}}{\tau_\mathrm{m}} + I_\mathrm{e} + I_\mathrm{syn}(t)\]

where:

• \(C_\mathrm{m}\) is the membrane capacitance

• \(\tau_\mathrm{m}\) is the membrane time constant

• \(I_\mathrm{e}\) is the constant external input current

• \(I_\mathrm{syn}(t)\) is the synaptic input current

The exact (analytic) integration over one time step \(h\) gives:

\[V_\mathrm{m}(t + h) = P_{33} \cdot V_\mathrm{m}(t) + P_{30} \cdot (I_0 + I_\mathrm{e}) + w_\mathrm{syn}\]

where:

• \(P_{33} = \exp(-h / \tau_\mathrm{m})\)

• \(P_{30} = \frac{\tau_\mathrm{m}}{C_\mathrm{m}}(1 - P_{33})\)

• \(I_0\) is the buffered current from the previous step (ring buffer)

• \(w_\mathrm{syn}\) is the sum of all incoming delta-shaped PSP jumps (in mV)

1.2. Transfer Function

The instantaneous firing rate is computed from the effective membrane potential \(V' = V_\mathrm{m} - E_\mathrm{sfa}\) using a flexible transfer function:

\[\text{rate}(t) = \text{Rect}\!\left[ c_1 \cdot V'(t) + c_2 \cdot \exp(c_3 \cdot V'(t)) \right]\]

where \(\text{Rect}(x) = \max(0, x)\) ensures non-negative rates.

By adjusting c_1, c_2, and c_3, the transfer function can be:

• Linear: Set c_3 = 0, c_1 > 0 – \(\text{rate} = c_1 V' + c_2\)

• Exponential: Set c_1 = 0 – \(\text{rate} = c_2 \exp(c_3 V')\)

• Mixed: All coefficients nonzero – linear + exponential

1.3. Spike-Frequency Adaptation

The adaptive threshold \(E_\mathrm{sfa}\) is the sum of multiple exponential kernels, each with its own time constant and jump size:

\[\tau_{\mathrm{sfa},i} \frac{dE_{\mathrm{sfa},i}}{dt} = -E_{\mathrm{sfa},i}\]

\[E_{\mathrm{sfa},i}(t) \to E_{\mathrm{sfa},i}(t) + q_{\mathrm{sfa},i} \quad \text{(on spike)}\]

\[E_\mathrm{sfa}(t) = \sum_{i=1}^{n} E_{\mathrm{sfa},i}(t)\]

The adaptation kernels decay exponentially with exact propagators:

\[E_{\mathrm{sfa},i}(t + h) = E_{\mathrm{sfa},i}(t) \exp(-h / \tau_{\mathrm{sfa},i})\]

1.4. Stochastic Spike Generation

• With dead time (dead_time > 0): At most one spike per time step. A uniform random number \(u \sim \mathcal{U}(0,1)\) is compared to the spike probability:

  \[P(\text{spike}) = 1 - \exp(-\text{rate} \cdot h \cdot 10^{-3})\]

  A spike is generated if \(u \le P(\text{spike})\).

• Without dead time (dead_time = 0): Multiple spikes per step are possible. The number of spikes is drawn from a Poisson distribution:

  \[n_{\text{spikes}} \sim \text{Poisson}(\text{rate} \cdot h \cdot 10^{-3})\]

The factor \(10^{-3}\) converts from Hz*ms to a dimensionless rate.

1.5. Dead Time (Refractory Period)

After each spike, the neuron enters a dead time during which it cannot spike:

• Fixed dead time: dead_time_random = False. The neuron is refractory for exactly dead_time milliseconds, converted to grid steps.

• Random dead time: dead_time_random = True. The dead time is drawn from a gamma distribution with shape dead_time_shape and mean dead_time.

If dead_time is nonzero but smaller than the simulation resolution \(h\), it is clamped to \(h\).

2. Numerical Integration and Update Order

The discrete-time update per simulation step follows this order:

1. Update membrane potential via exact propagator (including external current and synaptic delta inputs).

2. Decay adaptation elements and compute total \(E_\mathrm{sfa}\).

3. Spike check:

   • If not refractory: compute effective potential \(V' = V_\mathrm{m} - E_\mathrm{sfa}\), compute instantaneous rate, draw random number and potentially emit spike(s). If spike occurs:

     • Jump all adaptation elements by q_sfa

     • Optionally reset \(V_\mathrm{m}\) to 0 (if with_reset = True)

     • Set dead time counter

   • If refractory: decrement dead time counter

4. Buffer external current for the next step (ring buffer semantics).

3. Important Implementation Notes

• Relative membrane potential: The membrane potential \(V_\mathrm{m}\) is stored relative to the resting potential (resting potential = 0 mV). This differs from iaf_psc_delta, which uses absolute potentials.

• Stochastic reproducibility: Because spiking is stochastic (random number drawn each step), exact spike-time reproducibility requires matching the random number generator state. For deterministic testing, set rng_key explicitly.

• Dead time < dt clamping: If dead_time is nonzero but smaller than the simulation resolution, it is internally clamped to the resolution to match NEST behavior.

• Poisson mode performance: For non-refractory neurons (dead_time = 0), Poisson random draws are used, which are slower than uniform random draws. For typical firing rates (<1 spike/time_step), setting a small dead_time (e.g., 1e-8 ms) is faster and nearly equivalent.

4. State Variables

State Variable	Type	Description
V	HiddenState	Membrane potential (relative to rest)
refractory_step_count	ShortTermState	Remaining dead time grid steps
I_stim	ShortTermState	Buffered current applied in next step
last_spike_time	ShortTermState	Last spike time (for recording)
_q_elems	HiddenState	Adaptation kernel elements (internal)
_q_val	ShortTermState	Total \(E_\mathrm{sfa}\) (internal)
_rng_state	JAX PRNG key	Random number generator state (internal)

• Default parameter values match NEST C++ source for pp_psc_delta, which are based on Jolivet et al. (2006) [2].

• tau_sfa and q_sfa default to empty tuples (no adaptation). In NEST, the C++ defaults of tau_sfa=34.0 and q_sfa=0.0 are immediately cleared in the constructor, resulting in empty vectors.

• The recordable V_m in NEST corresponds to self.V.value in brainpy.state.

• The recordable E_sfa in NEST corresponds to self._q_val (the sum of all adaptation elements).

References

[1]

Cardanobile S, Rotter S (2010). Multiplicatively interacting point processes and applications to neural modeling. Journal of Computational Neuroscience 28(2):267-284. DOI: https://doi.org/10.1007/s10827-009-0204-0

[2]

Jolivet R, Rauch A, Luescher H-R, Gerstner W (2006). Predicting spike timing of neocortical pyramidal neurons by simple threshold models. Journal of Computational Neuroscience 21:35-49. DOI: https://doi.org/10.1007/s10827-006-7074-5

[3]

Pozzorini C, Naud R, Mensi S, Gerstner W (2013). Temporal whitening by power-law adaptation in neocortical neurons. Nature Neuroscience 16:942-948. DOI: https://doi.org/10.1038/nn.3431

[4]

Grytskyy D, Tetzlaff T, Diesmann M, Helias M (2013). A unified view on weakly correlated recurrent networks. Frontiers in Computational Neuroscience, 7:131. DOI: https://doi.org/10.3389/fncom.2013.00131

[5]

Deger M, Schwalger T, Naud R, Gerstner W (2014). Fluctuations and information filtering in coupled populations of spiking neurons with adaptation. Physical Review E 90:6, 062704. DOI: https://doi.org/10.1103/PhysRevE.90.062704

[6]

Gerstner W, Kistler WM, Naud R, Paninski L (2014). Neuronal Dynamics: From single neurons to networks and models of cognition. Cambridge University Press.

[7]

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

Examples

Basic usage with default parameters:

>>> import brainpy.state as bst
>>> import saiunit as u
>>> neurons = bst.pp_psc_delta(100)
>>> neurons.init_all_states()

Exponential transfer function (default):

>>> neurons = bst.pp_psc_delta(
...     100,
...     c_1=0.0,      # no linear part
...     c_2=1.238,    # exponential prefactor
...     c_3=0.25      # exponential coefficient
... )

Linear transfer function with offset:

>>> neurons = bst.pp_psc_delta(
...     100,
...     c_1=10.0,     # linear slope (Hz/mV)
...     c_2=5.0,      # offset rate (Hz)
...     c_3=0.0       # disable exponential
... )

With spike-frequency adaptation:

>>> neurons = bst.pp_psc_delta(
...     100,
...     tau_sfa=(100.0, 1000.0),  # two adaptation kernels
...     q_sfa=(5.0, 10.0)          # jump sizes in mV
... )

Poisson mode (no dead time):

>>> neurons = bst.pp_psc_delta(
...     100,
...     dead_time=0.0  # multiple spikes per step possible
... )

Random dead time:

>>> neurons = bst.pp_psc_delta(
...     100,
...     dead_time=2.0,           # mean dead time (ms)
...     dead_time_random=True,   # enable random dead time
...     dead_time_shape=2        # gamma distribution shape
... )

Reproducible stochastic behavior:

>>> import jax
>>> key = jax.random.PRNGKey(42)
>>> neurons = bst.pp_psc_delta(100, rng_key=key)

get_spike(V=None)

Compute surrogate gradient spike output for backpropagation.

This method is used for computing differentiable spike outputs during training. For a stochastic point process neuron, the true spike output is random and computed in update(). This method provides a surrogate gradient based on the membrane potential.

Parameters:

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

Returns:

spike – Differentiable spike output. Shape matches V.

Return type:

jax.Array

Notes

• This method is primarily used for gradient-based optimization.

• The surrogate gradient is computed by scaling the membrane potential and passing it through spk_fun (e.g., ReluGrad).

• The true stochastic spike output is computed in update() and is not directly differentiable.

init_state(batch_size=None, **kwargs)

Initialize all state variables.

Allocates and initializes membrane potential, spike times, refractory counters, buffered currents, adaptation kernels, and random number generator state.

Parameters:

• batch_size (int or None, optional) – If provided, states are created with shape (batch_size, *varshape) to support batched simulation. If None, states have shape varshape.

• **kwargs (dict, optional) – Additional keyword arguments (ignored).

Notes

• Membrane potential is initialized using V_initializer.

• Last spike time is initialized to -1e7 ms (sufficiently in the past).

• Refractory counter is initialized based on t_ref_remaining.

• Adaptation kernels (_q_elems) are initialized to zero.

• Random number generator state is initialized from rng_key or a default key.

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

Update neuron state for one simulation step.

Performs the complete update sequence: (1) updates membrane potential via exact propagator, (2) decays adaptation kernels, (3) computes instantaneous firing rate and stochastically generates spikes, (4) buffers input current for the next step.

Parameters:

x (Quantity, optional) – External current input (with current units). This input is added to the sum of all registered current inputs via projections. Default: 0.0 pA.

Returns:

spike – Binary spike output array. Shape: in_size. Values are 1.0 where spikes occurred, 0.0 otherwise. In Poisson mode (dead_time = 0), values can be integers > 1 representing multiple spikes per step.

Return type:

jax.Array

Notes

Update order per time step:

1. Membrane potential update: Apply exact propagator to update \(V_\mathrm{m}\) using buffered current from the previous step, constant external current, and delta-shaped synaptic inputs.

2. Adaptation decay: Decay all adaptation kernel elements using exponential propagators. Compute total \(E_\mathrm{sfa}\).

3. Spike generation:

   • If not refractory: compute effective potential \(V' = V_\mathrm{m} - E_\mathrm{sfa}\), compute instantaneous rate from transfer function, draw random number(s), and potentially emit spike(s).

   • If spike occurs: jump adaptation elements by q_sfa, optionally reset \(V_\mathrm{m}\) to 0, set dead time counter.

   • If refractory: decrement dead time counter.

4. Buffer input: Store external current input for the next step (ring buffer semantics, matching NEST).

Spike generation modes:

• With dead time (dead_time > 0): At most one spike per step. Uses uniform random numbers and spike probability.

• Without dead time (dead_time = 0): Poisson-distributed spikes. Multiple spikes per step are possible.

Failure modes:

• If C_m or tau_m contain invalid values (NaN, Inf), membrane potential update will fail silently (produces NaN).

• If c_3 * V' causes overflow in exp(), the exponential term will saturate to infinity. The rectifier ensures the rate remains non-negative.

• If random number generator state is corrupted, spike generation will produce undefined results.

Performance considerations:

• Poisson mode (dead_time = 0) is slower due to Poisson random draws.

• Setting a small dead_time (e.g., 1e-8 ms) uses faster uniform random numbers and is nearly equivalent for typical firing rates.

• Random dead time (dead_time_random = True) requires additional gamma distribution samples per spike.