Paper deep dive
SuperNeuroMAT: An Efficient Matrix-based Simulator for Spiking Neural Networks
Prasanna Date, Kevin Zhu, Shruti Kulkarni, Ashish Gautam, Chathika Gunaratne, Robert Patton, Tyler Nitzsche, Ian Mulet, Zachary Johnson-Scott, Addison Helms, Duncan Rowden, Simon Weston, Maryam Parsa, Catherine Schuman, Thomas Potok
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 92%
Last extracted: 8/12/2026, 1:47:46 AM
Summary
The paper introduces SuperNeuroMAT, an open-source, Python-based simulator for Spiking Neural Networks (SNNs) that utilizes a novel matrix-based approach to model leaky integrate-and-fire (LIF) neuron dynamics. It supports both dense and sparse execution modes, enabling efficient simulation of up to 10,000 neurons in dense mode and 100,000 in sparse mode on standard hardware. SuperNeuroMAT outperforms established simulators like NEST, Brian2, BindsNET, and snnTorch in execution speed and memory usage, and is applicable to machine learning benchmarks, neuromorphic vision tasks, and general-purpose algorithms.
Entities (30)
Relation Signals (20)
Maryam Parsa → affiliatedwith → George Mason University
confidence 95% · Maryam Parsa 2 ... 2 George Mason University
Simon Weston → affiliatedwith → Oak Ridge High School
confidence 95% · Simon Weston 4 ... 4 Oak Ridge High School
Duncan Rowden → affiliatedwith → Oak Ridge High School
confidence 95% · Duncan Rowden 4 ... 4 Oak Ridge High School
Catherine Schuman → affiliatedwith → University of Tennessee, Knoxville
confidence 95% · Catherine Schuman 3 ... 3 University of Tennessee, Knoxville
Addison Helms → affiliatedwith → University of Tennessee, Knoxville
confidence 95% · Addison Helms 3 ... 3 University of Tennessee, Knoxville
Zachary Johnson-Scott → affiliatedwith → University of Tennessee, Knoxville
confidence 95% · Zachary Johnson-Scott 1,3 ... 3 University of Tennessee, Knoxville
Ian Mulet → affiliatedwith → University of Tennessee, Knoxville
confidence 95% · Ian Mulet 3 ... 3 University of Tennessee, Knoxville
Tyler Nitzsche → affiliatedwith → University of Tennessee, Knoxville
confidence 95% · Tyler Nitzsche 3 ... 3 University of Tennessee, Knoxville
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Spiking neural networks (SNNs) offer a promising pathway to energy-efficient AI and brain-inspired computing. However, their widespread adoption is hindered by a lack of fast, accessible, and versatile simulation frameworks. In this paper, we introduce SuperNeuroMAT, an open-source, scalable, and highly efficient Python-based SNN simulator. We devise a novel matrix-based approach to model the leaky integrate-and-fire (LIF) neuron dynamics and natively support dense and sparse execution modes. This enables fast simulation of approximately 10,000 neurons in dense mode and 100,000 neurons in sparse mode on standard laptops and desktops without requiring specialized hardware. We demonstrate that SuperNeuroMAT consistently outperforms four established SNN simulators---NEST, Brian2, BindsNET, and snnTorch---on two performance metrics (execution speed and peak resident memory) and across various network sizes and connection probabilities. Furthermore, we demonstrate SuperNeuroMAT's applicability across a diverse set of problems. SuperNeuroMAT can efficiently handle conventional machine learning benchmarks such as the Digits and citation network datasets as well as neuromorphic event-based vision tasks such as N-CARS and ASL-DVS. Moreover, it can be extended beyond machine learning workloads and facilitate general-purpose workloads. We validated this by implementing the neuromorphic shortest path algorithm and two arithmetic primitives (addition and multiplication). SuperNeuroMAT can be installed via the Python Package Index (PyPI), thereby lowering the barrier to entry into the field of neuromorphic computing and accelerating the broader development of neuromorphic algorithms.
Tags
Links
- Source: https://arxiv.org/abs/2608.08479v1
- Canonical: https://arxiv.org/abs/2608.08479v1
Trouble viewing inline? Open PDF directly →
Full Text
84,810 characters extracted from source content.
Expand or collapse full text
SuperNeuroMAT: An Efficient Matrix-based Simulator for Spiking Neural Networks Prasanna Date 1* , Kevin Zhu 2 , Shruti Kulkarni 1 , Ashish Gautam 1 , Chathika Gunaratne 1 , Robert Patton 1 , Tyler Nitzsche 3 , Ian Mulet 3 , Zachary Johnson-Scott 1,3 , Addison Helms 3 , Duncan Rowden 4 , Simon Weston 4 , Maryam Parsa 2 , Catherine Schuman 3 , Thomas Potok 1 1 Oak Ridge National Laboratory, Oak Ridge, 37830, Tennessee, United States. 2 George Mason University, Fairfax, 22030, Virginia, United States. 3 University of Tennessee, Knoxville, Knoxville, 37996, Tennessee, United States. 4 Oak Ridge High School, Oak Ridge, 37830, Tennessee, United States. *Corresponding author(s). E-mail(s): datepa@ornl.gov; Abstract Spiking neural networks (SNNs) offer a promising pathway to energy-efficient AI and brain-inspired computing. However, their widespread adoption is hin- dered by a lack of fast, accessible, and versatile simulation frameworks. In this paper, we introduce SuperNeuroMAT, an open-source, scalable, and highly effi- cient Python-based SNN simulator. We devise a novel matrix-based approach to model the leaky integrate-and-fire (LIF) neuron dynamics and natively support dense and sparse execution modes. This enables fast simulation of approxi- mately 10,000 neurons in dense mode and 100,000 neurons in sparse mode on standard laptops and desktops without requiring specialized hardware. We demonstrate that SuperNeuroMAT consistently outperforms four established SNN simulators—NEST, Brian2, BindsNET, and snnTorch—on two performance metrics (execution speed and peak resident memory) and across various network sizes and connection probabilities. Furthermore, we demonstrate SuperNeuro- MAT’s applicability across a diverse set of problems. SuperNeuroMAT can efficiently handle conventional machine learning benchmarks such as the Digits and citation network datasets as well as neuromorphic event-based vision tasks such as N-CARS and ASL-DVS. Moreover, it can be extended beyond machine 1 arXiv:2608.08479v1 [cs.NE] 9 Aug 2026 learning workloads and facilitate general-purpose workloads. We validated this by implementing the neuromorphic shortest path algorithm and two arithmetic primitives (addition and multiplication). SuperNeuroMAT can be installed via the Python Package Index (PyPI), thereby lowering the barrier to entry into the field of neuromorphic computing and accelerating the broader development of neuromorphic algorithms. Keywords: Spiking Neural Network (SNN), Neuromorphic Computing, SNN Simulation, Energy-Efficient AI, Brain-Inspired Computing, Leaky Integrate and Fire (LIF) Neuron 1 Introduction Neuromorphic computing is a brain-inspired computing paradigm [1]. It performs com- putations by emulating mechanisms of the human brain [2]. Just like the human brain, which is known to consume as little as 20 W of power, neuromorphic computers are known to be extremely energy efficient. Some neuromorphic platforms, such as Intel Loihi 2, IBM North Pole, and NeuroCoreX, are known to consume thousands of times less energy than conventional hardware such as CPUs or GPUs [3–7]. The energy effi- ciency of neuromorphic computing stems from its event-driven operation, collocated processing and memory, and intrinsic parallelism [1]. Although neuromorphic com- puting is primarily explored on machine learning (ML) tasks, it is Turing-complete: that is, it is capable of executing all general-purpose computations that a CPU or a GPU can execute [8, 9]. Some examples of ML applications of neuromorphic comput- ing are graph learning [10, 11], autonomous racing [12], data filtering in high-energy physics [13], and classifying node failures on supercomputers [14]. Some examples of the non-ML applications of neuromorphic computing are graph algorithms [15–17], general-purpose computing [18–21], epidemiological simulations [22], and other scien- tific applications [23]. The primary ML model that runs on neuromorphic computers is called the spiking neural network (SNN) [24, 25]. An SNN is a type of neural network that more closely mimics the behavior and communication of the biological brain as compared to conventional artificial neural networks (ANNs). SNNs are con- sidered the third generation of neural networks and are the architectural foundation for neuromorphic computing. SNN simulators play a vital role in designing state-of-the-art SNN models. They are used to simulate computational neuroscience workloads, train SNNs off of the neuromorphic chip, co-design neuromorphic hardware through their integration with hardware simulators, and interface with domain-specific simulators (e.g., for self- driving car simulators). In the SNN literature, we observe two types of SNN simulators. On the one hand are simulators catered to neuroscience workloads such as Brian2, NEST (Neural Simulation Tool), FUGU, and STACS (Simulation Tool for Asyn- chronous Cortical Streams) [26–30]. These simulators try to solve systems of ordinary differential equations (ODEs) that govern the dynamics of biological neuron models such as the leaky integrate-and-fire (LIF) neuron model, the Izhikevich neuron model, and the Hodgkin–Huxley neuron model. Some, such as NEST, leverage discrete event 2 simulation techniques to simulate SNNs. In practice, these simulators take a long time to simulate the SNN operations. Although they work well on small problems, they soon become intractable as the number of neurons or synapses in the simulation increases. The other category of SNN simulators is those catered to ML workloads such as BindsNET, snnTorch, and spikingjelly [31–33]. These simulators are designed for speed, but they lack certain neuron or synapse parameters—for instance, synaptic delay. Moreover, neither the neuroscience-based nor the ML-based simulators are well suited for general-purpose computing tasks such as running graph algorithms, per- forming arithmetic operations for processing data, etc. Therefore, there is a need for a simulator that supports the desired functionality for neuroscience workloads and the speed required for ML workloads. Furthermore, such a simulator would have to be user-friendly, Python-based, scalable, open-source, and support a wide range of workloads, from neuroscience to ML [34]. In the work presented in this paper, we developed a fast and scalable SNN simulator called SuperNeuroMAT by leveraging a novel matrix-based approach to model the SNN computations. SuperNeuroMAT is designed for and meant to be used on personal computers such as laptops, desktops, and small-scale CPU clusters. SuperNeuroMAT can accelerate the off-chip training of SNNs on compute clusters with CPUs. It is lightweight, user friendly, Python-based, and open-source, and it can be installed using the pip command as it is available on the Python Package Installer (PyPI). The GitHub repository for SuperNeuroMAT is at the following link: https://github.com/ORNL/superneuromat. SuperNeuroMAT currently supports LIF neurons with four parameters: threshold, leak, reset state, and refractory period. Synapses in SuperNeuroMAT have weights, delays, and built-in learning based on the spike-timing=dependent plasticity (STDP) mechanism. SuperNeuroMAT supports both dense and sparse execution modes, and these modes can be selected automatically. SuperNeuroMAT further supports single precision (32-bit) and double precision (64-bit) floating-point operations on CPUs. The remainder of this paper presents the related work (Section 2), the notation (Section 3), the matrix-based approach used in SuperNeuroMAT (Section 5), scalability and benchmarking results (Section 6), and finally, the conclusion and future work (7). 2 Related Work Simulators for SNNs or neuromorphic computing can largely be divided into two cate- gories. The first category is those that are catered to neuroscience workloads, and the second is those catered towards deep learning workloads. Simulators such as NEST and Brian2 have defined the state-of-the-art in computational neuroscience [26, 27]. NEST has a Python frontend and a C++ backend, and it uses the Message Passing Interface (MPI) and is engineered for high-performance scalability. It aims to facilitate the sim- ulation of large-scale cortical structures on distributed high-performance computing (HPC) platforms. Brian2 provides users modeling flexibility through an equation-based domain-specific language. It employs runtime code generation to translate mathemat- ical differential equations into optimized C++ or Cython; thus, it enables the rapid prototyping of non-standard neural dynamics, though execution is limited to a single node. 3 Table 1: Summary of simulators for SNN and neuromorphic computing FrameworkComputational Paradigm Learning Mechanism Supported Parallelization Approach Brian2 [26]Equation-based, domain-specific language User-defined differential equations Single-Node Execution NEST [27]Discrete event sim- ulation Hebbian, plastic- ity rules Distributed MPI STACS [30]Discrete event sim- ulation Synaptic scalingCharm++ Parallel Pro- gramming Framework Fugu [28]Computation graphs Logic compositionSingle-Node Execution snnTorch [32]Frame-based, tensor-based Surrogate gradient,back- propagation through time GPU Acceleration SpikingJelly [33]Frame-based, tensor-based ANN-to-SNN, backpropagation through time CUDA/Triton Kernels BindsNET [31]Bio-inspired machine learning STDP, reinforce- ment learning PyTorch Autograd Fugu introduces a high-level algorithmic framework that represents SNNs as com- putational graphs [28]. This lets researchers compose complex spiking architectures without low-level manual tuning of synapse or neuron parameters. This abstraction layer is paired with scalable backends like STACS, which employs an asynchronous, event-driven execution model to manage high-dimensional neural communication [29]. Because these tools decouple the algorithmic specification from the execution engine, they provide a structured pathway for deploying robust spiking computations across heterogeneous computing resources, ranging from conventional CPU clusters to emerging neuromorphic hardware. The recent convergence of SNNs and deep learning has necessitated the develop- ment of frameworks that interface with automatic differentiation engines and leverage variants of the backpropagation algorithm for training SNNs. snnTorch and Spiking- Jelly are prominent PyTorch-based libraries that treat spiking neurons as recurrent primitives [32, 33]. These frameworks facilitate the training of deep SNNs via back- propagation through time (BPTT) using surrogate gradient descent to overcome the non-differentiability of the Heaviside step function (i.e., the spike function). Whereas snnTorch focuses on ease of integration with standard deep learning workflows, Spik- ingJelly provides a high-performance suite that features specialized CUDA kernels to minimize the overhead of temporal state updates. Complementing these is BindsNET, which is specifically built for biologically inspired learning paradigms, such as STDP and reinforcement learning, and provides a modular environment for unsupervised feature extraction [31]. 4 Table 2: Standard notation used to describe the LIF neuron model in the literature. SymbolDescription tTime; continuous functions of time are denoted using parentheses, (...) V (t)Membrane potential of the LIF neuron at time t V t Neuron threshold of the LIF neuron V 0 Reset membrane potential of the LIF neuron, typically around −70 mV τ m Time constant of the LIF neuron, dictates how fast the membrane poten- tial decays back to the reset state; larger means a slower leak RMembrane resistance of the LIF neuron I(t)Input current to the LIF neuron at time t t ref Refractory period of the LIF neuron Table 3: Notation used to describe the mathematical and computational approaches for neurons in SuperNeuroMAT. SymbolDescription tDiscretized time or time step; discrete functions of time are denoted using square brackets [...] NNumber of LIF neurons in the simulation 1, 2,...,NNeuron indices of all the neurons in the simulation v i [t]Membrane potential of LIF neuron i at time step t; v i [t]∈ R v[t]Vector of membrane potentials of all the neurons in the simulation at time step t; v∈ R N ν i [t]Threshold of neuron i at time step t; ν i [t]∈ R ν[t]Vector of thresholds of all neurons at time step t; ν ∈ R N λ i Leak of neuron i; SuperNeuroMAT supports a constant leak at each time step; λ i ∈ R, λ i ≥ 0 λVector of leaks of all the neurons; λ∈ R N ρ i Reset state of neuron i; ρ i ∈ R ρVector of reset states of all neurons; ρ∈ R N τ i Refractory period of neuron i; τ i ∈ W τVector of refractory periods of all neurons; τ ∈ W N τ ′ i [t]Refractory period of neuron i remaining at time step t; τ ′ i ∈ W τ ′ [t]Vector of refractory periods of all neurons remaining at time step t; τ ′ [t]∈ W N x i [t]Input spike of neuron i at time step t; x i [t] ∈ R; input spikes are given by the user and can have value x[t]Vector of input spikes of all neurons at time step t; x[t]∈ R N s i [t]Spike of neuron i at time step t; s i [t]∈ B s[t]Vector of spikes of all neurons in the simulation at time step t; s∈ B N 3 Notation We use B, N, W, and R to denote the sets of binary numbers (0, 1), natural numbers, whole numbers, and real numbers, respectively. In neuroscience, the LIF neuron is described using the notation presented in Table 2. In describing the mathematical and computational approaches pertaining to the neurons, synapses, and STDP parameters, 5 Table 4: Notation used to describe the mathematical and computational approaches for synapses in SuperNeuroMAT. SymbolDescription SNumber of synapses in the SNN (i,j)Index of synapse going from neuron i to neuron j; i,j ∈ 1, 2,...,N ω i,j [t]Synaptic weight of the synapse (i,j) at time t; w i,j [t]∈ R W [t]Synaptic weight matrix at time step t, where the value at i th row and j th column represents w ij [t]; W ∈ R N×N δ i,j Synaptic delay of synapse (i,j); δ i,j ∈ N DMatrix of synaptic delays, where the value at i th row and j th column represents δ ij [t]; W ∈ N N×N ε i,j A binary flag denoting whether STDP is enabled (1) or not enabled (0) for synapse (i,j); ε∈ B EMatrix of binary flags denoting whether STDP is enabled for all synapses; E ∈ B N×N Table 5: Notation used to describe the mathematical and computational approaches for STDP operations in SuperNeuroMAT. SymbolDescription T S Number of time steps over which STDP operations occur in the SNN; T S ∈ N α + A vector of coefficients for positive weight updates; the length of α + is T S ; α + i ∈ R ∀i∈1, 2,...,T S ; α + 1 ≥ α + 2 ≥ ...≥ α + T S ≥ 0 α − A vector of coefficients for negative weight updates; the length of α − is T S ; α − i ∈ R ∀i∈1, 2,...,T S ; α − 1 ≤ α − 2 ≤ ...≤ α − T S ≤ 0 we use the notation described in Table 3, Table 4, and Table 5 respectively. Throughout this paper, we use parentheses to denote functions of continuous variables (e.g., I(t) in Table 2) and rectangular braces to denote functions of discrete variables (e.g., v i [t] in Table 3 and Table 4). Note that we use t to denote time in both continuous and discrete domains. Whether t is continuous time or a discrete time step should be evident from the context of the equation; for instance, if an equation uses parentheses, then t in that equation is continuous time. 4 The Leaky Integrate-and-Fire Neuron Model The LIF model is one of the most widely used theoretical models to model biologi- cal neurons in computational neuroscience [35, 36]. It represents an abstraction of a biological neuron and focuses on the timing of action potentials (spikes) rather than the complex biophysics of the neuron’s ion channels. The dynamics of the LIF neu- ron can be modeled using a simple resistor–capacitor (RC) electrical circuit. The LIF dynamics can be broken down into three stages as follows: 6 1. Leak: The accumulated charge in the membrane potential of the neuron gradually leaks away and returns to its baseline resting state. This happens because the cell membrane of the neuron is not perfectly sealed; it has ion channels that allow charge to slowly leak over time. 2. Integrate: As synaptic inputs arrive from other neurons into the neuron under consideration, they deposit a charge into the neuron. This charge accumulates and increases the neuron’s membrane potential. 3. Fire: If the membrane potential continues to rise and crosses a specific threshold, then the neuron ‘fires,’ sending a spike along all of its outgoing synapses. The words ‘fire’ and ‘spike’ are used interchangeably in the literature. In this paper, ‘firing’ refers to the event that occurs when a neuron’s membrane potential reaches the threshold, whereas ‘spike’ refers to the resulting output signal. The dynamics of the LIF neuron before it fires are governed by the following linear differential equation: τ m dV (t) dt =−(V (t)− V 0 ) + RI(t)(1) The neuron spikes if the membrane potential V (t) is greater than the neuron threshold V t . In this case, the membrane potential is instantly reset to the reset state of the neuron V 0 . Although the basic LIF model is powerful, it is often modified to capture a bit more biological realism. For instance, some LIF variants introduce the notion of a refractory period such that after spiking, a neuron undergoes a brief period during which it cannot spike again. This is modeled by holding V (t) at V 0 for a fixed time interval t ref immediately after a spike. Some other variants of the LIF model use a constant leak for computational efficiency. The LIF neuron model is popular because it is computationally efficient. Other neuron models, such as Izhikevich and Hodgkin– Huxley, are more complex. Although small-scale simulations of these neuron models are possible on laptops and desktops, large-scale simulations are intractable. On the other hand, because the LIF model abstracts away ion channel dynamics, it cannot replicate complex single-neuron behaviors such as bursting and resonance or the detailed shape of an action potential, as the Izhikevich or Hodgkin–Huxly models can. 5 SuperNeuroMAT Our rationale behind designing the SuperNeuroMAT simulator was threefold: (1) Ease of use; (2) computational efficiency; and (3) open-source development. To ensure ease of use for the end user, we developed SuperNeuroMAT in Python, which is the most popular programming language in the neuromorphic community and has a vast ecosystem of libraries. Moreover, we wanted to lower the barrier to entry for new users. Therefore, we provide a very simple application programming interface (API) and a well-documented user guide. Users can install SuperNeuroMAT by running pip install superneuromat, assuming they have Python and pip installed already. We designed and engineered SuperNeuroMAT for laptops and desktops. To ensure high performance on these systems, we leveraged the numpy and scipy libraries 7 in Python [37, 38]. As a result, SuperNeuroMAT is fast and scalable on personal computers and supports a variety of workloads. It does not constrain the connec- tivity of neurons (fan-in/fan-out); the user can simulate a sparsely connected SNN with the same computational efficiency and speed as they can a densely connected SNN because SuperNeuroMAT supports both sparse and dense computations. Lastly, SuperNeuroMAT offers support for both single (32-bit) and double (64-bit) precision computations.This gives users the flexibility to achieve higher precision in their com- putation using double precision or to have more memory at their disposal and thus simulate larger SNNs using single precision. Lastly, we provide SuperNeuroMAT as a free, open-source library on GitHub to foster a collaborative research environment and software development. We envision SuperNeuroMAT being used by professors, scientists, postdocs, graduate students, undergraduate students, high school students, and SNN enthusiasts all over the world. We actively encourage contributions from the global SNN and neuromorphic computing community in further developing SuperNeuroMAT. 5.1 Neuron Parameters SuperNeuroMAT simulates LIF neurons. Each neuron has an associated index i∈ N, and the total number of neurons in the SNN is N ∈ N. On desktop and laptop computers, SuperNeuroMAT can support up to N = 10, 000 neurons with a dense implementation and around N = 100, 000 neurons with a sparse implementation. We denote the scalar versions of the neuron parameters for neuron i using Greek letters with a subscript i and the vector versions using bold Greek letters. We denote the time step in square brackets. For example, the membrane potential of neuron i at time t is denoted as v i [t]. The membrane potential of all N neurons at time t is denoted as v[t]. The membrane potential of a LIF neuron is the charge across the membrane of the neuron. In SuperNeuroMAT, the membrane potential represents the internal state of the neuron, and new charges are deposited at every time step. LIF neurons in SuperNeuroMAT have four parameters, as follows: 1. Threshold (ν): The threshold is a critical value associated with a neuron: if the membrane potential of the neuron becomes greater than the threshold, then the neuron can spike. The threshold of neuron i is denoted by ν i (ν i ∈ R ∀i), and the threshold of all neurons is denoted by the vector ν = [ν 1 ν 2 ... ν N ] T . The threshold is set when a neuron is created in SuperNeuroMAT and does not change during the course of the simulation. 2. Leak (λ): The leak is the constant amount by which the membrane potential of a neuron is either increased or decreased in order to bring it closer to the reset state of the neuron. If the membrane potential of neuron i is greater than its reset state, then the leak, denoted by λ i (λ i ∈ R ∀i), will be subtracted from its membrane potential. On the other hand, if the membrane potential of the neuron is less than its reset state, then the leak will be added to its membrane potential. The leak is added to or subtracted from the membrane potential until the membrane potential equals the reset state and no further. The leak of all N neurons in SuperNeuroMAT 8 is denoted by the vector λ = [λ 1 λ 2 ...λ N ] T (λ ∈ R N ). The leak is set when a neuron is created and does not change during the course of a simulation. 3. Reset State (ρ): The reset state of a neuron is the value that its membrane potential attains after the neuron spikes. It is also the value toward which the membrane potential of the neuron is brought by adding or subtracting leak in case the neuron does not spike. The reset state of neuron i is denoted by ρ i (ρ i ∈ R ∀i), and the reset state of all N neurons is denoted by the vector ρ = [ρ 1 ρ 2 ... ρ N ] T (ρ∈ R N ). 4. Refractory Period (τ ): The refractory period of a neuron is the number of time steps during which the neuron cannot spike after it spikes at a certain time step. For instance, if a neuron i spikes at time step t and its refractory period is τ i (τ i ∈ W ∀i), then the neuron cannot spike for the next τ i time steps. The vector of all refractory periods is denoted by τ = [τ 1 τ 2 ... τ N ] T (τ ∈ W N ). To track the refractory period remaining during the course of the simulation, we define a sister variable τ ′ = [τ ′ 1 τ ′ 2 ... τ ′ N ] such that τ ′ i ∈ W ∀i and τ ′ ∈ W N . Whereas τ does not change during the course of the simulation, τ ′ can change. In fact, τ ′ is checked and updated at every time step. A neuron i can spike at time step t only if its refractory period remaining equals 0 (i.e., τ ′ i [t] = 0). In SuperNeuroMAT, we store each neuron parameter as a Python list. When a neuron is created, we append all the neuron parameters to their corresponding lists. For instance, if a new neuron is created with threshold ν, leak λ, reset state ρ, and refractory period τ , then the lists of thresholds (ν), leaks (λ), reset states (ρ), and refractory periods (τ ) would be appended with ν, λ, ρ, and τ , respectively. 5.2 Synapse Parameters SuperNeuroMAT does not impose any restrictions on the connectivity or the fan- in/fan-out of the neurons. It can simulate all-to-all connectivity in an SNN. For an SNN with N neurons, the maximum number of synapses is N 2 , and SuperNeuroMAT indeed supports this. For computational efficiency, we provide two ways of representing the synapses in memory: dense and sparse representations. The user can select the representation that is the most apt for their application. SuperNeuroMAT also has the capability of automatically selecting dense or sparse representation based on the connectivity of the SNN. Each synapse has a unique index (i,j), where i and j are the neuron indices of the pre-synaptic and post-synaptic neurons, respectively. We support three synapse parameters as follows: 1. Weight (W): The synaptic weight represents the strength of the connection between two neurons. Each synapse multiplies the incoming spike by its weight. The weight of synapse (i,j) is denoted by ω i,j (ω i,j ∈ R). The matrix of all synaptic weights is denoted by W (W ∈ R N×N ). 2. Delay (D): The synaptic delay represents the number of time steps it takes for a spike to travel from a pre-synaptic neuron to a post-synaptic neuron. Each synapse has a default delay of 1, as we assume that it takes at least one unit of time for a synapse to propagate the spike. This assumption is fundamental to determine the computational complexity of neuromorphic algorithms [9]. The delay of synapse 9 (i,j) is denoted by δ i,j (δ i,j ∈ N). The matrix containing all synaptic delays is denoted by D (D ∈ N N×N ). 3. STDP Enabled Flag (E): We allow the option of enabling or disabling STDP learning on each synapse in SuperNeuroMAT. This is captured by the STDP enabled flag, which for synapse (i,j) is denoted by ε i,j (ε i,j ∈ B). STDP is enabled if the flag is set to 1; otherwise, it is disabled. The matrix of STDP enabled flags for all synapses is denoted by E (E ∈ B N×N ). While creating synapses, we store the synapse parameters in Python lists. When a synapse is created, we append the synapse parameters of the newly created synapse to their corresponding lists. Let us say a new synapse is to be created with weight ω, delay δ, and STDP enabled flag ε. If the delay is 1, then all synapse parameters (ω, δ, and ε) get appended to their respective lists: that is, weights list (ω), delays list (δ), and STDP enabled flags list (ε). If the delay δ is greater than 1, then SuperNeuroMAT creates a chain of additional δ− 1 neurons such that all synapses along this chain have a delay of 1. The last synapse in this chain is assigned the weight ω and STDP enabled flag ε. All synapse parameters are updated in their corresponding lists during the creation of these chained synapses. We also store the pre-synaptic and post-synaptic neuron indices for all synapses while creating them. While running the simulation, synapses are represented using dense or sparse arrays. In the dense representation, the synapse parameters are stored as O(N 2 ) numpy arrays. In the sparse representation, the synapse parameters are stored as O(S) sparse scipy arrays. 5.3 Input Spikes In SuperNeuroMAT, we enable the user to provide external input spikes to neurons at any time step. Unlike the ‘internal’ binary spikes propagated by the neurons in the SNN, the ‘external’ spikes are provided by the user and can be real valued. We denote the external spike given to neuron i at time step t by x i [t]. A vector of all input spikes at time step t is denoted by x[t]. These input spikes are ordered by the time step. For the most efficient simulation, it is better to provide details about input spikes before the simulate function is called. 5.4 Learning Parameters SuperNeuroMAT supports STDP learning with the following three parameters. 1. STDP Time Steps (T S ): This is the number of time steps over which STDP operations occur. Typically, T S ∈ W. If T S equals 0, then no STDP operations occur. 2. Coefficients for Positive Updates (α + ): This is a T S dimensional vector containing coefficients for positive weight updates. Typically, α + i ∈ R ∀i ∈ 1, 2,...,T S , α + ∈ R T S , and α + 1 ≥ α + 2 ≥ ...≥ α + T S ≥ 0. 3. Coefficients for Negative Updates (α − ): α − This is a T S dimensional vector of coefficients for negative weight updates in the STDP operation. Typically, α − i ∈ R ∀i∈1, 2,...,T S , α − ∈ R T S , and α − 1 ≤ α − 2 ≤ ...≤ α − T S ≤ 0. 10 If the user only requires positive weight updates during STDP, then α − i can be set to 0 for all i. Similarly, α + i can be set to 0 for all i if only negative weight updates are required during STDP. 5.5 Modeling the LIF Operations We model the LIF operations using a matrix-based approach. To the best of our knowledge, this approach is unprecedented in the literature. We use a matrix-based approach in order cater to the processors available on today’s laptops and desktops. The CPUs available on these machines are highly optimized to run matrix or vector operations using single instruction, multiple data (SIMD) computations. By modeling the LIF operations using matrices and vectors, and by leveraging the numpy and scipy libraries, which are based on the BLAS library in C [39], we were able to perform the LIF operations in a highly optimized and computationally efficient manner. In this section, we derive the discretized version of Equation 1 to make it amenable to computation on CPUs. Firstly, we start by multiplying Equation 1 by dt τ m on both sides. dV (t) = −(V (t)− V 0 ) + RI(t) τ m dt(2) Next, we discretize the equation in time and use square brackets to make this evident. Note that we switch to using square brackets for discrete time t according to our notation described in Section 3. ∆V [t] = −(V [t− ∆t]− V 0 ) + RI[t− ∆t] τ m ∆t(3) Next, we choose a small enough time step such that ∆t equals 1. In real physical units, this time difference can vary from application to application. Some applications may have a time step of 1 ms whereas others might have a time step of 1 ns or 1 ps. ∆V [t] = −(V [t− 1]− V 0 ) + RI[t− 1] τ m (4) Now, we move the τ m into the two terms on the right-hand side. ∆V [t] = −(V [t− 1]− V 0 ) τ m |z leak + R τ m |z weights I[t− 1] | z spikes (5) This separates the leak term from the weights and spikes term. We further split the second term on the right-hand side into two terms, one for the input spikes sent 11 by the user and the other for the spikes propagated within the network. ∆V [t] = −(V [t− 1]− V 0 ) τ m |z leak + RI in [t− 1] τ m |z input spikes + RI net [t− 1] τ m |z network spikes (6) Now, let us rewrite the above equation using the SuperNeuroMAT notation from Table 3 and Table 4. The membrane potential of neuron i at time step t in SuperNeu- roMAT is denoted using v i [t]. Therefore, the ∆V [t] term becomes v i [t]− v i [t− 1]. SuperNeuroMAT supports a constant leak at each time step. Therefore, the leak term for neuron i simply reduces to −λ i . The input spikes in SuperNeuroMAT are denoted by x i [t], which are not binary spikes but, rather, can have a real number value. There- fore, the input spikes term becomes x i [t]. Lastly, the network spikes term can be split into R τ m , which denotes the synaptic weights of all the incoming synapses and the I net [t], which denotes the binary spikes along all incoming synapses. Therefore, the above equation can be rewritten using the SuperNeuroMAT notation as follows: v i [t] − v i [t− 1] = −λ i + x i [t] + N X j=1 ω i,j [t− 1] s j [t− 1](7) Note that the summation of j is taken over all the neurons in the SNN. If a synapse from neuron j to neuron i does not exist in the SNN, then we set its weight to 0. From a computational efficiency point of view, it is more efficient to perform this operation over weights that are initialized to 0 rather than searching for only those synapses that exist in the SNN. We now move the v i [t− 1] term to the right-hand side to yield the equation to update the membrane potential of neuron i at time step t: v i [t] = v i [t− 1] − λ i + x i [t] + N X j=1 ω i,j [t− 1] s j [t− 1](8) We now vectorize this equation to compute the membrane potential of all N neurons in one shot: v[t] = v[t− 1] − λ + x[t] + W T [t− 1] s[t− 1](9) Note that in the matrix W , each element w i,j denotes the synaptic weight of the synapse going from neuron i to neuron j. In order to match the weights to the incoming spike vector s, we must take the transpose of W . Equation 9 gives us a vectorized equation to update the membrane potential of all neurons at time step t. We would like to point out a few considerations about computing the leak in Equation 9. The leak tries to bring the membrane potential of the neuron back to the rest state. If the membrane potential v i [t− 1] is greater than the reset states ρ i for neuron i, then the leak is subtracted as shown in Equation 1. However, if v i [t− 1] < ρ i for some neuron i, then the constant leak λ i will be added to the membrane potential 12 instead of being subtracted in order to bring the membrane potential closer to the reset state. Another point to keep in mind about the leak is that after subtracting/adding leak to the membrane potential, if the membrane potential becomes lower/higher than the reset state, then we set the membrane potential to the reset state. The neuron spikes if its membrane potential is greater than its threshold and its refractory period remaining at the current time step is 0. If the neuron spikes, then it enters its refractory period. We do this by setting its refractory period remaining τ ′ i to its refractory period τ i . If the neuron does not spike, we decrement its refractory period by 1 until it reaches 0. Accordingly, we have the two following equations to compute the spikes and the refractory periods. s i [t] = ( 1if v i [t] > ν i and τ ′ i [t] = 0 0otherwise (10) τ ′ i [t] = ( τ i if s i [t] = 1 max0, τ ′ i [t− 1]− 1otherwise (11) where the max... function selects the maximum element out of all the elements in a given set. 5.6 Modeling the STDP Operations The STDP algorithm is often colloquially summarized as follows: neurons that fire together, wire together. More formally, the STDP algorithm can be described as follows. Let us say that a neuron i spiked at time step t− 1 and another neuron j spiked at time step t (i.e., neuron j spiked immediately after neuron i). In this case, the weight of synapse (i,j) should be increased. The rationale behind this is that the spiking behavior of neuron i is positively correlated with that of neuron j. In some cases, neuron i could even be thought of as causing neuron j to spike, although causality is not always guaranteed. Conversely, if neuron j spiked at time step t− 1 and neuron i spiked at time step t—that is, if neuron j spiked before neuron i—then the weight of synapse (i,j) will be decreased. This is because neuron j spiked before neuron i, meaning their spiking behavior is negatively correlated. Now, in a series of steps, we will formulate the STDP weight update as matrix– vector operations. This is how we implement STDP operations in SuperNeuroMAT. First, we want to find the indices of the weights that need to be updated. To do so, we need to look at the spike vectors at times t− 1 and t. Since the spike vectors are binary (s[t− 1], s[t]∈ B N ), their outer product gives us an N×N binary matrix. The non-zero entries of this matrix represent the weights for which the pre-synaptic neuron spiked at time step t− 1 and the post-synaptic neuron spiked at time step t. This is a positive correlation, as described above, and these weights must be increased. Thus, the positive weight update matrix at time step t, ∆W + [t] can be written as follows. ∆W + [t] = s[t− 1] s T [t](12) 13 Note that if we subtract ∆W + [t] from an N × N matrix of ones, we will get the indices of the weights that should be decreased. Thus, the negative weight update matrix at time step t, ∆W − [t] can be written as follows. ∆W − [t] = 1 N×N − s[t− 1] s T [t](13) Typically, ∆W + and ∆W − are multiplied by coefficients α + 1 and α − 1 , respectively. These coefficients are called potentiation amplitude and depression amplitude, respec- tively, in neuroscience. They are described in Table 5. This is akin to the learning rate used in the backpropagation algorithm for training ANNs. Accounting for these coefficients, we can rewrite ∆W + [t] and ∆W − [t] as follows. ∆W + [t] = α + 1 s[t− 1] s T [t](14) ∆W − [t] = α − 1 1 N×N − s[t− 1] s T [t] (15) Furthermore, STDP operations occur over T S time steps. Accordingly, we have T S amplitudes of potentiation and depression, α + 1 ,α + 2 ,...,α + T S and α − 1 ,α − 2 ,...,α − T S . These are applied to the outer products of the spike vectors over previous T S time steps. The resulting matrices are summed together to yield an updated expression for the positive and negative weight updates as follows. ∆W + [t] = T S X t s =1 α + t s s[t− t s ] s T [t](16) ∆W − [t] = T S X t s =1 α − t s 1 N×N − s[t− t s ] s T [t] (17) Essentially, for positive weight updates, we apply the first potentiation amplitude α + 1 to the most recent spike vector, s[t− 1], the second potentiation amplitude α + 2 to the second most recent spike vector, s[t− 2], and so on. Typically, we have α + 1 ≥ α + 2 ≥ ...≥ α + T S ≥ 0, giving more importance to the recent spike vectors and less importance to older spike vectors. We eventually add all these matrices, and that becomes our positive weight update matrix at time step t. The same can be said for the negative weight updates, where we have α − 1 ≤ α − 2 ≤ ...≤ α − T S ≤ 0. Given W + [t] and W − [t], we can write the overall weight update ∆W as follows. Note that we do not need to subtract ∆W − because the depression coefficients are already assumed to be negative. ∆W [t] = ∆W + [t] + ∆W − [t](18) Next, SuperNeuroMAT allows users to enable STDP on certain synapses and dis- able it on others. We have a variable ε i,j associated with synapse (i,j) to capture this information. If STDP is enabled on synapse (i,j), then ε i,j is set to 1, else 0. This information is also represented in the matrix E, where the element in the i th row and j th column is ε i,j . E is called the STDP enabled matrix, and E ∈ B N×N . To allow 14 STDP operations on only those synapses for which STDP is enabled, we need to multi- ply the weight updates by E element-wise. We denote the element-wise multiplication operator for two matrices using the ⊙ symbol. The resulting matrix is added to the weight matrix W to update the weights at time step t. Therefore, the overall weight update can be written as follows. W [t] = W [t− 1] + E⊙ ∆W [t](19) 5.7 Setting Up the Simulation Environment All the neuron parameters, synapse parameters, and learning parameters are stored in memory as Python lists when the SNN is being created. The input spikes are stored as a Python dictionary, keyed on the time steps for efficient retrieval during simulation. When the simulation environment is being initialized, numpy arrays are created for all neuron and learning parameters. If the representation of matrices pertaining to the synaptic parameters is dense, then the synaptic weights and STDP enabled flags are stored as numpy arrays. If the representation is sparse, then they are stored as scipy sparse arrays. We also initialize variables to store the spikes and the spike trains. 5.8 Simulation Algorithm The simulation algorithm implements Equation 9, Equation 10, Equation 11, and Equation 19 at each time step. We leverage numpy and scipy arrays to store the vectors and matrices in memory. We also implement the computations in a manner that is amenable for SIMD computations on hardware. We assume here that the simulation is run for T time steps and determine the computational complexity of updating the membrane potential of the neurons (Equation 9) and updating the synaptic weights using STDP (Equation 19). These two equations are the dominant computations that must be performed at each simulation time step. For storing the vectors in Equation 9—v, λ, x, s etc.—we need O(N ) space. For storing the matrix W under a dense implementation, we need O(N 2 ) space because with N neurons, the maximum number of synapses possible is N 2 . Under a sparse implementation, storing W would requireO(S) space. Therefore, the space complexity in a dense implementation isO(N 2 ), and that in a sparse implementation isO(N +S). To compute the vector operations, we would needO(N ) time. To compute the matrix operations under a dense implementation, we would need O(N 2 ) time, and under a sparse implementation, we would need O(S) time. Therefore, the overall time com- plexity under a dense implementation is O(N 2 ), and that in a sparse implementation is O(N + S). The space complexity for computing Equation 19 is O(N 2 ) for a dense implemen- tation andO(S) for a sparse implementation. Since all the vectors are stored as dense vectors in SuperNeuroMAT, the time complexity of computing the outer products in Equation 16 and Equation 17 is O(N 2 ). We compute T S such outer products at each time step. Therefore, the overall time complexity of computing Equation 16 and Equation 17 is O(N 2 · T S ). The time complexity of computing Equation 19 is O(N 2 ) 15 for a dense implementation andO(S) for a sparse implementation. However, the over- all time complexity of STDP computations is dominated by Equation 16 and Equation 17. Therefore, the overall time complexity of STDP computations is O(N 2 · T S ), and the overall space complexity of STDP computations is O(N 2 ) for dense and O(S) for sparse implementations. We will also need O(N · T ) space to store the spike train. Therefore, the overall space complexity for dense computations is O(N 2 ) +O(N 2 · T S ) +O(N · T ). The overall space complexity for sparse computations is O(S) +O(N 2 · T S ) +O(N · T ). Taking the highest-order terms, we get overall space complexity asO(N 2 ·T S +N·T ). Similarly, the overall time complexity for each time step is O(N 2 ·T S ) for both sparse and dense computations. This is dominated by the STDP operations. The overall time complexity across T simulation time steps is O(N 2 · T S · T ). Note that since all operations are modeled using matrices and vectors, they can be parallelized to take advantage of multiple CPUs on a compute platform. Numpy and scipy inherently support this parallel execution because they are based on the BLAS library in C. Although the time complexity is a polynomial expression in the number of inputs and is thus considered efficient, we see even better performance in practice because of the parallel execution. 6 Results SuperNeuroMAT supports both sparse and dense representations of the SNN. We con- ducted a detailed study comparing the total run time and memory consumption for both sparse and dense representations of various SNNs generated uniformly at random, as detailed in in Section 6.1. We also compared the performance of SuperNeuroMAT with four SNN simulators: NEST, Brian2, BindsNET, and snnTorch. As shown in Section 6.2, we used SuperNeuroMAT in two conventional benchmark ML examples (Digits and Citation Networks), two neuromorphic ML examples (N-CARS and ASL- DVS), and two non-ML examples (shortest path and arithmetic primitives). These problems are crucial for general-purpose computing applications of neuromorphic computing. 6.1 Computational Efficiency All experiments in this section were conducted on a KVM virtual machine running Rocky Linux 9.8 with 32 vCPUs and 128 GiB of RAM, hosted on an AMD EPYC 7702 64-core (Zen 2) processor. Simultaneous multithreading was disabled and a single NUMA node was used. Storage consisted of a 60 GiB virtio block device holding the root file system as well as a 1.2 TiB shared NFS v4.2 volume mounted for user data. Transparent huge pages were enabled in always mode, and no swap was configured, so no measurement was affected by paging. 6.1.1 Performance of Sparse and Dense Simulation Modes We first compare the performance of the sparse and dense execution modes in SuperNeuroMAT. In this section, the connection probability is defined as S/N 2 , where 16 (a) 10 Neurons(b) 100 Neurons (c) 1,000 Neurons(d) 10,000 Neurons Fig. 1: Comparison of run times of sparse and dense modes in SuperNeuroMAT for (a) 10 neurons; (b) 100 neurons; (c) 1,000 neurons; and (d) 10,000 neurons. The connection probability is shown on the X-axis, and the total run time in seconds is shown on the Y-axis. S is the number of synapses, and N is the number of neurons in the SNN. Since N 2 is the maximum possible number of synapses in an SNN of N neurons, this quantity expresses the fraction of connections that are actually present. We evaluated a struc- tured grid of four network sizes (N = 10, 100, 1,000, and 10,000) and five connection probabilities ( S N 2 = 10 −4 , 10 −3 , 10 −2 , 10 −1 , and 1.0). Each of these twenty combina- tions was simulated in both modes, yielding forty distinct simulation configurations, each of which was executed ten times for a total of four hundred runs. The values shown in Figure 1 and Figure 2 are averaged over the ten runs of each configuration. Each SNN was constructed by first creating N neurons and then creating synapses uniformly at random according to the connection probability of the configuration. STDP was enabled on all synapses, and input spikes were assigned uniformly at ran- dom. Every network was simulated for 100 time steps. For each run, we recorded the network setup time, the simulation time, the total run time, and the peak resident 17 (a) 10 Neurons(b) 100 Neurons (c) 1,000 Neurons(d) 10,000 Neurons Fig. 2: Comparison of memory consumed by the sparse and dense modes in SuperNeu- roMAT for (a) 10 neurons; (b) 100 neurons; (c) 1,000 neurons; and (d) 10,000 neurons. The connection probability is shown on the X-axis, and the total memory consumed in megabytes (MB) is shown on the Y-axis. memory so that computational and memory costs could be compared directly across the two modes. Figure 1 presents the total run time (Y-axis) for both sparse and dense modes as a function of connection probability (X-axis) for each network size; because the Y- axis denotes run time, lower values are better. Figure 2 presents the peak memory consumption (Y-axis) for both modes as a function of connection probability (X-axis) for each of the four network sizes; once again, lower values are better. The results reveal a clear crossover the location of which depends jointly on network size and connectivity. For the smallest networks (10–100 neurons), dense mode was consistently faster at every connection probability, with simulation times of roughly 0.004–0.011 s compared with 0.07–0.11 s in sparse mode. Peak memory in this regime was essentially identical between the two modes (approximately 80 MB) and was dominated by run time overhead rather than by the network itself. At 1,000 neurons the ordering reversed at low and moderate connectivity: mean total run time in sparse mode was 5.1× faster 18 at 10 −4 (0.166 s vs. 0.851 s), 3.6× faster at 10 −3 (0.226 s vs. 0.821 s), and 2.4× faster at 10 −2 (0.400 s vs. 0.949 s), with correspondingly lower peak memory (84–88 MB vs. 119–122 MB). At a connection probability of 10 −1 , the two modes were effectively tied in terms of run time (1.72 s vs. 1.75 s), although the sparse mode retained a memory advantage (about 124 MB vs 155 MB), and at full connectivity the dense mode was again faster (3.23 s vs 9.82 s). The sparse mode performance was the most prominent at the largest scale. For 10,000 neurons, sparse mode reduced total run time by approximately 45× at a con- nection probability of 10 −4 (3.93 s vs. 179 s), 26× at 10 −3 (6.67 s vs. 176 s), and 8× at 10 −2 (23.0 s vs. 184 s), and it remained faster at 10 −1 (200 s vs. 273 s). Peak memory usage in these regimes was between 1.3× and 14× lower than in dense mode, with the largest reductions at the sparsest connectivities. At full connectivity, however, sparse mode required substantially more time (2,140 s vs. 508 s) and slightly more memory (27.2 GB vs. 25.6 GB) than the dense mode. The results above indicate that neither representation dominates across the param- eter space. Dense execution is preferable when the problem is small because its lower per-operation overhead outweighs any redundant computation, as well as when the network is close to fully connected because the sparse representation traverses almost as many elements as the dense one while paying additional indexing costs. Sparse execution is preferable for large, genuinely sparse networks, which is the regime most representative of application-scale SNNs; it yielded speedups of up to 45× and memory reductions of up to 14× in our experiments. Because users can present SuperNeuro- MAT with any of these configurations, it supports both execution modes. When the user does not explicitly select a mode, SuperNeuroMAT chooses one automatically from the number of neurons and the number of synapses in the SNN. Dense mode is selected for small problems and for large, densely connected problems. Sparse mode is selected for large, sparse problems. 6.1.2 Comparison with Other Simulators We now evaluate the computational efficiency of both sparse and dense modes of SuperNeuroMAT against four widely used SNN simulation frameworks: NEST, Brian2, BindsNET, and snnTorch. Rather than relying on a specific application-derived net- work, the comparison uses synthetic Erd ̋os–R ́enyi random graphs spanning four network sizes (10, 100, 1,000, and 10,000 neurons) and six connection probabilities (0.01, 0.1, 0.25, 0.5, 0.75, and 1.0), with 10 independently generated graphs for each combination of network size and connection probability. Each graph is converted into an equivalent SNN in every framework, and all neurons receive external spikes at every fifth time step. While each simulator operates differently, we ensured that the neuron and synapse models are matched as closely as possible within what is allowed under each framework’s API. In each simulator, we ensured that the neuron model being sim- ulated matches the default configuration of SuperNeuroMAT’s neuron and synapse: zero threshold, memoryless (infinite-leak) LIF neurons with instantaneous (“delta”) synapses carrying a unit delay. This was implemented via NEST’s iaf pscdelta, a custom-scheduled Brian2 equation set, BindsNET’s LIFNodes with near-zero decay, 19 (a) 10 Neurons(b) 100 Neurons (c) 1,000 Neurons(d) 10,000 Neurons Fig. 3: Run time comparison for SuperNeuroMAT sparse (green), SuperNeuroMAT dense (blue), NEST (orange), Brian2 (yellow), BindsNet (red), and snnTorch (purple) for (a) 10 neurons; (b) 100 neurons; (c) 1,000 neurons; and (d) 10,000 neurons. The connection probability is shown on the X-axis, and the time in seconds is shown on the Y-axis. and snnTorch’s Leaky neuron with zero decay. By equalizing simulation semantics across frameworks, differences in measured runtime and memory can be attributed to implementation efficiency rather than to differences in model complexity. For each run, setup time (data-structure preparation) and simulation time (100 time steps of execution) were recorded separately, along with peak resident memory usage, and 20 (a) 10 Neurons(b) 100 Neurons (c) 1,000 Neurons(d) 10,000 Neurons Fig. 4: Peak memory comparison for SuperNeuroMAT sparse (green), SuperNeuro- MAT dense (blue), NEST (orange), Brian2 (yellow), BindsNet (red), and snnTorch (purple) for (a) 10 neurons; (b) 100 neurons; (c) 1,000 neurons; and (d) 10,000 neu- rons. The connection probability is shown on the X-axis, and the time in seconds is shown on the Y-axis. SuperNeuroMAT was evaluated in both its dense (NumPy) and sparse (SciPy sparse matrix) modes. Figure 3 shows the comparison of run times in seconds along the Y-axis as a function of increasing connection probability along the X-axis for four network sizes. Across the tested range, SuperNeuroMAT outperforms all four simulators, and its performance in relation to other simulators improves significantly as the network size increases. At 1,000 neurons (averaged across all connection probabilities), SuperNeu- roMAT’s dense mode completes in 0.052 s, versus 0.56 s for Brian2, 2.52 s for snnTorch, 2.81 s for BindsNET, and 7.05 s for NEST—a 10–135× speedup. At 10,000 21 neurons under full connectivity (the most computationally demanding configuration tested), SuperNeuroMAT dense mode requires 7.18 s, compared with 70.0 s for Brian2 (9.8×), 578.7 s for snnTorch (81×), 613.3 s for BindsNET (85×), and 1,545.6 s for NEST (215×). Notably, at the smallest network size (10 neurons) NEST is marginally faster than SuperNeuroMAT, reflecting its low fixed per-run overhead. However, this crossover disappears by 100 neurons, since NEST’s per-neuron and per-edge object- creation calls scale poorly with network size, and SuperNeuroMAT’s advantage grows monotonically thereafter. Figure 4 shows the comparison of peak resident memory in MB along the Y-axis as a function of increasing connection probability along the X-axis for four network sizes. SuperNeuroMAT maintains a consistently smaller memory footprint. At small scale (10–100 neurons), its peak memory usage is negligible (well under 1 MB), whereas the other frameworks carry substantial fixed overhead—largely from their underly- ing runtime systems (PyTorch for BindsNET and snnTorch, the NEST kernel, and Brian2’s code-generation backend)—ranging from roughly 130 MB (NEST) to over 800 MB (BindsNET). At the largest scale tested (10,000 neurons, full connectivity), SuperNeuroMAT consumes approximately the same peak memory as Brian2 and is slightly higher than the remaining three simulators. However, the 2–3 orders of mag- nitude speedup obtained in the run times shown in Figure 3 by SuperNeuroMAT as compared to NEST, BindsNET, and snnTorch, more than compensates for the slightly higher peak memory consumption. Note that the memory values presented in Figure 4 include the memory consumed by the Erd ̋os–R ́enyi graphs, which were created using the networkx library in Python. 6.2 Benchmark Examples We provide six examples of how SuperNeuroMAT can be used in a variety of different settings. We focus on two conventional machine learning examples: the Digits dataset and citation networks datasets (Cora, Citeseer, and Pubmed). Next, we delve into two neuromorphic ML datasets: the N-CARS dataset and the ASL-DVS dataset. Finally, we look into two non-ML examples: the shortest path algorithm and the arithmetic primitives. 6.2.1 Digits We use SuperNeuroMAT to implement a two-layer SNN for the Digits dataset, which is available in the Scikit-learn library in Python. The Digits dataset is a benchmark handwritten digit classification dataset containing 1,797 grayscale images of decimal digits (0–9), with one class label per image. Each image has a spatial size 8× 8, so every sample is represented by 64 pixel features after flattening. Pixel intensities are quantized integer values in the range [0, 15], corresponding to 4-bit grayscale levels. We split the data into train and test subsets using a 80/20 partition without shuffling, yielding approximately 1, 438 training samples and 359 test samples. Our SNN contained 64 input neurons for the 8×8 pixels and 10 output neurons for the ten digit classes. We encoded the pixel intensities as spike amplitudes and 22 Fig. 5: Learned weights for the Digits dataset using a two-layer SNN in SuperNeuro- MAT. Fig. 6: Confusion matrix for the Digits dataset using a two-layer SNN in SuperNeu- roMAT. added input spikes to the correct output neuron during training to facilitate STDP- based synaptic updates. Figure 5 shows the learned synaptic weights after the STDP operations. After training, the learned weight matrix was used for inference with STDP disabled. Predictions were decoded by selecting outputs with the highest spike counts and then resolving ties using earliest spike times. We obtained 93% one-versus-all accuracy, which is at par with the state-of-the-art for such a two-layer SNN. Figure 6 shows the confusion matrix for this dataset. The Digits example demonstrates that a 23 compact two-layer SNN can be successfully trained in SuperNeuroMAT using built-in STDP. 6.2.2 Citation Networks We used SuperNeuroMAT for node classification on three citation network datasets (Cora, Citeseer, and Pubmed). These are well-established benchmarks for evaluating graph learning algorithms. In these datasets, each node corresponds to a scientific pub- lication, directed edges denote citation relationships, and every publication belongs to one of several research topics [40, 41]. Conventional approaches formulate this problem using graph neural networks (GNNs), where node attributes are iteratively aggregated through neighborhood message passing and optimized using backpropaga- tion [42]. In contrast, we develop a neuromorphic formulation in SuperNeuroMAT such that publications are represented as neurons and citation relationships as synapses. This approach enables classification to emerge from event-driven spike propagation and local synaptic plasticity. This approach relies exclusively on graph topology and does not utilize node features, demonstrating that meaningful classification can be achieved through biologically plausible local learning operating directly on the graph structure. Cora consists of 2,708 publications connected by 5,429 citation links and catego- rized into seven research topics. CiteSeer contains 3,327 publications, 4,732 citation links, and six topic classes, while PubMed comprises 19,717 publications connected by 44,338 citation links across three research areas. Standard Planetoid train/valida- tion/test splits were adopted for all experiments [43]. Each citation graph was mapped directly onto an SNN, with every publication represented by a neuron and every cita- tion represented by a synapse. To suppress multiple spikes at a given neuron arising from cycles within the citation graph, neurons were assigned a sufficiently long refrac- tory period such that each neuron emitted at most one spike during a simulation. A dedicated output neuron was allocated to each topic class, and every publication neuron in the test set was connected to all topic neurons through trainable synapses initialized with small random weights. During classification, a spike injected into a test publication neuron propagated through the citation graph according to its connectiv- ity. Synaptic weights between publication neurons and topic neurons were updated online using STDP learning. At the end of the simulation, the predicted class was assigned to the topic neuron whose incoming synapse from the test publication neuron attained the largest weight. Table 6 summarizes the classification performance—these are state-of-the-art results for spiking graph neural network models that do not incorporate features. Although the proposed method does not employ node features, global message passing, or gradient-based optimization, it successfully extracts discriminative infor- mation directly from graph topology through local spike-based computation. Our results demonstrate that structurally meaningful representations can emerge from purely local synaptic adaptation. This highlights the flexibility of SuperNeuro- MAT for graph-structured learning tasks beyond conventional neuromorphic sensory benchmarks. 24 Table 6: Node classification performance on citation net- work benchmarks using graph topology alone. No node features were used during learning or classification. DatasetNodesCitation LinksClassesAccuracy (%) Cora2,7085,429762.7 CiteSeer3,3274,732642.4 PubMed19,71744,338354.7 Fig. 7: Mean spike activity of the excitatory neurons associated with false positive and false negative classifications across simulation time steps. 6.2.3 N-CARS Approach 1: Diehl–Cook-Inspired STDP with Neuron-Level Fine-Tuning We used SuperNeuroMAT to evaluate unsupervised SNN approaches for car-versus- background classification on the event-based N-CARS dataset [44]. The first approach employed a Diehl–Cook-inspired architecture comprising an input layer with 1,080 input neurons to encode the event camera voxel grid, a 64-neuron excitatory population, and a 64-neuron inhibitory population implementing winner-take-all com- petition [45]. Input-to-excitatory synapses were trained without class labels using STDP, after which each excitatory neuron was assigned to the class that produced its highest average training set spike response. SuperNeuroMAT’s spike train was then used to compare activity associated with true positive, true negative, false positive, and false negative predictions. As shown in Figure 7, two neurons exhibited firing rates at least 30% higher during misclas- sifications, whereas six additional neurons displayed slowly decaying error-associated activity. The two neurons were therefore deactivated, and the leak values of the six remaining neurons were increased from 0.05 to 0.15. These targeted interventions improved the excitatory neuron accuracy by 4.3 percentage points, from 60.7% to 65.0%. This example demonstrates how SuperNeuroMAT can be used to interpret and enhance SNNs in an unsupervised learning setting. 25 (a) Population-level neural activity.(b) Logistic regression classification results. Fig. 8: Analysis of the biologically inspired SuperNeuroMAT representation on the N-CARS dataset. Approach 2: Biologically Inspired Neuromorphic Feature Extraction We used SuperNeuroMAT in a custom biologically inspired feature extraction pipeline to perform car-versus-background classification on the event-based N-CARS dataset [44]. Unlike conventional approaches that rely on established event represen- tations, such as Histograms of Averaged Time Surfaces (HATS) [44], Hierarchies of Time Surfaces (HOTS) [46], or related time-surface methods, the proposed approach learned a neuromorphic representation directly from the asynchronous event stream. Events from the 100× 120 pixel sensor were spatially mapped to a 10× 12 lattice of ON- and OFF-polarity relay neurons and processed by a sparse 960-neuron SNN. The network incorporated center–surround filtering, fast and slow LIF traces, lateral inhi- bition, cross-polarity inhibition, and STDP at the excitatory trace-to-output synapses. The use of competitive lateral inhibition and unsupervised STDP builds on estab- lished SNN feature learning principles, including the competitive learning architecture introduced by Diehl and Cook [45]. STDP training was performed in an unsupervised learning setting using the 80% training dataset. After training, the SNN was frozen, and features derived from output spike counts, first-spike latencies, and population-level activity (Figure 8a) were com- bined into a 2,880-dimensional representation. This representation was then evaluated using both a regularized logistic regression classifier and a nonlinear multilayer percep- tron with a 2,880→ 128→ 64→ 2 architecture. On the untouched 20% test partition, comprising 4,806 recordings, the regularized logistic regression achieved 87.7% accu- racy, an F1 score of 88.0%, and a receiver operating characteristic area under the curve (ROC-AUC) of 0.938 (Figure 8b). When the same SuperNeuroMAT representation was paired with the nonlinear MLP, performance improved to 92.4% accuracy, an F1 score of 92.7%, and a ROC-AUC of 0.975, demonstrating that the learned neuromorphic features were particularly effective when combined with a nonlinear classifier. 6.2.4 ASL-DVS We further evaluated the proposed SuperNeuroMAT feature-extraction approach on the ASL-DVS dataset [47]. ASL-DVS contains 100,800 event camera recordings repre- senting 24 American Sign Language letters from A–Y, excluding J and Z. Each class contains 4,200 recordings captured at a spatial resolution of 240× 180 pixels. An illus- trative image is shown in Figure 9. A stratified 80/20 partition was used in the present 26 Fig. 9: ASL-DVS letter classification results. study, producing 80,640 training recordings and 20,160 held-out test recordings, with 840 test samples per class. SuperNeuroMAT was used to construct and simulate the SNN and perform STDP learning. Events were mapped to separate ON- and OFF-polarity pathways within a 10× 12 receptive field lattice and processed using a 2 ms simulation time step. The resulting 960-neuron network incorporated center–surround filtering, fast- and slow- memory LIF neurons, lateral inhibition, and cross-polarity competition. STDP was restricted to the 480 excitatory trace-to-output synapses, which were trained without class labels using only the training partition. Following unsupervised STDP training, the network parameters were frozen, and each recording was transformed into a 1,548-dimensional feature representation composed of output spike counts, first-spike latencies, and population-level activity measurements. The resulting features were standardized and supplied to a regular- ized logistic regression classifier. On the test set, the SNN–STDP pipeline achieved 85.3% accuracy, with a 95% confidence interval of 84.8%–85.7%. Macro-averaged pre- cision, recall, and F1 score were 85.3%, 85.3%, and 85.2%, respectively, while the multiclass receiver operating characteristic area under the curve (ROC-AUC) reached 0.992. Classification performance was strongest for the G, P, Q, and H signs, whereas M, S, and T were among the most difficult classes to distinguish. 6.2.5 Shortest Path Algorithm To evaluate SuperNeuroMAT’s ability to solve graph problems natively through spike timing, we implemented the neuromorphic shortest path algorithm [15] for a directed, weighted graph of 8 nodes and 15 edges as shown in Figure 10 [48]. Each node in the graph is represented by a neuron with a threshold of zero and a large refractory period such that it spikes exactly once during the course of the simulation. Each directed edge is represented by a synapse whose delay is proportional to the corresponding edge weight scaled by a factor of 100 to preserve two decimal digits of precision. Under this encoding, the arrival time of a spike at a given neuron corresponds to the cumulative weight of the path it traveled. A single input spike was injected into a given source node at time zero, and the network was simulated for 300 time steps. This was repeated for all eight source nodes. 27 Fig. 10: The graph containing eight nodes and fifteen edges used for comparing the neuromorphic shortest path algorithm implemented in SuperNeuroMAT with Dijk- stra’s algorithm. The spike times recorded for each neuron were decoded back into path distances by dividing by 100. The resulting values were compared against shortest distances com- puted using Dijkstra’s algorithm. For all pairs of source and destination nodes, the neuromorphic shortest distances matched those from Dijkstra’s algorithm. These are shown in Table 7. For instance, for the source neuron 0, neuron 2 spiked at time step 26, corresponding to a distance of 0.26, and neuron 6 spiked at time step 151, correspond- ing to a distance of 1.51. These results demonstrate the usability of SuperNeuroMAT in solving the shortest path problem without requiring explicit priority queues or relaxation steps as required in Dijkstra’s algorithm. 6.2.6 Arithmetic Primitives To perform arithmetic operations using SNNs in SuperNeuroMAT, we first encode each operand as a fixed-width binary integer, an n-bit integer across n input neurons: neuron i emits a spike if bit i is 1. For addition, bit i of each operand feeds a triad of neurons with thresholds 0, 1, and 2 that reproduce the sum and carry logic of 28 Table 7: Shortest paths obtained from each source node to each destination node. These were obtained by implementing the neuromorphic shortest path algo- rithm in SuperNeuroMAT. When scaled down by 100, these shortest paths match the shortest paths obtained using the Dijkstra’s shortest path algorithm. SourceDestination 01234567 001052699387315160 113701202817318280154 218294073976212534 310918692014515452126 4185671687603512837 5169321526035011228 6571344011393102074 714860131396328910 Fig. 11: 4-bit spiking adder. Takes two 4-bit operands (A and B) encoded using binary encoding as inputs, performs the addition operation through the hidden neurons (H) and computes a 5-bit output (C). a full adder; the carry neuron at position i is connected to the triad at position 29 i + 1 with a one time step delay, so that as the input spike wavefront propagates through the network, the carry ripples forward automatically. This results in a static, feed-forward ripple-carry adder whose control logic is the pattern of synaptic delays. Figure 11 shows a 4-bit adder as an illustration of this approach. Because the SNN is purely feed-forward, its single instantiation can be pipelined such that successive operand pairs are injected one time step apart and are processed concurrently as they move through the network. We validated our approach by building an 8-bit ripple- carry adder. Including the delay neurons, this adder SNN contained 293 neurons and 336 synapses. We streamed 100 randomly generated pairs of 8-bit unsigned integers through it and decoded the resulting spike trains back into decimal values. Comparing against conventional integer addition showed exact match for all 100 input pairs. We also implemented an 8-bit multiplication SNN. Each partial product A i ∧B j was computed with a dedicated spiking neuron. The resulting eight shifted partial product rows were summed using a binary merge tree of delay-based ripple-carry adders rather than a single sequential accumulation: at each step, the two partial sums that are ready earliest are combined first, so most merges require no additional resynchronization delay and the SNN depth scales with log 2 (N ), where N is the number of bits. Testing this SNN on 20 randomly generated pairs of 8-bit unsigned integers, the correct 16-bit product was recovered from the output spike times in all 20 cases, using 2,238 neurons and 2,802 synapses. 7 Conclusion This paper introduced SuperNeuroMAT, a fast, scalable, and open-source SNN simu- lator that leverages a novel matrix-based approach to model the LIF operations, our simulator achieves exceptional computational efficiency on standard CPUs. Both dense and sparse execution modes efficiently simulate up to 10,000 neurons under a dense implementation and approximately 100,000 neurons under a sparse implementation directly on personal computers. Our results indicate that SuperNeuroMAT consistently outperforms four estab- lished SNN frameworks—NEST, Brian2, BindsNET, and snnTorch—in both total run time and peak resident memory across a wide range of connection probabilities and network sizes. Furthermore, we showcased the versatility of SuperNeuroMAT across six diverse examples. We applied SuperNeuroMAT to conventional ML benchmarks (Digits and citation network datasets) as well as to neuromorphic event-based vision datasets (N-CARS and ASL-DVS). Finally, we validated its capability for non-ML, general-purpose computing by accurately implementing a neuromorphic shortest path algorithm and pipelined arithmetic primitives. Our guiding rationale was to create a framework that is both highly performant and widely accessible; therefore, we developed it to be Python-based, user-friendly, and easily installable via PyPI. By lowering the barrier to entry and providing a powerful tool for SNN simulation and off-chip training, we intend for SuperNeuroMAT to foster a collaborative global research environment and accelerate the broader development of SNNs and neuromorphic computing. 30 In the future, we would like to implement SuperNeuroMAT in a distributed fashion to simulate billions of neurons and trillions of synapses. Such an implemen- tation will have to leverage both CPUs and GPUs. Unlike desktops or laptops, this implementation would cater to exascale HPC systems such as Frontier, Lux, and Discovery. Acknowledgements This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, under contract number DE-AC05-00OR22725. This manuscript has been co-authored by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725 with the US Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world- wide license to publish or reproduce the published form of this manuscript or allow others to do so, for United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy. gov/downloads/doe-public-access-plan). The authors would like to thank John Baston I, Technical Editor at the Oak Ridge National Laboratory (ORNL), for reviewing and editing this manuscript. The authors would like to thank Nicholas Quentin Haas, Software Developer at ORNL, for providing access to the virtual machine on which the computational efficiency results were obtained in this manuscript. References [1] Schuman, C.D., Kulkarni, S.R., Parsa, M., Mitchell, J.P., Date, P., Kay, B.: Opportunities for neuromorphic computing algorithms and applications. Nature Computational Science 2(1), 10–19 (2022) [2] Aimone, J.B., Date, P., Fonseca-Guerra, G.A., Hamilton, K.E., Henke, K., Kay, B., Kenyon, G.T., Kulkarni, S.R., Mniszewski, S.M., Parsa, M., et al.: A review of non-cognitive applications for neuromorphic computing. Neuromorphic Computing and Engineering 2(3), 032003 (2022) [3] Orchard, G., Frady, E.P., Rubin, D.B.D., Sanborn, S., Shrestha, S.B., Sommer, F.T., Davies, M.: Efficient neuromorphic signal processing with loihi 2. In: 2021 IEEE Workshop on Signal Processing Systems (SiPS), p. 254–259 (2021). IEEE [4] Cassidy, A.S., Arthur, J.V., Akopyan, F., Andreopoulos, A., Appuswamy, R., Datta, P., Debole, M.V., Esser, S.K., Otero, C.O., Sawada, J., et al.: 11.4 ibm northpole: An architecture for neural network inference with a 12nm chip. In: 2024 IEEE International Solid-State Circuits Conference (ISSCC), vol. 67, p. 214–215 (2024). IEEE 31 [5] Gautam, A., Date, P., Kulkarni, S., Patton, R., Potok, T.: NeuroCoreX: An open- source fpga-based spiking neural network emulator with on-chip learning. arXiv preprint arXiv:2506.14138 (2025) [6] Maheshwari, D., Young, A., Date, P., Kulkarni, S., Witherspoon, B., Miniskar, N.R.: An fpga-based neuromorphic processor with all-to-all connectivity. In: 2023 IEEE International Conference on Rebooting Computing (ICRC), p. 1–5 (2023). IEEE [7] Miniskar, N.R., Young, A.R., Asifuzzaman, K., Kulkarni, S., Date, P., Bean, A., Vetter, J.S.: Neuro-spark: A submicrosecond spiking neural networks architec- ture for in-sensor filtering. In: 2024 International Conference on Neuromorphic Systems (ICONS), p. 63–70 (2024). IEEE [8] Date, P., Potok, T., Schuman, C., Kay, B.: Neuromorphic computing is Turing- complete. In: Proceedings of the International Conference on Neuromorphic Systems 2022, p. 1–10 (2022) [9] Date, P., Kay, B., Schuman, C., Patton, R., Potok, T.: Computational complex- ity of neuromorphic algorithms. In: International Conference on Neuromorphic Systems 2021, p. 1–7 (2021) [10] Cong, G., Lim, S.-H., Kulkarni, S., Date, P., Potok, T., Snyder, S., Parsa, M., Schuman, C.: Semi-supervised graph structure learning on neuromorphic comput- ers. In: Proceedings of the International Conference on Neuromorphic Systems 2022, p. 1–4 (2022) [11] Cong, G., Kulkarni, S., Lim, S.-H., Date, P., Snyder, S., Parsa, M., Kennedy, D., Schuman, C.: Hyperparameter optimization and feature inclusion in graph neural networks for spiking implementation. In: 2023 International Conference on Machine Learning and Applications (ICMLA), p. 1541–1546 (2023). IEEE [12] Patton, R., Schuman, C., Kulkarni, S., Parsa, M., Mitchell, J.P., Haas, N.Q., Stahl, C., Paulissen, S., Date, P., Potok, T., et al.: Neuromorphic computing for autonomous racing. In: International Conference on Neuromorphic Systems 2021, p. 1–5 (2021) [13] R. Kulkarni, S., Young, A., Date, P., Rao Miniskar, N., Vetter, J., Fahim, F., Parpillon, B., Dickinson, J., Tran, N., Yoo, J., et al.: On-sensor data filtering using neuromorphic computing for high energy physics experiments. In: Proceedings of the 2023 International Conference on Neuromorphic Systems, p. 1–8 (2023) [14] Date, P., Carothers, C.D., Hendler, J.A., Magdon-Ismail, M.: Efficient classifi- cation of supercomputer failures using neuromorphic computing. In: 2018 IEEE Symposium Series on Computational Intelligence (SSCI), p. 242–249 (2018). IEEE 32 [15] Kay, B., Date, P., Schuman, C.: Neuromorphic graph algorithms: Extracting longest shortest paths and minimum spanning trees. In: Proceedings of the 2020 Annual Neuro-inspired Computational Elements Workshop, p. 1–6 (2020) [16] Kay, B., Schuman, C., O’connor, J., Date, P., Potok, T.: Neuromorphic graph algorithms: Cycle detection, odd cycle detection, and max flow. In: International Conference on Neuromorphic Systems 2021, p. 1–7 (2021) [17] Hamilton, K., Mintz, T., Date, P., Schuman, C.D.: Spike-based graph centrality measures. In: International Conference on Neuromorphic Systems 2020, p. 1–8 (2020) [18] Date, P., Kulkarni, S., Young, A., Schuman, C., Potok, T., Vetter, J.: Encoding integers and rationals on neuromorphic computers using virtual neuron. Scientific Reports 13(1), 10975 (2023) [19] Date, P., Kulkarni, S., Young, A., Schuman, C., Potok, T., Vetter, J.S.: Vir- tual neuron: A neuromorphic approach for encoding numbers. In: 2022 IEEE International Conference on Rebooting Computing (ICRC), p. 100–105 (2022). IEEE [20] Wurm, A., Seay, R., Date, P., Kulkarni, S., Young, A., Vetter, J.: Arithmetic primitives for efficient neuromorphic computing. In: 2023 IEEE International Conference on Rebooting Computing (ICRC), p. 1–5 (2023). IEEE [21] Schuman, C.D., Kay, B., Date, P., Kannan, R., Sao, P., Potok, T.E.: Sparse binary matrix-vector multiplication on neuromorphic computers. In: 2021 IEEE Inter- national Parallel and Distributed Processing Symposium Workshops (IPDPSW), p. 308–311 (2021). IEEE [22] Hamilton, K., Date, P., Kay, B., Schuman D, C.: Modeling epidemic spread with spike-based models. In: International Conference on Neuromorphic Systems 2020, p. 1–5 (2020) [23] Patton, R., Date, P., Kulkarni, S., Gunaratne, C., Lim, S.-H., Cong, G., Young, S.R., Coletti, M., Potok, T.E., Schuman, C.D.: Neuromorphic computing for scientific applications. In: 2022 IEEE/ACM Redefining Scalability for Diversely Heterogeneous Architectures Workshop (RSDHA), p. 22–28 (2022). IEEE [24] Date, P.: Combinatorial Neural Network Training Algorithm for Neuromorphic Computing. Rensselaer Polytechnic Institute, ??? (2019) [25] Date, P., Carothers, C.D., Mitchell, J.E., Hendler, J.A., Magdon-Ismail, M.: Training deep neural networks with constrained learning parameters. In: 2020 International Conference on Rebooting Computing (ICRC), p. 107–115 (2020). IEEE 33 [26] Stimberg, M., Brette, R., Goodman, D.F.: Brian 2, an intuitive and efficient neural simulator. elife 8, 47314 (2019) [27] Gewaltig, M.-O., Diesmann, M.: Nest (neural simulation tool). Scholarpedia 2(4), 1430 (2007) [28] Aimone, J.B., Severa, W., Vineyard, C.M.: Composing neural algorithms with fugu. In: Proceedings of the International Conference on Neuromorphic Systems, p. 1–8 (2019) [29] Wang, F., Kulkarni, S., Theilman, B., Rothganger, F., Schuman, C., Lim, S.-H., Aimone, J.B.: Scaling neural simulations in stacs. Neuromorphic Computing and Engineering 4(2), 024002 (2024) [30] Wang, F.: Simulation tool for asynchronous cortical streams (stacs): interfacing with spiking neural networks. Procedia Computer Science 61, 322–327 (2015) [31] Hazan, H., Saunders, D.J., Khan, H., Patel, D., Sanghavi, D.T., Siegelmann, H.T., Kozma, R.: Bindsnet: A machine learning-oriented spiking neural networks library in python. Frontiers in neuroinformatics 12, 89 (2018) [32] Eshraghian, J.K., Ward, M., Neftci, E., Wang, X., Lenz, G., Dwivedi, G., Ben- namoun, M., Jeong, D.S., Lu, W.D.: Training spiking neural networks using lessons from deep learning. Proceedings of the IEEE 111(9), 1016–1054 (2023) [33] Fang, W., Chen, Y., Ding, J., Yu, Z., Masquelier, T., Chen, D., Huang, L., Zhou, H., Li, G., Tian, Y.: Spikingjelly: An open-source machine learning infrastructure platform for spike-based intelligence. Science Advances 9(40), 1480 (2023) [34] Date, P., Gunaratne, C., R. Kulkarni, S., Patton, R., Coletti, M., Potok, T.: SuperNeuro: A fast and scalable simulator for neuromorphic computing. In: Pro- ceedings of the 2023 International Conference on Neuromorphic Systems, p. 1–4 (2023) [35] Lapicque, L.: Recherches quantitatives sur l’excitation ́electrique des nerfs trait ́e comme une polarisation. Journal de Physiologie et de Pathologie G ́en ́erale 9, 620–635 (1907) [36] Gerstner, W., Kistler, W.M.: Spiking Neuron Models: Single Neurons, Popula- tions, Plasticity. Cambridge University Press, Cambridge, UK (2002) [37] Harris, C.R., Millman, K.J., Walt, S.J., Gommers, R., Virtanen, P., Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N.J., Kern, R., Picus, M., Hoyer, S., Kerkwijk, M.H., Brett, M., Haldane, A., R ́ıo, J.F., Wiebe, M., Peterson, P., G ́erard-Marchant, P., Sheppard, K., Reddy, T., Weckesser, W., Abbasi, H., Gohlke, C., Oliphant, T.E.: Array programming with NumPy. Nature 585(7825), 357–362 (2020) https://doi.org/10.1038/s41586-020-2649-2 34 [38] Virtanen, P., Gommers, R., Oliphant, T.E., Haberland, M., Reddy, T., Courna- peau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S.J., Brett, M., Wilson, J., Millman, K.J., Mayorov, N., Nelson, A.R.J., Jones, E., Kern, R., Larson, E., Carey, C.J., Polat, ̇ I., Feng, Y., Moore, E.W., Vander- Plas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero, E.A., Harris, C.R., Archibald, A.M., Ribeiro, A.H., Pedregosa, F., van Mulbregt, P., SciPy 1.0 Contributors: SciPy 1.0: Fundamental Algorithms for Scientific Com- puting in Python. Nature Methods 17, 261–272 (2020) https://doi.org/10.1038/ s41592-019-0686-2 [39] Lawson, C.L., Hanson, R.J., Kincaid, D.R., Krogh, F.T.: Basic linear algebra subprograms for fortran usage. ACM Transactions on Mathematical Software (TOMS) 5(3), 308–323 (1979) [40] Sen, P., Namata, G., Bilgic, M., Getoor, L., Gallagher, B., Eliassi-Rad, T.: Collective classification in network data. In: AI Magazine, vol. 29, p. 93–106 (2008) [41] Namata, G., London, B., Getoor, L., Huang, B., EDU, U.: Query-driven active surveying for collective classification. In: 10th International Workshop on Mining and Learning with Graphs (MLG) (2012) [42] Kipf, T.N., Welling, M.: Semi-supervised classification with graph convolutional networks. In: International Conference on Learning Representations (ICLR) (2017) [43] Yang, Z., Cohen, W.W., Salakhutdinov, R.: Revisiting semi-supervised learning with graph embeddings. In: Proceedings of the 33rd International Conference on Machine Learning (ICML), p. 40–48 (2016) [44] Sironi, A., Brambilla, M., Bourdis, N., Lagorce, X., Benosman, R.: Hats: His- tograms of averaged time surfaces for robust event-based object classification. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recog- nition, p. 1731–1740. IEEE, ??? (2018). https://doi.org/10.1109/CVPR.2018. 00186 [45] Diehl, P.U., Cook, M.: Unsupervised learning of digit recognition using spike- timing-dependent plasticity. Frontiers in Computational Neuroscience 9, 99 (2015) https://doi.org/10.3389/fncom.2015.00099 [46] Lagorce, X., Orchard, G., Galluppi, F., Shi, B.E., Benosman, R.B.: Hots: A hierarchy of event-based time-surfaces for pattern recognition. IEEE Transac- tions on Pattern Analysis and Machine Intelligence 39(7), 1346–1359 (2017) https://doi.org/10.1109/TPAMI.2016.2574707 [47] Bi, Y., Chadha, A., Abbas, A., Bourtsoulatze, E., Andreopoulos, Y.: Graph- based object classification for neuromorphic vision sensing. In: Proceedings of the 35 IEEE/CVF International Conference on Computer Vision, p. 491–501 (2019). https://doi.org/10.1109/ICCV.2019.00058 [48] Sedgewick, R., Wayne, K.: DijkstraSP.java: Dijkstra’s Shortest Path Algo- rithm. https://algs4.cs.princeton.edu/code/edu/princeton/cs/algs4/DijkstraSP. java.html. Part of the Algorithms, 4th Edition codebase. Accessed: 2026-08-02 36