siegert_neuron

class brainpy.state.siegert_neuron(*args, **kwargs)

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:

>>> from brainpy import state as bp
>>> import brainunit 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:

>>> from brainpy import state as bp
>>> import numpy as np
>>> import brainunit 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)

Advance the Siegert rate by one step (NEST non-WFR semantics).

Drift and diffusion are read from the dual-channel substrate seam that a diffusion_connection deposits into (goal 15c, design A):

• drift \(\mu = \mathrm{sum\_current\_inputs}(x, r) + \mathrm{drift\_input} + \mathrm{sum\_delta\_inputs}(0,\ \text{label}=\) 'diffusion_mu' \()\),

• diffusion \(\sigma^2 = \mathrm{diffusion\_input} + \mathrm{sum\_delta\_inputs}(0,\ \text{label}=\) 'diffusion_sigma2' \()\).

The two channels carry distinct labels so a single diffusion_connection making two seam deposits (drift_factor * rate and diffusion_factor * rate) never cross-contaminates \(\mu\) and \(\sigma^2\). The rate then relaxes by the exact exponential propagator \(r \leftarrow P_1 r + P_2(\mathrm{mean} + \Phi)\).

The Siegert transfer is evaluated with the JAX port _siegert_phi_jax(), so the whole step lowers under brainstate.transform.for_loop / jit (drive it with those, not a bare Python loop).

Parameters:

• x (ArrayLike, optional) – External drive forwarded to sum_current_inputs (dimensionless).

• drift_input (ArrayLike, optional) – Direct drift contribution added to the total drift before the Siegert evaluation.

• diffusion_input (ArrayLike, optional) – Direct diffusion (variance) contribution; must be non-negative.

Returns:

rate_new – Updated firing rate in Hz (shape self.rate.value.shape).

Return type:

jax.Array