siegert_neuron

class brainpy.state.siegert_neuron(in_size, tau=Quantity(1., 'ms'), tau_m=Quantity(5., 'ms'), tau_syn=Quantity(0., 'ms'), t_ref=Quantity(2., 'ms'), mean=0.0, theta=15.0, V_reset=0.0, rate_initializer=Constant(value=0.0), name=None)

NEST-compatible siegert_neuron mean-field rate model.

1. Overview

Mean-field rate model using the Siegert gain function of a noisy LIF neuron. This model computes the population-averaged firing rate from drift-diffusion input statistics via an analytic transfer function, enabling efficient large-scale network simulation without explicit spike generation.

2. Mathematical Description

The rate dynamics follow a first-order ODE:

\[\tau\,\frac{dr(t)}{dt} = -r(t) + \text{mean} + \Phi(\mu, \sigma^2),\]

where:

• \(r(t)\) is the population firing rate (Hz)

• \(\tau\) is the rate time constant

• \(\text{mean}\) is a constant baseline drive

• \(\Phi(\mu, \sigma^2)\) is the Siegert transfer function

• \(\mu\) is the total drift input (mean membrane potential shift)

• \(\sigma^2\) is the total diffusion input (variance)

The Siegert function \(\Phi\) analytically computes the steady-state firing rate of a leaky integrate-and-fire neuron receiving white noise with drift \(\mu\) and diffusion \(\sigma^2\), subject to threshold \(\theta\), reset \(V_{\text{reset}}\), refractory period \(t_{\text{ref}}\), and membrane time constant \(\tau_m\) [2].

For colored noise (finite \(\tau_{\text{syn}} > 0\)), a threshold shift correction is applied [3]:

\[\Delta_{\text{th}} = \frac{\alpha}{2} \sqrt{\frac{\tau_{\text{syn}}}{\tau_m}},\]

where \(\alpha = |\zeta(1/2)| \sqrt{2} \approx 2.0653\).

The integration is performed via exact exponential propagators:

\[r(t + \Delta t) = e^{-\Delta t / \tau} r(t) + \left(1 - e^{-\Delta t / \tau}\right) \left(\text{mean} + \Phi(\mu, \sigma^2)\right).\]

3. NEST-Compatible Update Ordering (Non-WFR Path)

For each simulation step:

1. Collect delayed and instantaneous diffusion-event buffers from queues.

2. Sum all drift and diffusion contributions (delayed, instant, direct inputs).

3. Evaluate Siegert transfer function \(\Phi(\mu_{\text{total}}, \sigma^2_{\text{total}})\).

4. Update rate via exact exponential step: \(r \leftarrow P_1 r + P_2 (\text{mean} + \Phi)\).

5. Publish updated rate to delayed_rate and instant_rate buffers for outgoing connections.

This mirrors NEST’s non-waveform-relaxation update_ semantics where emitted diffusion coefficients are overwritten with the post-update rate.

4. Diffusion Event Handling

Runtime diffusion events modulate drift and diffusion inputs. Events can be supplied via two channels:

• instant_diffusion_events: applied in the current step (delay = 0)

• delayed_diffusion_events: scheduled by integer delay_steps (default 1)

Event format supports dicts, tuples, or lists. Dict keys:

• coeff (or rate/value): base coefficient

• drift_factor: multiplier for drift contribution

• diffusion_factor: multiplier for diffusion contribution

• weight: connection weight (default 1)

• multiplicity: event count (default 1)

• delay_steps (or delay): integer delay in steps

Tuple/list format: (coeff, drift_factor, diffusion_factor, delay_steps, weight, multiplicity). Shorter tuples use default values for trailing fields.

Drift and diffusion contributions are computed as:

\[\begin{split}\mu &= \text{coeff} \times \text{weight} \times \text{multiplicity} \times \text{drift\_factor}, \\ \sigma^2 &= \text{coeff} \times \text{weight} \times \text{multiplicity} \times \text{diffusion\_factor}.\end{split}\]

5. Siegert Transfer Function Computation

The Siegert function is evaluated element-wise for array inputs. For each population element, the computation handles three regimes:

• Deterministic (σ² ≤ 0): If μ > θ, returns LIF firing rate; else 0.

• Very subthreshold (θ - μ > 6σ): Returns 0 (Brunel 2000 fast path).

• General diffusive: Computes via integral of scaled complementary error function (erfcx) and Dawson’s integral, using either SciPy (if available) or custom Gauss-Legendre quadrature with asymptotic expansions.

Numerical integration uses 64-point Gauss-Legendre quadrature for erfcx and adaptive segmentation for Dawson’s integral, ensuring relative accuracy ~1.5e-8.

Parameters:

• in_size (Size) – Population shape. Tuple of ints or single int for 1D populations. Determines the spatial structure of the rate model. For example, (10, 10) creates a 10×10 grid of mean-field neurons.

• tau (Quantity[ms], optional) – Time constant of the first-order rate dynamics (must be > 0). Controls the rate of convergence to the steady-state Siegert value. Smaller values produce faster tracking of input changes. Default: 1 ms.

• tau_m (Quantity[ms], optional) – Membrane time constant used in the Siegert gain function (must be > 0). Represents the passive membrane time constant of the modeled LIF neurons. Default: 5 ms.

• tau_syn (Quantity[ms], optional) – Synaptic time constant for colored-noise threshold correction (must be ≥ 0). When tau_syn > 0, applies a threshold shift to account for finite synaptic rise time [3]. Use 0 ms for white noise (no correction). Default: 0 ms.

• t_ref (Quantity[ms], optional) – Refractory period in the Siegert gain function (must be ≥ 0). Represents the absolute refractory period during which the neuron cannot spike. Increases the interspike interval and reduces firing rates. Default: 2 ms.

• mean (float, optional) – Constant additive baseline drive in the rate ODE (dimensionless). Shifts the firing rate upward without affecting dynamics. Can be scalar or array matching in_size. Default: 0.0.

• theta (float, optional) – Spike threshold relative to resting potential (dimensionless, corresponds to mV in NEST). Must be > V_reset. Defines the firing threshold in the Siegert transfer function. Default: 15.0.

• V_reset (float, optional) – Reset potential relative to resting potential (dimensionless, corresponds to mV in NEST). Must be < theta. Neuron is reset to this value after spiking in the underlying LIF model. Default: 0.0.

• rate_initializer (Callable, optional) – Initializer function for the rate state variable. Called as rate_initializer(shape, batch_size) during init_state(). Default: braintools.init.Constant(0.0) (all neurons start at 0 Hz).

• name (str, optional) – Unique identifier for this module. If None, auto-generated. Used for logging and debugging.

Parameter Mapping

NEST Parameter	brainpy.state Parameter	Description
tau	tau	Rate dynamics time constant
tau_m	tau_m	Membrane time constant (Siegert)
tau_syn	tau_syn	Synaptic time constant (threshold shift)
t_ref	t_ref	Refractory period (Siegert)
mean	mean	Constant baseline drive
theta	theta	Spike threshold (relative to rest)
V_reset	V_reset	Reset potential (relative to rest)
rate	rate	Current firing rate (Hz)

Raises:

• ValueError – If tau ≤ 0, tau_m ≤ 0, tau_syn < 0, t_ref < 0, or V_reset ≥ theta.

• ValueError – If instant_diffusion_events contains non-zero delay_steps.

• ValueError – If delayed_diffusion_events contains negative delay_steps.

• ValueError – If event tuples have length > 6 or < 1.

Notes

Computational Complexity

• Siegert evaluation is the primary bottleneck (O(N) per neuron).

• Without SciPy, custom quadrature adds ~10× overhead.

• Delayed event queues are sparse dicts (O(1) insertion, O(K) retrieval for K active delays).

Numerical Stability:

• Uses erfcx(x) = exp(x²) erfc(x) to avoid overflow for large x.

• Asymptotic expansions for erfcx and Dawson’s integral when x > 8.

• Exact exponential propagators (exp and expm1) prevent drift accumulation.

Batch Dimensions:

States support an optional leading batch dimension for parallelizing multiple network realizations. Initialize with init_state(batch_size=B) to create shape (B, *in_size).

Integration with NEST:

This implementation reproduces NEST 3.9+ behavior for siegert_neuron in non-waveform-relaxation mode. Key differences:

• NEST uses precise spike times; brainpy.state uses fixed-step integration.

• NEST’s WFR mode (iterative delay resolution) is not implemented.

• Event formats are compatible but may differ in edge cases (consult NEST docs).

References

[1]

Hahne J, Dahmen D, Schuecker J, Frommer A, Bolten M, Helias M, Diesmann M (2017). Integration of continuous-time dynamics in a spiking neural network simulator. Frontiers in Neuroinformatics, 11:34. DOI: 10.3389/fninf.2017.00034.

[2] (1,2)

Fourcaud N, Brunel N (2002). Dynamics of the firing probability of noisy integrate-and-fire neurons. Neural Computation, 14(9):2057-2110. DOI: 10.1162/089976602320264015.

[3] (1,2,3)

Schuecker J, Diesmann M, Helias M (2015). Modulated escape from a metastable state driven by colored noise. Physical Review E, 92:052119. DOI: 10.1103/PhysRevE.92.052119.

Examples

Basic usage with constant input:

>>> import brainpy.state as bp
>>> import saiunit as u
>>> import brainstate
>>> model = bp.siegert_neuron(in_size=10, tau=2*u.ms, tau_m=10*u.ms)
>>> model.init_all_states()
>>> with brainstate.environ.context(dt=0.1*u.ms):
...     for _ in range(100):
...         rate = model.update(drift_input=12.0, diffusion_input=4.0)
>>> print(rate)  # Steady-state firing rate in Hz

Using diffusion events for network coupling:

>>> model.init_all_states()
>>> event = {'coeff': 50.0, 'drift_factor': 0.1, 'diffusion_factor': 0.05, 'delay_steps': 1}
>>> with brainstate.environ.context(dt=0.1*u.ms):
...     rate = model.update(delayed_diffusion_events=event)
>>> print(model.rate.value)  # Rate after delayed event delivery

Mean-field network with recurrent connections:

>>> exc = bp.siegert_neuron(in_size=800, tau=1*u.ms, theta=15.0)
>>> inh = bp.siegert_neuron(in_size=200, tau=1*u.ms, theta=15.0)
>>> exc.init_all_states()
>>> inh.init_all_states()
>>> # Simulate recurrent network (conceptual; requires projection setup)
>>> with brainstate.environ.context(dt=0.1*u.ms):
...     for t in range(1000):
...         exc_drive = exc.rate.value.sum() * 0.01
...         inh_drive = inh.rate.value.sum() * -0.02
...         exc.update(drift_input=exc_drive + inh_drive, diffusion_input=2.0)
...         inh.update(drift_input=exc_drive, diffusion_input=1.0)

init_state(**kwargs)

State initialization function.

siegert_rate(mu, sigma_square)

Evaluate the NEST-compatible Siegert transfer function.

Computes the steady-state firing rate \(\Phi(\mu, \sigma^2)\) of a noisy LIF neuron with drift mu and diffusion sigma_square, using the analytic Siegert formula [2] with optional colored-noise correction [3].

The computation is vectorized over population elements. Inputs are broadcast with model parameters (theta, tau_m, etc.) to produce an output array matching the broadcast shape.

Parameters:

• mu (ArrayLike) – Drift input (mean membrane potential shift, dimensionless). Scalar or array broadcastable with sigma_square and model parameters. Positive values depolarize the neuron. Typically in the range [0, 30] for physiological parameters.

• sigma_square (ArrayLike) – Diffusion input (membrane potential variance, dimensionless squared). Scalar or array broadcastable with mu and model parameters. Must be non-negative. Typical values: 0.1–10 for moderate noise. Zero produces deterministic LIF behavior.

Returns:

rate – Firing rate in Hz (shape matches broadcast of inputs and model parameters). Values ≥ 0. Returns 0 for subthreshold inputs (μ < θ with low noise). Maximum rate is approximately 1000 / t_ref Hz (refractory limit).

Return type:

ndarray

Notes

Special Cases:

• If sigma_square ≤ 0: deterministic LIF (returns 0 if μ ≤ θ, else fires at constant-input rate).

• If (θ - μ) > 6σ: deep subthreshold (returns 0, Brunel 2000 fast path).

• If t_ref = 0: no refractory limit (rate can diverge for μ >> θ).

Performance:

• Without SciPy: uses 64-point Gauss-Legendre quadrature (~10× slower).

• With SciPy: uses scipy.integrate.quad and scipy.special (faster).

Broadcasting Rules:

Output shape is np.broadcast(mu, sigma_square, theta).shape. For example, if model has in_size=(10,), mu is scalar, and sigma_square has shape (10,), output shape is (10,).

Examples

Single neuron with varying drift:

>>> import brainpy.state as bp
>>> import numpy as np
>>> import saiunit as u
>>> model = bp.siegert_neuron(in_size=1, tau_m=10*u.ms, t_ref=2*u.ms, theta=15.0)
>>> mu_vals = np.linspace(0, 25, 50)
>>> rates = model.siegert_rate(mu=mu_vals, sigma_square=2.0)
>>> print(rates.shape)  # (50,)
>>> print(rates.max())  # Maximum firing rate in Hz

Population with heterogeneous noise:

>>> model = bp.siegert_neuron(in_size=100, tau_m=10*u.ms)
>>> sigma_sq = np.linspace(0.1, 5.0, 100)
>>> rates = model.siegert_rate(mu=15.0, sigma_square=sigma_sq)
>>> print(rates.shape)  # (100,)

2D grid with spatially varying input:

>>> model = bp.siegert_neuron(in_size=(10, 10), tau_m=10*u.ms)
>>> mu_grid = np.random.uniform(10, 20, size=(10, 10))
>>> rates = model.siegert_rate(mu=mu_grid, sigma_square=3.0)
>>> print(rates.shape)  # (10, 10)

update(x=0.0, drift_input=0.0, diffusion_input=0.0, instant_diffusion_events=None, delayed_diffusion_events=None, _precomputed_drive=None)

Advance the rate dynamics by one simulation timestep.

Integrates the first-order rate ODE using exact exponential propagators, incorporating drift/diffusion inputs from multiple sources (direct inputs, current/delta hooks, and diffusion events). Updates internal state variables and publishes the new rate to delayed_rate and instant_rate buffers for outgoing connections.

Update Sequence:

1. Retrieve timestep dt from brainstate.environ.

2. Drain delayed event queues for the current step index.

3. Schedule incoming delayed events into future queue slots.

4. Accumulate instant events (must have delay_steps=0).

5. Sum all drift and diffusion contributions:

   • Delayed events (from queue)

   • Scheduled delayed events with delay_steps=0

   • Instant events

   • Direct inputs (drift_input, diffusion_input)

   • Dynamics hooks (current_inputs, delta_inputs)

6. Evaluate Siegert transfer function \(\Phi(\mu_{\text{total}}, \sigma^2_{\text{total}})\).

7. Update rate: \(r \leftarrow P_1 r + P_2 (\text{mean} + \Phi)\), where \(P_1 = e^{-\Delta t / \tau}\) and \(P_2 = 1 - P_1\).

8. Copy new rate to delayed_rate and instant_rate (NEST non-WFR semantics).

9. Increment internal step counter.

Parameters:

• x (ArrayLike, optional) – External input passed to sum_current_inputs() hook (dimensionless). Scalar or array broadcastable to in_size. Used for compatibility with standard Dynamics input API. Default: 0.0.

• drift_input (ArrayLike, optional) – Direct drift contribution (dimensionless). Scalar or array broadcastable to in_size. Added to total drift before Siegert evaluation. Positive values increase firing rate. Default: 0.0.

• diffusion_input (ArrayLike, optional) – Direct diffusion contribution (dimensionless squared). Scalar or array broadcastable to in_size. Added to total diffusion (variance) before Siegert evaluation. Must be non-negative. Default: 0.0.

• instant_diffusion_events (None, dict, tuple, list, optional) –

  Diffusion events applied in the current step (delay = 0). Can be:

  • None: no events

  • Single dict: {'coeff': float, 'drift_factor': float, ...}

  • Tuple/list of event dicts

  • Tuple of (coeff, drift_factor, diffusion_factor, …)

  All events must have delay_steps=0 (implicit or explicit). Raises ValueError if non-zero delay is specified. Default: None.

• delayed_diffusion_events (None, dict, tuple, list, optional) – Diffusion events scheduled for future delivery. Format identical to instant_diffusion_events, but delay_steps can be any non-negative integer (default 1). Events with delay_steps=0 are applied immediately. Negative delays raise ValueError. Default: None.

Returns:

rate – Updated firing rate in Hz (shape matches in_size or (batch_size, *in_size)). Also stored in self.rate.value. Values are non-negative.

Return type:

ndarray

Raises:

• ValueError – If instant_diffusion_events contains events with delay_steps != 0.

• ValueError – If delayed_diffusion_events contains events with delay_steps < 0.

• ValueError – If event tuples have invalid length (must be 1–6 elements).

Notes

State Updates:

The following state variables are modified in-place:

• self.rate: current firing rate (Hz)

• self.delayed_rate: rate for delayed connections (copy of rate)

• self.instant_rate: rate for instant connections (copy of rate)

• self._step_count: internal step counter (int64)

Event queues (_delayed_drift_queue, _delayed_diffusion_queue) are updated: delivered events are removed, new events are added.

Numerical Properties:

• Exact integration: exponential propagators ensure no drift accumulation.

• Stability: unconditionally stable for all tau > 0 and dt > 0.

• Precision: limited by Siegert evaluation accuracy (~1.5e-8 relative error).

Broadcasting:

All inputs are broadcast to a common state_shape, which is the maximum of self.rate.value.shape and any batch dimension. Scalar inputs are automatically tiled.

NEST Compatibility:

Reproduces NEST’s non-waveform-relaxation update semantics:

• Delayed events use integer step delays (not continuous time).

• Outgoing diffusion coefficients are updated post-integration (not mid-step).

• No iterative waveform relaxation (NEST’s WFR mode is not implemented).

Examples

Single step with constant input:

>>> import brainpy.state as bp
>>> import saiunit as u
>>> import brainstate
>>> model = bp.siegert_neuron(in_size=10, tau=2*u.ms)
>>> model.init_all_states()
>>> with brainstate.environ.context(dt=0.1*u.ms):
...     rate = model.update(drift_input=12.0, diffusion_input=3.0)
>>> print(rate.shape)  # (10,)
>>> print(model.rate.value)  # Updated firing rates

Using delayed events for network coupling:

>>> model.init_all_states()
>>> event = {'coeff': 50.0, 'drift_factor': 0.1, 'diffusion_factor': 0.05, 'delay_steps': 5}
>>> with brainstate.environ.context(dt=0.1*u.ms):
...     for step in range(10):
...         rate = model.update(delayed_diffusion_events=event if step == 0 else None)
...         if step == 5:
...             print(f"Event delivered at step {step}, rate = {rate[0]:.2f} Hz")

Batch simulation with heterogeneous parameters:

>>> model = bp.siegert_neuron(in_size=100, tau=1*u.ms)
>>> model.init_all_states(batch_size=32)  # 32 independent realizations
>>> drift = np.random.uniform(10, 20, size=(32, 100))
>>> with brainstate.environ.context(dt=0.1*u.ms):
...     rate = model.update(drift_input=drift, diffusion_input=2.0)
>>> print(rate.shape)  # (32, 100)

Multiple simultaneous instant events:

>>> events = [
...     {'coeff': 10.0, 'drift_factor': 1.0, 'diffusion_factor': 0.0},
...     {'coeff': 5.0, 'drift_factor': 0.5, 'diffusion_factor': 0.1}
... ]
>>> model.init_all_states()
>>> with brainstate.environ.context(dt=0.1*u.ms):
...     rate = model.update(instant_diffusion_events=events)
>>> print(f"Combined event effect: {rate.mean():.2f} Hz")