Fast and Flexible Spiking Network Simulations on Massively Distributed On-Chip Memory Accelerators

Fast and Flexible Spiking Network Simulations on Massively Distributed On-Chip Memory Accelerators

arXiv:2610.18217v1 [cs.PF] 13 Oct 2026

Catherine M. Schöfmann1,2*, Jan Finkbeiner1,2 , Jan Vogelsang1,2 , Susanne Kunkel1 1

Neuromorphic Software Ecosystems (PGI-15), Jülich Research Centre, Jülich, Germany. 2 RWTH Aachen University, Aachen, Germany.

*Corresponding author(s). E-mail(s): c.schoefmann@fz-juelich.de; Abstract Studying brain function and dysfunction in spiking networks with realistic synaptic densities pushes simulation technology beyond its current limits. It requires simulations substantially faster than biological real-time to cover long model times and enable efficient iterative exploration, particularly when connectivity is plastic. The modelling community aspires to a 100× speedup. Yet on conventional many-core systems this remains elusive, bounded by inter-node communication and memory access. So far, significant speedup has only been shown on a custombuilt FPGA system. We demonstrate the potential of commercial systems with large distributed on-chip memory to accelerate brain research. To exploit the local SRAM, we introduce a group-based mapping and communication strategy tackling spike-routing challenges on massively distributed systems. On Graphcore IPUs, the cortical microcircuit runs at a record 22.7× faster than biological real-time with negligible build times and a reusable compiled graph. By design, the algorithm makes no model- or hardware-specific assumptions, providing a general strategy for advancing neuronal simulation technology.

1

Fig. 1 Group-based algorithm. Centre: the cores (small grey squares) of all compute nodes (large grey squares) form a grid of source groups (numbered pink rows) and target groups (numbered cyan columns). The pink number on each core denotes the source group whose synapses it stores. Surrounding panels show one simulation cycle. Update (red): cores advance their neurons by one time step and record which spike. Collocate (yellow): the source group local IDs of spiking neurons are compacted into one buffer per source group, with the spike count in the first slot. Communicate (green): each source group’s buffer is sent to its matching core in every target group. Deliver (blue): each core looks up the synapses of every received ID, adds their weights to its copy of the target group’s spike buffer at the slot given by the delay, and the copies are summed within the target group.

1 Introduction Comprehensive simulation studies of cortical networks depend on compute systems capable of running models significantly faster than biological real-time. Any improvement in simulation efficiency directly benefits neuroscientists in modelling tasks, allowing for more extensive parameter scans and larger scale experiments without the need for massive clusters. For example, capturing learning processes such as long-term plasticity requires tracking network activity over biological hours or days [1], demanding simulation speeds that remain difficult to achieve at scale on current hardware. For 2

a large variety of neural simulations, CPUs remain the most used target platform [1–3]. Increasingly efficient algorithm design [4] has led to a consistent increase of simulators’ speedup relative to biological simulated time Tbio (real-time factor, RTF = Twall /Tbio ) over the years. Across CPU, GPU and SpiNNaker implementations of the same benchmark, however, these gains largely follow hardware improvements. The achieved RTF approximately scales with the square of the process node and linearly with the energy per synaptic event, the only exception being FPGA-based systems [5]. One major bottleneck limiting further speedups lies in the very nature of these simulations: low compute requirements with frequent and often irregular memory access, along with unstructured communication needs between neuron-hosting processes. Due to the separation of memory and compute in the von-Neumann architecture, runtimes are dominated by cores waiting for data movements to complete. This is especially true for simulations of point-neuron models, since the computational load for simple update equations (e.g. Leaky Integrate-and-Fire, LIF) is minimal. GPU-based simulators [6, 7] offer a more compact compute-memory layout than CPUs, enabling slightly more graceful scaling, but ultimately suffer from the same problem. Moreover, especially larger neural models tend to exhibit extremely sparse connectivity. A single neuron is unlikely to have a large number of post-synaptic connections on any given process [4], leading to spike delivery suffering the most from disjoint and uncoalesced memory access. The key to breaking this scaling barrier, therefore, lies in efficient routing. A recent breakthrough was made by the FPGA-based neuroAIx simulator [8], which uses a dedicated router unit and algorithm to simulate the cortical microcircuit benchmark at 20× biological real-time. This comes with flexibility trade-offs typical of FPGA-based implementations, where adopting a different neuron model or network often requires re-synthesising the design. We aim to further improve on this level of performance without trading off flexibility. State-of-the-art technology for distributed simulations of spiking networks on conventional hardware typically advances in discrete time steps and periodically exchanges spikes between compute nodes to preserve causality [9]. The interval over which each node may progress independently is set by the minimum synaptic transmission delay; one such interval constitutes a simulation cycle. In line with NEST terminology (e.g., [10]), we distinguish four phases. In the update phase, every node advances the dynamical state of its local neurons over the minimum-delay interval, optionally in smaller steps, and records which of them cross the firing threshold. In the collocate phase, the resulting spikes are compacted into communication buffers. In the communicate phase, these buffers are exchanged between nodes through a collective operation. In the deliver phase, the synaptic weights are added to the corresponding target neuron’s spike buffer, contributing to the state update at the point in simulated time set by the remaining synaptic transmission delay. We evaluate our approach using the fullscale Potjans-Diesmann cortical microcircuit model (PD14), a standard benchmark comprising approximately 77,000 neurons and 3 × 108 synapses (Sec. 3.1) [5, 11]. For our solution, we make use of distributed memory architectures, systems that allocate fast volatile memory (i.e., SRAM) to every computing core instead of organising it in cache hierarchies or sharing between parallel units. This reduces memory

3

access latency of each processor to its local memory to a small constant duration, regardless of the access pattern, offering orders-of-magnitude lower latency and higher bandwidth compared to external memory (DRAM). Neuromorphic systems often provide dedicated hardware support for spike communication through event-driven or mixed-signal designs, but constrain the model to fit it. SpiNNaker was designed for 1 ms time steps and fixed-point arithmetic, and required a 20× slow-down to simulate PD14 at 0.1 ms resolution [12], with real time reached only through model-specific redistribution of synaptic processing [13]. Loihi [14] and TrueNorth [15] restrict weight resolution, delays and neuron dynamics, while the analog BrainScaleS-2 system [16] operates with limited neuron fan-in, provides 256 synapses per neuron circuit and 131,072 synapses per chip. PD14, with 3880 synapses per neuron and 3 × 108 in total, would therefore require over 2000 chips. In parallel, machine-learning accelerators such as Cerebras [17], NorthPole [18] and the IPU [19] demonstrate the value of colocating memory and compute at scale. We choose the IPU as reference hardware for our implementation, as it offers a large amount of local memory per core and provides a flexible low-level software stack with its custom poplar SDK. In order to make use of low access latency for local memory, we implement a parallel spike delivery algorithm that maximises local computation by structuring spike exchange between groups of cores. Our algorithm defines communication pathways for a configurable group size at build time, thus reducing the required filtering as every post-synaptic target builds a connection lookup table for its pre-synaptic source groups ahead of simulation. Each core is part of both a source and a target group. The source group defines which buffer its spikes will be aggregated into, while a core’s position in a target group defines which source-target group pair it will hold the synapse data for. By tweaking the ratio of groups, we can change the communication granularity and which stages of the simulation benefit most from parallelisation. The optimal configuration will depend on the hardware and implementation of inter-core communication. Combined with the fast random access times of distributed memory, this architecture underpins our state-of-the-art performance. At the same time, the IPU retains the flexibility of a higher-level software-programmable accelerator. Changes to the model or implementation only require recompiling C++ codelets, which takes seconds as opposed to hardware re-synthesis. We achieve simulation times 22.7× faster than biological real-time, meaning a 1.11× speedup and approximately 3.6× better energy efficiency compared to the previous state-of-the-art neuroAIx [8]. Part of this gap reflects that our implementation requires no external memory, but the values also differ in measurement level (IPU device power vs. FPGA board power, Sec. 2.4). Furthermore, we demonstrate that, next to our 32-bit default, the architecture can utilise 16-bit precision synaptic weights to minimise memory consumption while still delivering accurate simulation output. Verification is done through comparison of the same model run on CPU-based simulator NEST [20], see Appendix A.

4

2 Results 2.1 Group-based Algorithm The time-driven simulation advances by repeating the same four phases every cycle (Fig. 1): evolving the neuron states (update ), collecting emitted spikes into communication buffers (collocate ), exchanging these buffers between cores (communicate ), and delivering the spikes to their target neurons (deliver ). The work within the update, collocate, and deliver phases is parallelised, but cores must synchronise across node boundaries in the communicate phase. The key idea of our approach is a group-based mapping of synapses to cores. Functional groups partition the cores such that each core is assigned to one source group and one target group. The C cores form a Gsrc × Gtgt grid, where each source group spans Gtgt cores (one per target group) and vice versa each target group spans Gsrc cores (one per source group). All neurons of the simulated network are distributed round-robin across the source groups (interleaved per neuron population) and, within each source group, across its cores (see Appendix B). The groups determine the mapping of synapses to cores: a target group hosts all incoming synapses of its neurons, while each core of a target group hosts only the synapses originating from a specific source group. This mapping enables a group-based communication scheme and shapes all four phases of the simulation cycle. In the update phase, neuron dynamics evolve, integrating spike contributions and external input as well as detecting threshold crossings of the neuronal membrane potential. Each core is responsible for the update of its local neurons NC = N/C , where N is the total number of neurons in the network. The benchmark network model PD14 (Sec. 3.1) consists of simple current-based leaky integrate-and-fire neurons, allowing for exact integration of the model dynamics [21]. A per-neuron buffer records the emitted spikes of the current time step as the local neuron IDs of the respective source group. In the collocate phase, spikes are hierarchically gathered into communication buffers. A source group manages the joint communication of the spikes emitted by its corresponding Nsrc = NC Gtgt neurons. To this end, every neuron j is assigned a local index nsrc j ∈ {0, 1, 2, . . . , Nsrc − 1}. Each source group has one or multiple root cores, which gather all spikes emitted by the source group in a given time step. This is done over several intermediate aggregation steps, resulting in a final payload containing the number of spikes followed by a contiguous list of local indices. In the communicate phase, the packed spike buffers are exchanged between source and target groups according to the group-based mapping. Each target group receives the buffers of all source groups, while each core within the target group receives only the buffer of its corresponding source group (Fig. 1). In the deliver phase, each core processes the spikes in its received buffer and routes them to their target neurons. Target group cores locally hold all synapse data for the matching source group and each core has a copy of the entire target group’s ring buffer (C, dmax , Ntgt ), a data structure meant to hold spike weights until they can be processed by the target neuron with a simulated delay up to a maximum number of steps, dmax . The copy is necessary to enable parallel routing of spikes. Synapse data comprises weights, delays, and the target group local index ntgt ∈ {0, 1, 2, . . . , Ntgt − 1}, which j

5

identifies the target neuron within the target group. Synapse data (weights, delays, targets) are stored per core in three contiguous arrays sorted by source. This information is accessed using an index range (begin, end) into the synapse data, retrieved from indexing into a lookup table with the incoming spikes’ source group local index nsrc j . We then iterate over the index range and add each synapse weight in its corresponding ring buffer slot, identified by the local index ntgt and an offset from the j current time step index given by the synapse delay. Finally, each target group performs a reduce over its Gsrc ring buffers. This ensures that every target group neuron has a complete view of its received spikes. ntgt j

2.2 Memory Usage and Workload Because the group mapping is independent of the instantiated connectivity, its expected per-core memory footprint can be predicted before graph construction as shown in Sec. 3.6.2 and Fig. 2d). Increasing the number of source groups reduces the lookup-table footprint as each core represents fewer potential source neurons, but enlarges the target groups and therefore the neuronal spike buffers. The synapse-data contribution remains approximately constant. Measured memory usage closely follows this prediction and is minimized at an intermediate group ratio of Gsrc = 16. Small deviations across connectivity seeds arise because, although the connection count between each population pair is fixed, the individual pre- and post-synaptic neuron pairs are sampled randomly, producing seed-dependent differences in the synapse load assigned to each core (see Sec. 3.1). On M = 4 IPUs, configurations with Gsrc = 4 and Gsrc ≥ 64 exceed the memory budget, as minor overhead from runtime structures (see Sec. 3.6.2) pushes slightly over the predicted value in Fig. 2d. Doubling to M = 8 IPUs halves the per-core load and admits Gsrc = 64. Even though grouping determines how synapses are distributed, it does not change the mean synaptic work per thread. Each thread is predicted to process on average 2.7 active synapses per time step for all group configurations (Sec. 3.6.3, Fig. 2e). Scaling groups nonetheless impacts performance. With more source groups, a thread receives fewer spikes, but each spike activates more synapses on it, so that the predicted maxto-mean ratio of active synapses rises from roughly 5 at Gsrc = 4 to 14 at Gsrc = 92, compared with 4.4 for an ungrouped distribution (Fig. 2f). Conversely, the number of spikes each thread has to look up decreases with Gsrc . At worst, a thread scans about 4 spikes per step at Gsrc = 16, compared with 10.4 in the ungrouped case even with ideal target filtering (see Tab. C1). Together with the reduction cost, which grows with the target group size and is similar to the memory scaling, this trade-off between lookup work and synaptic load imbalance favours intermediate group configurations for performance (Fig. 2b).

2.3 Performance We achieve state-of-the-art acceleration of 22.7× over biological real-time, 1.11× faster than the previous record [8] (Fig. 2a). This is done on M = 4 IPUs, the smallest number of nodes that can fit the full-scale microcircuit in memory, and with Gsrc = 16 due to the load balancing and memory trade-offs discussed in Sec. 2.2. Doubling to

6

target groups Gtgt (4 IPUs / 8 IPUs) 736/1472 368/736 184/368 92/184

acceleration factor (1/RTF)

100

4 IPUs 8 IPUs

22.67× 10

0.50

92

64

neuronal spike buffers

predicted 4 IPU nodes

0.25

e

60

synapse data + lookup table

40

deliver

20

communicate

10 expected workload

total

collocate update

spikes active synapses

1

0.1

f

1 CPU 0.1

GPU FPGA (neuroAix)

0.01

FPGA (Lindqvist et al.)

max / mean workload

wall-clock time (ms)

target groups Gtgt 368 184

0.00

0 energy per synaptic event ( J)

736

0.75

1

b

c

1472

d memory per core (MiB)

a

1 8

16 32 source groups Gsrc

64

spikes active synapses

10

4

8

16 32 source groups Gsrc

64

92

Fig. 2 Performance, memory and spike statistics across group configurations. All plots show data for different numbers of source groups Gsrc on the bottom axis, whereas top is the corresponding number of target groups Gtgt . Whenever applicable, error bars are shown for three differing seeds run for Tbio = 15 min each. a. shows speedup factors (1/RTF); the horizontal lines show the best-performing setup for different types of hardware, on both M = 4 and 8 IPU configurations. Top to bottom: FPGA cluster [8], GPU [22] and CPU [23]. b. Breakdown of simulation time into algorithm phases on four IPUs for Tbio = 1 s: Update, Collocate, Communicate and Deliver. The marker indicates total simulation time. c. shows the energy per synaptic event (see Sec. 3.6.5) for 4 and 8 IPU configurations. For FPGAs, the level of measurement (e.g. board, system, device) differs depending on available data and tooling (see Sec. 2.4). d. Per-core memory usage with stacked contributions of different data structures (see legend). Available memory per core on the reference system is indicated by the horizontal dotted grey line. For prediction (open circle) see Sec. 3.6.2. e. Expected spikes per source group (open) and active synapses per thread (filled) per time step. f. Max-to-mean ratio of active synapses per thread (filled, Monte Carlo) and of spikes per source group (open, measured). Dotted line: ungrouped placement (Gsrc = 1). Values in Table C1.

M = 8 IPUs lowers the speedup to at most 15.6×. The benchmarked microcircuit model is computationally lightweight enough that it above all benefits from locality, with neuron updates being the smallest contributor to simulation time. The error bars in Fig. 2a show minimal standard deviation due to the static graph compilation on the IPU and semi-fixed connectivity scheme of the microcircuit (Sec. 3.1).

7

Table 1 Performance (RTF) and energy-efficiency (Esyn ) benchmark results for full-scale PD14 simulations, ordered by publication year and showing only the best reported performance or energy-efficiency results for each architecture. Energy values depend on the measurements reported in the original studies and therefore may be device-, board-, or system-level. For the nearest competitors, footnotes specify which hardware components are included in the power measurements. RTF

Simulator

Ref.

Esyn (µJ)

1.00 SpiNNaker 1 [13] Rhodes et al. (2019) 0.601 0.53 NEST CPU [23] Kurth et al. (2022) 0.48 0.27 GeNN [22] Golosio et al. (2023) 0.074 0.05 neuroAIx [8] Kauth et al. (2023) 0.048† 0.79 horiz/mono/s [24] Lindqvist et al. (2024) 0.021 / 0.009‡ 0.044 IPU This work 0.013 † ‡

# Nodes

System

12 2 1 35 1 4

48x18xARM-968 2 AMD EPYC Rome 7702 RTX 4090 NetFPGA SUME 1 Agilex 7 FPGA Graphcore Bow IPU

Board-level power, including external memory and peripherals. Device power only, excludes host and system overhead.

2.4 Energy For the fastest configuration, three runs with different connectivity seeds consume (11,361.77 ± 5.34) Joule of device energy per run and produce an average of 222,580,247 spikes over Tbio = 15 min. Using the mean neuronal fanout KN = 3880, we obtain 13.16 nJ per synaptic event, where one event denotes the delivery of one emitted spike to one post-synaptic target. At the device level, this is competitive with the record 9 nJ FPGA-fabric-only estimate reported by Lindqvist and Podobas [24]. While energy usage is 3.6× lower than neuroAIx’s reported 47.55 nJ, neuroAIx’s measurement includes external memory and board peripherals, while ours covers aggregate per-IPU device power. This ratio therefore compares reported values with different measurement boundaries. Published results across architectures are compared in Tab. 1 and Fig. 2c. An extended version of the table, including process node and external-drive configuration, is provided in Tab. C2.

2.5 Build and Compilation For the full-scale microcircuit on M = 4 IPUs, most of the build time is spent generating the lookup table used for spike routing, which takes (129.92 ± 2.54) s. The lookup table records how many local targets each incoming source index has for a given source-target group pair. Its construction requires local sorting of target indices to avoid threads accessing the same ring buffer slot during routing, along with alignment steps to account for padding. The Graphcore software stack also requires an internal graph compilation step, which takes (4.21 ± 1.03) s for the same samples. Once built, the computational graph can be stored and loaded from file in subsequent runs, allowing parameters that are not tied to tensor sizing or logic flow, such as simulation duration, to change without repeating these build steps.

8

3 Methods 3.1 Benchmark Network Model Our primary benchmark target is the well-established microcircuit model by Potjans and Diesmann [11]. The model represents 1 mm2 of cortical tissue incorporating experimental connectivity data from rat and cat neocortex, spanning around 300 million synapses. It is subdivided into four layers, each with both excitatory and inhibitory neuron populations. Lateral connectivity is random, meaning there is no inherent locality that can be exploited during spike propagation. Synapse data, meaning primarily weights and delays, are drawn from population-specific distributions. Neuron pairs (pre, post) are randomly picked within population pairs up to a fixed total number derived from the inter-population probabilities. The model is considered a building block of cortical function and features as a basic unit in larger models [25–27]. Due to these characteristics, the microcircuit model has become a stand-in for measuring neural simulator efficiency across platforms and software packages. Senk et al. [5] define a recipe for benchmarking and validation, which we make use of. The performed tests, comparing population mean rates, correlation, and inter-spike intervals, create fundamental computational requirements in maintaining the desired accuracy. Common standards to ensure this are:

• Individual representation of synapse data, which imposes a minimum memory requirement. This includes the synapse target (connectivity), weights, and delays. Although previous work has established that single-float (32-bit) precision weights still produce statistically correct dynamics [28], this per-synapse data still dominates the overall memory requirement for models. Especially on accelerators and other memory-constrained devices, this limits the options for parallel algorithms that would require duplicated data structures per thread. Procedural connectivity generation can reduce this footprint by drawing connections on-the-fly from seeded distributions [29], allowing reproducible simulations of larger static networks at competitive performance. Here, we use explicit synapse data to follow the full-density benchmark specification and to keep the routing and memory requirements directly comparable across simulator backends. • Strictly temporally sequential simulation with a ∆t of 0.1 ms. This significantly restricts the room for parallelisation, as the computational burden of updating exactly integrated [21] current-based (CUBA) LIF neurons is already rather small. The bottleneck of delivery and communication oftentimes makes scaling up the number of processes naively more costly, due to more disjoint memory and physical separation, outweighing gains from update parallelism. Given these requirements, frequent memory access to look up and route incoming spikes and the corresponding synapse data (deliver phase) is unavoidable. On conventional hardware, this requires carefully chosen data layouts to ensure predictable access to allow pre-fetching and minimise cache misses [30]. Even with these optimisations applied, delivery time takes up a significant portion of simulation time on a microcircuit scale. We simulate the full-scale PD14 model at a ∆t = 0.1 ms simulation step: N = 77,169 neurons, approximately 3 × 108 synapses (KN = 3880 per

9

neuron on average), mean firing rate νN = 3.2 Hz, and synaptic delays starting at dmin = 0.1 ms. Unless stated otherwise, runs use M = 4 Bow IPUs (one Bow-2000), providing C = M · CM = 5888 (CM = 1472 per IPU) cores with TC = 6 hardware threads and 624 KiB local SRAM per core [31].

3.2 Reference Distributed Memory System Under the hood, Graphcore’s software stack implements a bulk synchronisation pattern [32]. This scheme divides computation into a local update step, a communication phase where cores can send and request non-local data, and finally a global synchronisation barrier. Random access to core-local SRAM comes at a constant cost of 6 cycles. To hide this instruction latency, the IPU employs a barrel processor with 6 threads, where each is scheduled in round-robin fashion each cycle. Every core is capable of MIMD, meaning every thread can execute a different context. Any IPU program is compiled into a static graph, where codelets (akin to CUDA kernels) represent vertices and edges are data dependencies. The implications of these paradigms are that programs must be reasonably well load balanced and that data, represented by tensors at runtime, is of static shape. The latter requirement creates an algorithmic challenge for efficient implementation of structures that are usually represented dynamically in neural simulations, such as spike buffers. IPUs support single (32-bit) and half (16-bit) precision floating point arithmetic. The software stack extensively makes use of low and mixed precision operations and offers additional utility to alleviate precision loss, for instance via stochastic rounding [33]. The IPU however, both in terms of hardware and software stack, primarily targets machine learning use-cases. Previous work has used it to train and run spiking neural networks in a deep-learning setting [34, 35], whereas our simulation of biologically constrained networks with explicit synapses and delays relies on fine-grained control over data layout and movement. Fortunately, Graphcore provides both a high-level Python front-end and a lower level C++ ecosystem, poplar, which largely fulfils these requirements. To achieve peak performance, we nonetheless rely on optimised handcrafted codelets and algorithms to cover our specific simulation steps and work around certain IPU design choices, further described in Sec. 3.3.

3.3 Algorithm Implementation We implement the four simulation phases introduced in Sec. 2.1 (and visualised in Fig. 1) as a static poplar computational graph using custom C++ codelets. The following subsections describe the data layouts and operations used in each phase. During graph construction, a layout checker verifies that tensors intended to remain core- or group-local are not accessed outside their assigned scope. The neuron state update step of the algorithm is straightforward and fully described by the outline given in Sec. 2.1.

3.3.1 Collocate Collocation is split into local on-core and source group-wide neuron aggregation steps. Aggregation steps can be viewed as a parallel prefix-sum algorithm [36] with a fixed 10

number of stages. Indices of spiking neurons are recorded in integer precision to comply with the minimum word size and to leverage the control-flow-optimised execution pipeline. Core-local Aggregation. First, a single thread per core gathers the source group local index nsrc j of all NC locally hosted neurons that crossed the firing threshold into a local buffer sized bcore . The expected number of spikes per time step and per core for a given model with an average neuronal firing rate of νN is specified by νN NC ∆t. For the microcircuit, for most realistic group combinations, this value falls below one. The buffer must be sized to accommodate the maximum per-time step spike volume. Load balancing measures described in Appendix B help minimise this peak capacity requirement. Source Group Aggregation. In the next stage, all Gtgt bcore potential spike source indices are aggregated into a combined source group-wide buffer of size bsrc . For this, we use all available threads on cagg cores in the centre of the group (assuming linear layout), each writing into fixed offsets bsrc /TC within the output. Every thread denotes the number of elements found si in the 0th slot of its allocation. Buffer Packing. Finally, we run another multithreaded codelet on the same aggregation cores to write the indices into a contiguous buffer for all grouped aggregation cores cagg . This is referred to as the send buffer and forms the actual payload we will communicate to targets in the next step. To do so, every thread calculates the cumuPt−1 lative sum of spikes up to its own ID, defined as i=0 si , and then uses the resulting offset into the final output buffer to copy over indices. A count of all spikes in the group aggregation, ncount , is written to index 0. Viewed across all groups, this results in a buffer of shape (Gsrc , cagg , nfinal ) where nfinal is the upper limit for the count of src spikes per group and each contains Ninc ⊂ Nsrc source group indices.

3.3.2 Communicate The communicate phase transfers the packed send buffers produced during the collocate phase from source groups to target-group cores. Within each target group, a core with kjtgt , j ∈ {0, . . . , Gsrc − 1}, receives the send buffers of the corresponding source group Gsrc,j . The target group collectively receives the buffers of all source groups, while each of its cores receives data from only one source group. On the IPU, this exchange is specified implicitly. Target-side access to a source-group send buffer establishes the required data dependency, and poplar performs the data movement before the deliver computation.

3.3.3 Deliver The deliver phase first maps the source group local neuron indices contained in the received buffers to locally stored synapses and accumulates the weighted events in the target neurons’ ring buffers. Each target-group core stores the synapse data originating from its corresponding source group and maintains a private ring buffer contribution covering all neurons in the target group. In a second step, the duplicate ring buffer views are reduced to ensure each neuron has a full view of all exchanged spikes in a time step.

11

Spike Delivery. Translation of the received buffers’ source group local indices src Ninc to synapse data as well as parallelisation of the whole process is enabled by the lookup table. During build time, we create it as a tensor of shape (Gtgt , Gsrc , Nsrc , TC + 1). The first two dimensions identify a unique target group core by the source-target group pair, the third spans the maximum number of incoming source group local indices and the last splits the table by threads. It thus defines a cumulative range of indices into local synapse data. This thread split must account for neuron boundaries, as splitting the same synapse targets across threads could lead to race conditions. Depending on our choice of (Gsrc , Gtgt ), the lookup table’s maximum range value for local synapse data will be small enough to encode it as 16-bit integers. The TC dimension partitions the local target neurons Ntgt into contiguous blocks, one per thread, rather than partitioning synapses directly. Due to the incoming send buffers’ shape, (Gsrc , cagg , nfinal ), an outer and an inner loop are required for delivery of all spikes. The effective size of the group aggregation dimension ncount is stored at index 0, and serves as the limit for the inner loop to avoid iterating over empty buffer slots. Finally, we can use the index range retrieved from the lookup table to store synapse weights at position (tnow + d) mod dmax in the ring buffer, where tnow is a running index in range [0, dmax ) signalling the current time step and d is the synapse delay used to offset. In the microcircuit model, the synaptic time constant τsyn is the same for excitatory and inhibitory connections, which allows us to merge both types of possible weights into one dimension in the ring buffer. Reduce. Each target group core handles a unique (Gsrc,i , Gtgt,j ) combination in terms of spike routing. However, since one source group’s send buffer is sent to every target group, the buffers inevitably contain source group indices matching neuron targets beyond what is local on any one core. Conversely, one target group core only handles spikes from one source group, but must be able to receive spikes from all. For this reason, every core must be able to accept spikes with targets located on any other core within its target group. This is reflected in the shape of the ring buffer (C, dmax , Ntgt ), which must be reduced along the core dimension for d = tnow each time step. The Gsrc per-core copies of shape (dmax , Ntgt ) are summed element-wise at slot tnow , yielding one input value per target neuron. The operation is implemented using a poplar library method which leverages manual assembly codelets. Since group sizes are aligned with the IPU node boundary, this reduce operation is always IPU-local and therefore incurs no global communication.

3.4 Modelling Choices In the standardised version of the microcircuit model1 , synaptic delays are drawn from a truncated normal distribution starting at the minimum delay and unbounded above [dmin , ∞). When building the computational graph and compiling the codelets on the IPU, only fixed-size data structures are allowed. As the maximum delay dmax impacts the sizing of the ring buffer (C, dmax , Ntgt ), the maximum delay has a strong effect on the total memory usage. No delays are truncated in our implementation, but dmax directly scales the ring buffer memory (Sec. 3.6.2). Bounding the delay distribution offers a way to fit models slightly exceeding the available memory, at the cost of a small 1

https://github.com/INM-6/microcircuit-PD14-model/

12

deviation from the model’s specification. Synaptic weights could also be stored in half precision, which halves their memory footprint and passes the statistical validation (Appendix A). However, since the routing codelet processes spikes scalarly and the IPU’s minimum word size is 32 bits, the required realignment cancels out the speedup of the reduction, so all reported results use full precision.

3.5 Codelet Optimisation Graphcore’s codelet compiler popc is based on LLVM [37], which allows for interpretation and iterative optimisation of performance-critical code by compiling to an intermediate representation and checking the issued instructions. IPU memory is 64 bits wide, with a 2-way interleaved region that additionally supports 128-bit loads and stores [38], and we vectorise wherever the data layout permits. All input and output data is defined as tensors via poplar and mapped to a one-dimensional naturally aligned vector in memory for processing in codelets. The neuron update processes neuron pairs with 64-bit float2 arithmetic, and the reduce step reads and clears the current ring buffer slice with 128-bit float4 accesses before summing it across source groups with poplar’s vectorised reduction. Aggregation and routing codelets benefit less. Spike indices are integers, and the IPU’s integer pipeline lacks vector instructions [38]: the aggregation step can scan its input in int4 blocks to skip empty entries, but compacting the remaining spikes is scalar. Routing performs data-dependent scatter-adds into the ring buffer, indexed by target neuron and delay, and therefore remains scalar regardless of data type. Spike delivery is consequently the most expensive phase (Fig. 2b). Since the scatter-add cannot be vectorised, its cost per step is set by the number of active synapses on the most loaded thread, which motivates the workload analysis in Sec. 3.6.3.

3.6 Benchmarking and Analysis To assess whether the benchmark model can be simulated on the reference system and to characterise the expected workload and imbalance, we provide analytical predictions of the per-core memory footprint and the per-thread computational load. The analysis is supported by Monte-Carlo estimates of the load imbalance as well as measured data. The predictions are evaluated as a function of the source group count (see Fig. 2), which controls how neurons and their synapses are partitioned across cores (Sec. 3.3).

3.6.1 Setup The reference system (Sec. 3.2) comprises M = 4 IPU nodes, each with CM = 1472 cores and TC = 6 hardware threads per core, giving a total of C = M CM = 5888 cores and T = C TC = 35,328 threads. The benchmark network is a cortical microcircuit model (Sec. 3.1), with N = 77,169 neurons and on average KN = 3880 incoming (equivalently outgoing) synapses per neuron. We assume a mean firing rate of νN = 3.2 Hz and a simulation resolution of ∆t = 0.1 ms. The minimum synaptic delay equals the simulation step size, which means that spikes need to be communicated every step. The predictions do not take into account the different neuronal populations and their specific firing rates and connectivity. 13

3.6.2 Memory Usage Prediction. On average NC = ⌈N/C⌉ neurons and KC = NC KN synapses are placed on each core. The cores are partitioned into Gsrc source groups and Gtgt = C/Gsrc target groups. We consider Gsrc ∈ {4, 8, 16, 32, 64, 92}. A source group comprises Nsrc = Gtgt NC neurons and a target group holds Ntgt = Gsrc NC neurons. For the prediction of the per-core memory usage, we take into account the contributions of synapse data (weights, delays, targets), the lookup table storing synapse-data positions, and the neuronal spike buffer (Fig. 2d). We assume that each synapse takes up msyn = 7 Byte of memory, which is composed of a mw = 4 B weight, a 1 B delay, and a 2 B target group local neuron ID. The lookup table, viewed per-core, has two dimensions and holds TC + 1 local connection indices of mlcid = 2 B for each neuron of the sending source group. Altogether, synapse data and lookup table require Msyn = Nsrc (TC + 1) mlcid + KC msyn .

(1)

The two-dimensional neuronal spike buffer reserves one weight-sized slot per target group local neuron and per potential delay step, where we assume a maximum of dmax = 60 steps. This yields a memory usage of

Mbuf = Ntgt dmax mw

(2)

for the neuronal spike buffers. In total, the predicted per-core memory amounts to

M = Msyn + Mbuf ,

(3)

which we compare against the available budget of approximately 624 KiB (Fig. 2d). Measurement. The actual memory measurements include a few quantities omitted from the prediction because they are difficult to predict ahead of time. This includes, for instance, code storage, live exchange buffers (created automatically when data is communicated) and padding for alignment purposes. These contributions are relatively minor, which is why the error bars in Fig. 2d are shown to line up closely with the prediction. Memory is recorded using the PopVision Analysis Library (libpva), with which we programmatically read the mean memory of Gsrc cores (one aggregation core per group). We then take the mean and standard deviation across three connectivity seeds.

3.6.3 Workload Prediction. Each source group has target synapses on one core of each target group. Hence, the synapses of one source group are distributed over Ttgt = Gtgt TC

(4)

destination threads. The probability for a neuron to spike in a time step is defined as

p = νN ∆t. 14

(5)

The expected spike count per source group, as shown in Fig. 2e, then is

Ssrc =

Np . Gsrc

(6)

Increasing the number of source groups reduces the expected number of spikes per source group, Ssrc , but concentrates each neuron’s KN synapses on fewer cores, since a source group spans only Gtgt = C/Gsrc cores. Both effects cancel out exactly, so that the expected number of active synapses per thread and step can be written as:

A = Ssrc

KN N p KN Gsrc N KN p = · = ≈ 2.7 Ttgt Gsrc T T

(7)

and is independent of the group configuration (see Fig. 2e). Each thread also scans all Ssrc spikes of its source group in the deliver loop, so the lookup work per thread decreases with Gsrc . Since the IPU’s cores must sync after each operation (see Sec. 3.2), a simulation step is only as fast as its slowest parallel operation. So the relevant performance bottleneck is the expected maximum load relative to the mean: E[maxthread Athread ]/A. Although the mean number of active synapses per thread is constant, a larger Gsrc means that each target thread receives fewer spikes (fewer cores per corresponding source group), but each spike activates more synapses on it (more cores in the target group it belongs to). The per-thread load therefore fluctuates more strongly, and its maximum over all threads grows, as shown in Fig. 2f. The figure’s ratio of active synapses E[maxthread Athread ]/A is estimated by Monte Carlo simulation over 100 iterations of 1000 time steps each, drawing each neuron’s KN target threads uniformly from the Ttgt threads of its source group and its spikes independently with probability p. As a reference, the same procedure is applied to an ungrouped distribution in which synapses are drawn from all T threads (Gsrc = 1). Measurement. The analytical solutions as well as Monte Carlo simulations make some idealised assumptions, such as a homogeneous, uncorrelated firing rate and a uniform split of synapses over target threads. In our simulation code, in order to avoid overlapping write access, each target core thread k ∈ {0, . . . , TC − 1} is assigned a fixed slice of possible local neuron IDs:



 Gsrc NC nk,start = k . TC

(8)

Since neurons of each population are interleaved across source groups, each block contains a mix of populations, so the impact of this difference is expected to be small. Empirical spike counts per source group (and the resulting max-mean ratios) are obtained from recorded spike trains (spike output enabled, untimed runs of 15 min biological time, 3 seeds), binned per time step and source group. Per-thread active synapse counts were not manually recorded.

15

3.6.4 Simulation Time All results are measured by recording cycle counts using IPU native methods over the whole simulation duration. Cycle counts were converted to wall-clock time using the IPU device clock frequency f = 1.85 GHz. Real-time factors are reported as RTF = Twall /Tbio , where values below one indicate faster-than-real-time simulation. The primary codelets involved in a simulation cycle are additionally timed to allow for a time breakdown. Performance measurements generally exclude copying spikes to host. Spike-train output is enabled only for (untimed) validation runs, as is standard practice for acceleration comparison [5]. Static buffer sizes (e.g. bcore , nfinal ) and the aggregation granularity cagg were tuned per configuration by grid search. Each configuration was first run in a debug build that checks all aggregation buffers for overflow and aborts if one occurs. Timed runs use the same seeds and durations with these checks disabled.

3.6.5 Energy Energy usage was measured using the IPU’s built-in gcipuinfo command, which reads out the momentary per-IPU power consumption. In order to periodically record values during a simulation we use the JSC power tool2 described in [39], which spawns a monitoring thread on the host and aggregates the power measurements. Running this over 3 seeds for full-scale simulations of Tbio = 15 min biological time for the fastest configuration and querying power every ∆tsample = 0.1 s of runtime, we can calculate total energy Er consumed by integrating per-IPU power P over IPUs k and individual samples i for one run r: nr X M X Er = ∆tsample Pr,k,i . (9) i=1 k=1

Based on standard practice [12, 40], the number of synaptic events was estimated from the number of emitted spikes Sr and the mean neuronal fanout KN , with which we calculate the energy per synaptic event:

Esyn,r =

Er . Sr KN

(10)

Reported values are the mean and standard deviation of Esyn,r across the three runs.

3.6.6 Build Time Simulation parameters and model data are specified through a Python frontend. Synapse and neuron inputs are provided in a structure of arrays (SoA) format, reorganised for the chosen C++ layout, and converted to runtime tensors. As a result, changes to hyperparameters or input models do not require C++ re-compilation. This keeps the edit-run cycle short compared with static simulators that achieve comparable simulation speedups only after substantially longer compilation steps, such as FPGA-based systems [8, 41]. Because the model implementation is decoupled from the 2

https://github.com/FZJ-JSC/jpwr

16

hardware-specific graph construction, new synapse or neuron models only require software re-compilation, adding seconds of build time at most. Build times were recorded for the IPU’s host system, featuring an AMD EPYC 7742 processor.

4 Discussion In this work, we demonstrate that distributed on-chip memory can move spiking network simulation beyond the performance envelope of conventional many-core systems. Our group-based routing strategy maps spike exchange onto structured source-target communication, reducing runtime filtering while preserving the flexibility of a programmable accelerator. Applied to the full-scale cortical microcircuit, the resulting IPU simulator reaches 22.7× biological real-time, validates against NEST results, and combines high speed with short build times and competitive energy efficiency. These results establish distributed-memory hardware as a compelling platform for largescale computational neuroscience, and provide an algorithmic template for exploiting locality in future neural simulators. While the microcircuit achieves peak performance at just M = 4 IPUs, this can mainly be attributed to the low computational load, which makes a compact fit in memory and short communication pathways crucial. Since on-chip memory grows linearly with the number of IPUs (about 0.9 GB each [31]), larger systems provide tens of gigabytes of SRAM. This capacity may benefit models that make more use of update parallelism or have spatially constrained activity, for which the additional communication cost of spreading across more nodes, as seen here when moving from M = 4 to M = 8, is less dominant. When respecting the modelling requirement that synapse data must be explicitly stored on the device, the next largest benchmark models, such as the mesocircuit [27] or the multi-area model [42], are at least an order of magnitude larger in neuron count and therefore hard to fit into purely local memory. This is a capacity limit of the present explicit-connectivity implementation rather than a limitation of the groupbased routing strategy itself. For long-running studies with evolving connectivity, such as structural plasticity experiments, the IPU is nevertheless a promising target: its achieved speedup directly reduces the wall-clock cost of biologically long simulations, while hardware-supported random number generation [43] can make stochastic target selection and connectivity updates comparatively cheap. The IPU’s static tensor model means that such dynamics would still need bounded data structures, for example through over-allocation or plasticity rules with explicit upper limits, but these constraints are compatible with a broad class of structural plasticity experiments. At 22.7× real time, a full day of biological time of the microcircuit is simulated in about an hour on M = 4 IPUs. This makes the platform a natural candidate for plastic networks, whose biologically relevant timescales of days to weeks remain prohibitive for conventional simulators.

17

Acknowledgements The presented conceptual and algorithmic work is part of our long-term collaborative project to provide the technology for neural systems simulations3 . This research used resources of the Argonne Leadership Computing Facility, a U.S. Department of Energy (DOE) Office of Science user facility at Argonne National Laboratory and is based on research supported by the U.S. DOE Office of Science-Advanced Scientific Computing Research Program, under Contract No. DE-AC02-06CH11357. This work is partly funded by Volkswagen Foundation - 0071143 (STRUCTICITY) and the German Research Foundation (DFG) - 545776403/FOR5880 (Mod4Comp). We would further like to thank Jette Oberländer for her help in reviewing this manuscript.

Competing interests The authors declare no competing interests.

References [1] Zenke, F. & Gerstner, W. Limits to high-speed simulations of spiking neural networks using general-purpose computers. Frontiers in Neuroinformatics Volume 8 - 2014 (2014). URL https://doi.org/10.3389/fninf.2014.00076. [2] Diesmann, M. & Gewaltig, M.-O. NEST: An environment for neural systems simulations. Forschung und wissenschaftliches Rechnen, Beiträge zum HeinzBilling-Preis 58, 43–70 (2001). [3] Hines, M. L. & Carnevale, N. T. The neuron simulation environment. Neural computation 9, 1179–1209 (1997). URL https://doi.org/10.1162/neco.1997.9.6. 1179. [4] Jordan, J. et al. Extremely scalable spiking neuronal network simulation code: from laptops to exascale computers. Frontiers in neuroinformatics 12, 2 (2018). URL https://doi.org/10.3389/fninf.2018.00002. [5] Senk, J. et al. Constructive community race: full-density spiking neural network model drives neuromorphic computing. Neuromorphic Computing and Engineering 6, 012001 (2026). URL https://doi.org/10.1088/2634-4386/ae379a. [6] Knight, J. C., Komissarov, A. & Nowotny, T. PyGeNN: a python library for GPUenhanced neural networks. Frontiers in Neuroinformatics 15, 659005 (2021). URL https://doi.org/10.3389/fninf.2021.659005. [7] Golosio, B. et al. Fast simulations of highly-connected spiking cortical models using GPUs. Frontiers in Computational Neuroscience Volume 15 - 2021 (2021). URL https://doi.org/10.3389/fncom.2021.627620. 3

https://www.nest-initiative.org

18

[8] Kauth, K., Stadtmann, T., Sobhani, V. & Gemmeke, T. neuroAIx-Framework: design of future neuroscience simulation systems exhibiting execution of the cortical microcircuit model 20× faster than biological real-time. Frontiers in Computational Neuroscience Volume 17 - 2023 (2023). URL https://doi.org/ 10.3389/fncom.2023.1144143. [9] Morrison, A., Mehring, C., Geisel, T., Aertsen, A. & Diesmann, M. Advancing the boundaries of high connectivity network simulation with distributed computing. Neural Computation 17, 1776–1801 (2005). URL https://doi.org/10.1162/ 0899766054026648. [10] Lober, M., Diesmann, M. & Kunkel, S. Exploiting network topology in brainscale simulations of spiking neural networks. Neuromorphic Computing and Engineering 6, 024024 (2026). URL https://doi.org/10.1088/2634-4386/ae762e. [11] Potjans, T. C. & Diesmann, M. The Cell-Type Specific Cortical Microcircuit: Relating Structure and Activity in a Full-Scale Spiking Network Model. Cerebral Cortex 24, 785–806 (2012). URL https://doi.org/10.1093/cercor/bhs358. [12] van Albada, S. J. et al. Performance comparison of the digital neuromorphic hardware SpiNNaker and the neural network simulation software NEST for a fullscale cortical microcircuit model. Frontiers in Neuroscience Volume 12 - 2018 (2018). URL https://doi.org/10.3389/fnins.2018.00291. [13] Rhodes, O. et al. Real-time cortical simulation on neuromorphic hardware. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 378, 20190160 (2019). URL https://doi.org/10.1098/rsta. 2019.0160. [14] Davies, M. et al. Loihi: A neuromorphic manycore processor with on-chip learning. IEEE Micro 38, 82–99 (2018). URL https://doi.org/10.1109/MM.2018. 112130359. [15] Akopyan, F. et al. TrueNorth: Design and tool flow of a 65 mW 1 million neuron programmable neurosynaptic chip. IEEE Transactions on ComputerAided Design of Integrated Circuits and Systems 34, 1537–1557 (2015). URL https://doi.org/10.1109/TCAD.2015.2474396. [16] Pehle, C. et al. The brainscales-2 accelerated neuromorphic system with hybrid plasticity. Frontiers in Neuroscience Volume 16 - 2022 (2022). URL https: //doi.org/10.3389/fnins.2022.795876. [17] Lie, S. Cerebras architecture deep dive: First look inside the hw/sw co-design for deep learning: Cerebras systems. 2022 IEEE Hot Chips 34 Symposium (HCS), 1–34 (IEEE Computer Society, 2022). URL https://doi.org/10.1109/HCS55958. 2022.9895479.

19

[18] Modha, D. S. et al. Neural inference at the frontier of energy, space, and time. Science 382, 329–335 (2023). URL https://doi.org/10.1126/science.adh1174. [19] Jia, Z., Tillman, B., Maggioni, M. & Scarpazza, D. P. Dissecting the graphcore IPU architecture via microbenchmarking. arXiv preprint arXiv:1912.03413 (2019). URL https://doi.org/10.48550/arXiv.1912.03413. [20] Schöfmann, C. M. et al. NEST 3.10. Zenodo (2026). URL https://doi.org/10. 5281/zenodo.18268837. [21] Rotter, S. & Diesmann, M. Exact digital simulation of time-invariant linear systems with applications to neuronal modeling. Biological cybernetics 81, 381– 402 (1999). URL https://doi.org/10.1007/s004220050570. [22] Golosio, B. et al. Runtime construction of large-scale spiking neuronal network models on gpu devices. Applied Sciences 13 (2023). URL https://doi.org/10. 3390/app13179598. [23] Kurth, A. C., Senk, J., Terhorst, D., Finnerty, J. & Diesmann, M. Sub-realtime simulation of a neuronal network of natural density. Neuromorphic Computing and Engineering 2, 021001 (2022). URL https://doi.org/10.1088/2634-4386/ ac55fc. [24] Lindqvist, B. A. & Podobas, A. Algorithms for fast spiking neural network simulation on FPGAs. IEEE Access 12, 150334–150353 (2024). URL https: //doi.org/10.1109/ACCESS.2024.3479933. [25] Schuecker, J., Schmidt, M., van Albada, S. J., Diesmann, M. & Helias, M. Fundamental activity constraints lead to specific interpretations of the connectome. PLOS Computational Biology 13, e1005179 (2017). URL https://doi.org/10. 1371/journal.pcbi.1005179. [26] Schmidt, M. et al. A multi-scale layer-resolved spiking network model of restingstate dynamics in macaque visual cortical areas. PLOS Computational Biology 14, e1006359 (2018). URL https://doi.org/10.1371/journal.pcbi.1006359. [27] Senk, J., Hagen, E., van Albada, S. J. & Diesmann, M. Reconciliation of weak pairwise spike–train correlations and highly coherent local field potentials across space. Cerebral Cortex 34, bhae405 (2024). URL https://doi.org/10.1093/cercor/ bhae405. [28] Knight, J. C. & Nowotny, T. GPUs outperform current HPC and neuromorphic solutions in terms of speed and energy when simulating a highly-connected cortical model. Frontiers in Neuroscience Volume 12 - 2018 (2018). URL https://doi.org/10.3389/fnins.2018.00941.

20

[29] Knight, J. C. & Nowotny, T. Larger GPU-accelerated brain simulations with procedural connectivity. Nature Computational Science 1, 136–142 (2021). URL https://doi.org/10.1038/s43588-020-00022-7. [30] Pronold, J. et al. Routing brain traffic through the von neumann bottleneck: Efficient cache usage in spiking neural network simulation code on general purpose computers. Parallel computing 113, 102952 (2022). URL https://doi.org/10. 1016/j.parco.2022.102952. [31] Graphcore. Bow-2000 IPU-Machine datasheet (2023). URL https://docs. graphcore.ai/projects/bow-2000-datasheet/en/latest/. Accessed: 2026-10-11. [32] Valiant, L. G. A bridging model for parallel computation. Communications of the ACM 33, 103–111 (1990). URL https://doi.org/10.1145/79173.79181. [33] Forsythe, G. E. Reprint of a note on rounding-off errors. SIAM review 1, 66 (1959). URL https://doi.org/10.1137/1001011. [34] Sun, P.-S. V. et al. Exploiting deep learning accelerators for neuromorphic workloads. Neuromorphic Computing and Engineering 4, 014004 (2024). URL https://doi.org/10.1088/2634-4386/ad2373. [35] Finkbeiner, J., Gmeinder, T., Pupilli, M., Titterton, A. & Neftci, E. Harnessing manycore processors with distributed memory for accelerated training of sparse and recurrent models. Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 38, 11996–12005 (2024). URL https://doi.org/10.1609/aaai.v38i11. 29087. [36] Hillis, W. D. & Steele Jr, G. L. Data parallel algorithms. Communications of the ACM 29, 1170–1183 (1986). URL https://doi.org/10.1145/7902.7903. [37] Lattner, C. & Adve, V. LLVM: A compilation framework for lifelong program analysis and transformation. CGO, 75–88 (San Jose, CA, USA, 2004). URL https://doi.org/10.1109/CGO.2004.1281665. [38] Graphcore. Tile vertex instruction set architecture (ISA): Ipu21. Tech. Rep. Version 1.3.1, Graphcore (2022). URL https://docs.graphcore.ai/projects/isa/ en/latest/ static/TileVertexISA-IPU21-1.3.1.pdf. Accessed: 2026-05-05. [39] John, C. M., Nassyr, S., Penke, C. & Herten, A. Performance and power: Systematic evaluation of ai workloads on accelerators with caraml. SC24-W: Workshops of the International Conference for High Performance Computing, Networking, Storage and Analysis, 1164–1176 (IEEE, 2024). URL https://doi.org/10.1109/ SCW63240.2024.00158.

21

[40] Sharp, T., Galluppi, F., Rast, A. & Furber, S. Power-efficient simulation of detailed cortical microcircuits on SpiNNaker. Journal of Neuroscience Methods 210, 110–118 (2012). URL https://doi.org/10.1016/j.jneumeth.2012.03.001. Special Issue on Computational Neuroscience. [41] Heittmann, A. et al. Simulating the cortical microcircuit significantly faster than real time on the IBM INC-3000 neural supercomputer. Frontiers in neuroscience 15, 728460 (2022). URL https://doi.org/10.3389/fnins.2021.728460. [42] Schmidt, M., Bakker, R., Hilgetag, C. C., Diesmann, M. & Van Albada, S. J. Multi-scale account of the network structure of macaque visual cortex. Brain Structure and Function 223, 1409–1435 (2018). URL https://doi.org/10.1007/ s00429-017-1554-4. [43] Hanlon, J. & Felix, S. A fast hardware pseudorandom number generator based on xoroshiro128. IEEE Transactions on Computers 72, 1518–1524 (2022). URL https://doi.org/10.1109/TC.2022.3204226.

Appendix A

Verification of implementation

We verify the correctness of spike trains by applying the established statistical tests described in [5] to the output of a 15 min biological time full-scale simulation after a five second warm-up period, during which spikes are discarded. We generate the microcircuit connectivity for a random seed on the IPU host side and plot the probability distributions of the resulting metrics in comparison to a random seed run performed on the CPU version of NEST 3.10 [20], shown in Fig. A1. Although we do not use 16-bit floating point weights for our simulation results, we include half precision with stochastic rounding to show that it is a viable choice for neural simulations on accelerators with dedicated support for 16-bit arithmetic. We also include the Earth Mover’s Distance (EMD) between the full precision IPU run and the NEST baseline. Per population and their respective neurons, we plot the firing rate over the whole simulation duration, the coefficient of variation between spike timings as well as the Pearson correlation coefficient. For the latter, we subsample 200 randomly chosen neurons in each population and analyse spike trains binned at 2 ms intervals.

Appendix B

Load Balancing

The bulk-synchronous computation scheme of the IPU requires all cores to sync after computation, meaning cores that consistently see a higher workload will cause undue idle time. The aggregation and route step are both dependent on the number of spikes per source group and therefore stand to benefit from neuron placement that takes the group-based routing scheme into account. Outgoing connection density and firing rate are strongly correlated with neuron populations, which is why the naive round-robin placement of neurons often achieves suboptimal load balancing when a small (Npop ≪ C ) but highly active population is primarily mapped into one or very few source 22

A

L23E

EMD: 0.0184

L23I

EMD: 0.0380

B

EMD: 0.0014

L23E

EMD: 0.0012

L4E

EMD: 0.0024

L5E

EMD: 0.0070

L6E

EMD: 0.0024

2

0.5

L23I

C

EMD: 0.0004

L23E

EMD: 0.0001

L4E

EMD: 0.0003

L23I

IPU (fp32) IPU (fp16)

p(a.u.)

NEST CPU

0

5

10

0

10

L4E

EMD: 0.0900

L5E

EMD: 0.3254

20

0

L4I

EMD: 0.0008

1

0.75

1.00

1.25

L4I

0.00

0.02

EMD: 0.0001

0.00

0.01

L4I

p(a.u.)

EMD: 0.0277

0

20

0

20

0.5

1.0

L5I

EMD: 0.0045

L6I

EMD: 0.0028

0.5

1.0

1.5

0.00

L5I

EMD: 0.0002

0.01

L6I

EMD: 0.0001

L5E

0.00

0.01

L5I

EMD: 0.0002

p(a.u.)

EMD: 0.2845

0

20

40

L6E

0

20

EMD: 0.0807

0.5

1.0

0.5

1.0

0.00

0.02

L6E

0.00 EMD: 0.0003

0.01

L6I

p(a.u.)

EMD: 0.0215

0

10

0

Firing Rate [Hz]

20

40

0

1

Coefficient of Variation

1.0

0.000

0.025

0.050

0.00

0.01

Pearson Correlation

Fig. A1 Probability density plot of the three metrics (A: firing rates, B: coefficient of variation and C: Pearson correlation) chosen to validate simulation results plotted for every cortical layer and excitatory/inhibitory population. The baseline run from an established CPU simulator, NEST, is shown in black along with our data for both full- (red) and half (blue, dotted) precision. Earth Mover’s Distance (EMD) between full precision IPU and NEST is labelled in the top right for each population.

groups. Furthermore, because each core in our delivery scheme holds the synapse data for a unique (Gsrc,i , Gtgt,j ) group combination, highly heterogeneous models in terms of outgoing connectivity may run out of memory on select cores using global round-robin. A population- and group-aware approach therefore yields fewer stalls, a more even memory distribution across cores and, in some cases, allows for smaller buffer sizes in aggregation due to a lower maximum spike count per source group in each time step. We apply a simple source group round-robin scheme, in which neurons are placed interleaved among source groups rather than filling up cores one after another (Fig. B2). This massively improves source group spike rate balance and reduces the spread in the number of synapses accessed across cores during routing. This enables smaller static buffer sizes, which is the primary source of speedup from the source group round-robin scheme. While the balanced placement would suggest a more noticeable acceleration factor, actual performance objectives show that spike count extrema are preserved in both versions and merely more evenly distributed in the source group balanced case. This most likely comes from the populations of the cortical microcircuit model already

23

Table C1 Monte Carlo estimates of spike and synaptic workload per simulation step for the grouped placement and an ungrouped reference (Gsrc = 1). Mean and maximum are averaged over 100 iterations of 1000 time steps each. For grouped configurations, spike counts are given per source group, which every thread of the matching target-group cores scans. Active synapses are counted per thread assuming a uniform split. Spikes

∗

Active synapses

Gsrc

Gtgt

mean

max

max/mean

mean

max

max/mean

1∗

–

2.57

10.40

4.05

2.71

11.80

4.35

4 8 16 32 64 92

1472 736 368 184 92 64

6.18 3.09 1.54 0.77 0.39 0.27

8.79 5.77 4.05 3.03 2.39 2.17

1.42 1.87 2.62 3.92 6.20 8.09

2.71 2.72 2.71 2.71 2.71 2.71

13.46 15.29 18.31 23.31 31.62 38.11

4.96 5.63 6.75 8.60 11.66 14.06

Hypothetical ungrouped synapses baseline, drawn uniformly from all T threads.

being close in mean firing rate, along with single-neuron-level bursts and differences in outgoing connectivity that were not balanced for in this scheme.

Appendix C

Supporting Performance Data

Tab. C1 shows the predicted values for active synapses and spike counts across group configurations as used in Fig. 2e–f. Tab. C2 lists additional RTF runs, including process node design as well as the external drive provided. Replacing the Poisson background current with an equivalent DC compensation preserves output statistics in the PD14 model [11].

24

mean neuron spikes / s

Standard round-robin, firing rate CV=0.101

A

2 0

0

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

source group 11.0 0 4 8 12

0-91

0 4 8 12

92-183

0 4 8 12

184-275

0 4 8 12

276-367

10.5

276

298

321

344

367

mean neuron spikes / s

target group (wrapped) Source-group round-robin, firing rate CV=0.015

B

10.0

weighted mean addr span

source group

Weighted mean active synapse count, CV=0.036

2 0 0

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

source group

source group

Weighted mean active synapse count, CV=0.024 0 4 8 12

0-91

0 4 8 12

92-183

0 4 8 12

184-275

0 4 8 12

276-367

9.5

276

298

321

344

target group (wrapped)

367

9.0

Fig. B2 Distribution of mean per-neuron firing rate across source groups and the mean active synapse count (the expected number of targets) for a source neuron on any (Gsrc,i , Gtgt,j ) combination. This is shown in A. for round-robin across all cores and B. round-robin within source group boundaries, respectively. All data measured over one second of biological time on identical simulation seeds.

25

26

Simulator NEST CPU GeNN SpiNNaker 1 SpiNNaker 1 NEST GPU GeNN CsNN NEST CPU NEST CPU NEST GPU GeNN neuroAIx horiz/mono/s IPU

RTF 2.46 1.83 20.00 1.00 1.05 0.70 0.25 0.53 0.67 0.38 0.27 0.05 0.79 0.044

[12] [28] [12] [13] [7] [6] [41] [23] [23] [22] [22] [8] [24] This work

Ref. 9.94 0.47 4.4 0.601 0.25 0.284 0.48 0.33 0.104 0.074 0.048 0.021 / 0.009 0.013

Esyn (µJ) 12 1 6 12 1 1 345 2 1 1 1 35 1 4

# Nodes

2 Intel Xeon E52680v3 Tesla V100 48x18xARM-968 48x18xARM-968 RTX 2080 Ti Titan RTX IBM INC-3000 2 AMD EPYC Rome 7702 2 AMD EPYC Rome 7702 RTX 4090 RTX 4090 NetFPGA SUME 1 Agilex 7 FPGA Graphcore Bow IPU

System

22 12 130 130 12 12 28 7 7 5 5 28 10 7

Process (nm)

DC Poisson DC DC Poisson Poisson Poisson DC DC DC Poisson DC Poisson DC

Drive

Table C2 Extended version of published benchmark results, adapted from Senk et al. [5]. Colours highlight similar hardware: CPU (pink), GPU (blue), FPGA (light blue), SpiNNaker (pastel blue) and IPU (dark blue).

Related documents

• Fast and Flexible Spiking Network Simulations on Massively Distributed On-Chip Memory Accelerators#1343733
• SupraSNN: Exploiting Synapse-Level Parallelism in Spiking Neural Network Accelerators through Co-Optimized Mapping and Scheduling#271789
• NeuroTrain: Surveying Local Learning Rules for Spiking Neural Networks with an Open Benchmarking Framework#187353
• Integrated Neuromorphic Photonic Computing for AI Acceleration: Emerging Devices, Network Architectures, and Future Paradigms.#671216
• The nanoscale mobility of calcium channels is driven by readily releasable synaptic vesicles to support precise neurotransmission in live C. elegans#146488