Simulators and fidelity

A simulator turns equations into numbers one step at a time, and every step involves choices, about where an input lands, what is held constant and how the arithmetic rounds. This chapter shows those choices, how far each one moves the answer, and how to tell whether two simulators agree.

Equations do not run themselves

Every model in this course is a differential equation, and no computer solves one directly. A simulator replaces it with a rule that moves the state forward by a step $\mathrm{d}t$, and applies the rule again and again. Two simulators given the same equations can disagree for three reasons. They can order the events of a step differently: when an input lands, when the threshold is tested, how long the reset is held. They can integrate differently, holding different things constant over a step. And they can round differently, in different precision or in a different order of operations. Chapter 9 showed what chaos does with the smallest difference: in a large network, any of the three becomes different spikes.

What happens in a step

sparx follows NEST’s order. In each step, the inputs due are delivered, every membrane is integrated, membranes at threshold fire and reset, refractory neurons stay held, and the new spikes are sent along their connections, to land after their delays. A spike sent in step $m$ over a delay of $D$ steps lands at the end of step $m + D$, and a synapse it lands on shapes the membrane from the next step on.

Simulators disagree about the details, and each difference moves spikes:

• A delta synapse, a jump of the voltage, lands in NEST before the threshold test, so it can fire the neuron in the step it arrives. Brian2’s on_pre="v += w" lands it after the test: it takes effect from the next step, decaying over it first, and is lost if the neuron fired. sparx does either: Delta() is NEST’s, Delta(after_threshold=True) Brian2’s.
• Brian2 counts the refractory period from the start of the step in which the neuron crossed threshold, so it holds the neuron one step less than NEST. sparx follows NEST, which never holds it for less than the stated period, and its parity test with Brian2 shortens t_ref by one step.
• snnTorch’s LIF by default subtracts the threshold one step after the spike, without decaying it. That is not a discretization of the continuous LIF, whose reset happens at the spike and then decays, so sparx resets in the spiking step and pins the difference in a test.

None of these is a bug in anyone’s code. Each is a choice, and a model’s results depend on it.

How accurate a step is

For a linear membrane with current input, chapter 2 found the step’s exact solution, so the step size changes nothing at the instants the steps share. NEST’s iaf_psc_exp integrates the same way, following Rotter and Diesmann, and sparx agrees with it to $10^{-13}$ mV.

Conductances break this. Within a step a synaptic conductance decays, so the membrane’s time constant changes during the step, and the equation has no simple exact solution. NEST integrates it adaptively, in as many small steps as its error bound needs. Brian2’s exponential_euler and sparx both hold the conductance constant over the step and solve the membrane exactly for that value, and differ in the value they hold: Brian2 holds the conductance at the start of the step, sparx, by default, its exact average over the step. The difference is the order of the scheme: the error of the first falls in proportion to $\mathrm{d}t$, that of the second in proportion to $\mathrm{d}t^2$.

Orange holds each conductance at its value at the start of the step, Brian2's exponential_euler: from 7.6e-1 mV at 1/2 ms to 9.2e-2 mV at 1/16 ms. Blue holds it at its exact mean over the step, sparx's default: from 7.9e-3 mV to 1.2e-4 mV. The truth is the mean scheme at 1/1024 ms.

A membrane receives AMPA input every 3 ms and GABA-A input every 7 ms, landing on millisecond edges so that every step size sees the same inputs. sparx runs it at four step sizes with each hold, against the mean scheme at 1/1024 ms:

from itertools import pairwise

import jax
import jax.numpy as jnp
import numpy as np

import sparx
from sparx.dynamics import (
    Arrivals,
    Exponential,
    LeakyIntegrateAndFire,
    PointNeuron,
    Receptor,
)

jax.config.update("jax_enable_x64", True)
ms = np.arange(1, 201)                           # 200 ms
weights = {"ampa": 6.0 * (ms % 3 == 0),          # nS, every 3 ms
           "gaba_a": 20.0 * (ms % 7 == 0)}       # every 7 ms


def membrane(dt, hold):
    """The membrane at the end of each ms; spikes land on ms edges."""
    cell = PointNeuron(
        LeakyIntegrateAndFire(v_th=jnp.inf),
        {"ampa": Receptor(Exponential(2.0), "conductance"),
         "gaba_a": Receptor(Exponential(5.0), "conductance")},
        hold=hold)
    per = round(1 / dt)
    spikes = {}
    for name, w in weights.items():
        spikes[name] = np.zeros(200 * per)
        spikes[name][per - 1::per] = w
    (_, v), _ = sparx.run(cell, Arrivals(0.0, spikes), dt=dt,
                          record=lambda s: s.neuron.v)
    return np.asarray(v)[per - 1::per]


truth = membrane(1 / 1024, "mean")
for hold in ("start", "mean"):
    errors = [np.abs(membrane(dt, hold) - truth).max()
              for dt in (1 / 2, 1 / 4, 1 / 8, 1 / 16)]
    ratios = [a / b for a, b in pairwise(errors)]
    print(f"hold={hold}:", ", ".join(f"{e:.1e}" for e in errors),
          "mV; each halving divides it by",
          ", ".join(f"{r:.1f}" for r in ratios))

Run it in ColabSource

hold=start: 7.6e-01, 3.7e-01, 1.8e-01, 9.2e-02 mV; each halving divides it by 2.0, 2.0, 2.0
hold=mean: 7.9e-03, 2.0e-03, 4.9e-04, 1.2e-04 mV; each halving divides it by 4.0, 4.0, 4.0

At half a millisecond the start-of-step hold is off by 0.76 mV, because a conductance that has just jumped is held at its peak for the whole step. The mean is a hundred times closer, and each halving of the step makes it four times closer, against two for the start.

Nonlinear neurons are harder again. AdEx’s exponential upswing is stiff, so sparx integrates it by RK4 in substeps of 0.01 ms. Hodgkin and Huxley’s gates are stiff under hyperpolarization: at −128 mV the sodium gate’s closing rate $\beta_m$ is about 130 per ms, and RK4 diverges there; a shorter substep only lowers the voltage at which it does. sparx’s default for it, Strang splitting, solves each gate with the voltage held and the voltage with the gates held, each exactly, which keeps every substep stable.

Three kinds of agreement

So when two simulators are compared, the right claim depends on what they compute. Where they integrate the same equations the same way in the same precision, they should agree to rounding. Where they integrate differently, the claim is a measured tolerance, with the reason for it. Where the system is chaotic, the last bit grows into different spikes, and the claim is statistics over many runs.

Below, sparx runs the inputs of the reference fixtures its parity tests run: 300 ms of random excitatory and inhibitory spikes and a constant current, in float64. Orange is the reference, blue is sparx drawn over it, and the dots below are their difference on a log scale.

recorded runs

LIF · NEST iaf_psc_expLIF · NEST iaf_cond_expLIF · Brian2 exponential_eulerIzhikevich · NEST izhikevich

Zoom

1 ×

LIF, exponential current synapses against NEST iaf_psc_exp: same numbers the same 10 spikes; voltages within 9.9e-14 mV.

Both integrate the membrane and the synaptic currents exactly over each step. Test: test_current_synapses_match_nest_to_rounding.

The current-based LIF agrees with NEST to $10^{-13}$ mV, and sparx’s conductance scheme with Brian2’s to $7 \times 10^{-14}$ mV once it holds the conductance as Brian2 does. Against NEST’s adaptive integration of the same conductance model, which lands within about $10^{-9}$ mV of the truth, sparx’s mean hold is within $7 \times 10^{-4}$ mV and fires the same 33 spikes. The Izhikevich neuron is NEST’s to the last bit, but only run operation by operation: compiled, XLA fuses the quadratic’s arithmetic and rounds its last bit differently, and the membrane amplifies that until a spike moves by a step.

Izhikevich’s twenty patterns

Izhikevich’s 2004 paper shows twenty behaviours of his model, from tonic spiking to inhibition-induced bursting, made by his MATLAB script, each panel with its own current, step and starting state. sparx runs each panel as the script does, beside the script itself run in GNU Octave.

1. A. Tonic spiking5 spikes on his steps · 3e-12 mV
2. B. Phasic spiking1 spike on his steps · 6e-11 mV
3. C. Tonic bursting28 spikes on his steps · 2e-11 mV
4. D. Phasic bursting6 spikes on his steps · 2e-10 mV
5. E. Mixed mode6 spikes on his steps · 5e-12 mV
6. F. Spike frequency adaptation6 spikes on his steps · 8e-13 mV
7. G. Class 1 excitable10 spikes on his steps · 2e-11 mV
8. H. Class 2 excitable14 spikes on his steps · 1e-2 mV
9. I. Spike latency1 spike on his steps · 7e-10 mV
10. J. Subthreshold oscillations1 spike on his steps · 6e-12 mV
11. K. Resonator1 spike on his steps · 5e-11 mV
12. L. Integrator1 spike on his steps · 7e-14 mV
13. M. Rebound spike1 spike on his steps · 1e-9 mV
14. N. Rebound burst7 spikes on his steps · 1e-9 mV
15. O. Threshold variability1 spike on his steps · 6e-12 mV
16. P. Bistability5 spikes on his steps · 3e-9 mV
17. Q. Depolarizing after potential1 spike on his steps · 4e-13 mV
18. R. Accommodation1 spike on his steps · 3e-13 mV
19. S. Inhibition induced spiking3 spikes on his steps · 3e-9 mV
20. T. Inhibition induced bursting12 spikes on his steps · 3e-11 mV

Blue is sparx's izhikevich_2004(pattern), orange his figure1.m run in GNU Octave, each with its own current, step and start. The voltages differ by at most the figure under each panel, between spikes. In class 2 excitability, Octave's V^2 rounds differently from v * v in about one value in a thousand, and the slow ramp through the bifurcation grows that to 0.012 mV.

Chaos and precision

For chaotic networks, spikes cannot be the test. sparx checks Brunel’s network at 2,500 neurons in each of his four regimes by the statistics chapter 9 used: one sparx run’s excitatory rate, mean CV and Fano factor must fall within NEST’s spread over eight seeds. The cortical microcircuit is checked both ways: over the reference’s own random draws at 15 seeds, each population’s mean rate is within four standard deviations of NEST’s, and on the one network sparx draws, NEST and sparx fire the same 12,689 spikes for 300 ms.

Precision is part of the comparison. sparx runs in float32 unless JAX’s 64-bit mode is on, and its parity checks run in float64. A float32 run of a chaotic network is a different run, not a less accurate copy, since chaos grows its rounding differences into different spikes, and it can only be compared by its statistics.

Speed

On a 4-core CPU, with NEST and Brian2 using four threads and sparx running in float32, a simulated second takes:

Network	Neurons, synapses	sparx	NEST	Brian2 standalone
Brunel	12,500, 15.6 million	9.6 s	7.5 s	11.8 s
CUBA	4,000, 320,000	0.71 s	0.42 s	0.33 s
COBA	4,000, 320,000	1.01 s	3.85 s	0.54 s
Microcircuit, a fifth	15,435, 12 million	8.76 s	2.89 s

NEST integrates COBA’s conductances adaptively, which is more accurate and costs it that row; sparx and Brian2 hold the conductance over the step. No GPU or TPU numbers exist yet. The fidelity ledger lists every model’s checks and every known difference, and performance the measurements.

Try this

1. In the parity figure, choose the conductance case against NEST and zoom in. Where does the difference grow, and where does it shrink?
2. From the chart, which step does each scheme need for the membrane to be within 0.001 mV of the truth?
3. A neuron’s membrane is 0.5 mV below threshold, and a delta input of 1 mV arrives at the end of the step. When does it fire in NEST, and in Brian2?
4. Why is a float32 run of Brunel’s network not a less accurate version of a float64 run?

Answers

1. It is exactly zero for 2 ms after each spike, while both neurons are held at the same reset voltage. Then it grows again as the membrane integrates its inputs, to about 5 × 10⁻⁴ mV before the next spike.
2. The mean hold is at $2.0 \times 10^{-3}$ mV at 1/4 ms and $4.9 \times 10^{-4}$ mV at 1/8 ms, so 1/8 ms does it. The start-of-step hold halves its error with each halving of the step, from 0.092 mV at 1/16 ms, so it needs about a hundred times shorter, under a thousandth of a millisecond.
3. In NEST in that step: the jump lands before the threshold test, and the membrane is 0.5 mV above threshold. In Brian2 the test comes first, so the jump lands after it, and the neuron fires in the next step if the membrane, decaying over it, is still above threshold.
4. A chaotic network grows any difference into different spikes, and float32’s rounding is a difference at every step. The float32 run is a different trajectory of the same network, as valid as the float64 one, and the two can only be compared by their statistics.

Summary

A simulator chooses where inputs land in a step, what it holds constant over the step and how it rounds, and each choice moves spikes. For linear membranes sparx integrates exactly and agrees with NEST to rounding; for conductances its mean hold is second order, a hundred times closer than the start-of-step hold at half a millisecond; for chaotic networks the comparison is statistics. The next chapter leaves simulation for hardware: event cameras, neuromorphic chips, and the format that moves a network between them.

References

• S. Rotter and M. Diesmann, “Exact digital simulation of time-invariant linear systems with applications to neuronal modeling”, Biological Cybernetics 81, 1999, doi:10.1007/s004220050570.
• M.-O. Gewaltig and M. Diesmann, “NEST (NEural Simulation Tool)”, Scholarpedia 2, 2007, doi:10.4249/scholarpedia.1430.
• M. Stimberg, R. Brette and D. F. M. Goodman, “Brian 2, an intuitive and efficient neural simulator”, eLife 8, 2019, doi:10.7554/eLife.47314.
• R. Brette et al., “Simulation of networks of spiking neurons: a review of tools and strategies”, Journal of Computational Neuroscience 23, 2007, doi:10.1007/s10827-007-0038-6.
• S. Rush and H. Larsen, “A practical algorithm for solving dynamic membrane equations”, IEEE Transactions on Biomedical Engineering 25, 1978, doi:10.1109/TBME.1978.326270.
• G. Strang, “On the construction and comparison of difference schemes”, SIAM Journal on Numerical Analysis 5, 1968, doi:10.1137/0705041.
• E. M. Izhikevich, “Which model to use for cortical spiking neurons?”, IEEE Transactions on Neural Networks 15, 2004, doi:10.1109/TNN.2004.832719.
• N. Brunel, “Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons”, Journal of Computational Neuroscience 8, 2000, doi:10.1023/A:1008925309027.
• T. C. Potjans and M. Diesmann, “The cell-type specific cortical microcircuit: relating structure and activity in a full-scale spiking network model”, Cerebral Cortex 24, 2014, doi:10.1093/cercor/bhs358.