1 Introduction↩︎

The human brain, recognized as one of the most sophisticated computational systems on the planet, demonstrates unparalleled energy efficiency and cognitive capabilities. Spiking Neural Networks (SNNs) have been proposed as one of the most representative brain-inspired computational models, aiming to mimic the efficient spatiotemporal information processing in the brain [1]. In contrast to traditional artificial neural networks (ANNs) that rely on real-valued neural representations and continuous activation functions, SNNs utilize discrete spike trains to represent information and inherently support efficient event-driven computation. Furthermore, spiking neurons can incorporate rich neuronal dynamics for effective temporal processing [2]. At present, SNNs have demonstrated competitive performance and substantial energy savings compared to traditional ANNs [3]–[5] in a wide range of applications, such as image classification [6], [7], audio processing [8], [9], and robotic control [10], [11].

While SNNs have gained increasing attention in recent years, their main applications have been primarily confined to computer vision tasks that involve limited temporal dynamics, such as classifying static images or data collected from dynamic vision sensors (DVS) with artificially added saccade motion [12], [13]. However, real-world scenarios are often more challenging as they involve sensory signals with information that varies across different timescales. For instance, speech recognition tasks often necessitate models to establish dependencies across various timescales, encompassing phonemes, words, and sentences. While several algorithms have been developed to enhance the temporal processing capacity of SNNs [14]–[16], most of them still struggle to establish diverse scales of temporal representations, resulting in poor performance in complex temporal processing tasks.

Considering the remarkable temporal processing capabilities observed in biological neurons, it is crucial to carefully examine their computational mechanisms. In the context of large-scale SNNs, the majority of existing spiking neurons are modeled as single-compartment systems, such as the Leaky Integrate-and-Fire (LIF) model [17]. In these single-compartment models, the neuronal dynamics are simplified to first-order dynamics of a single state variable - the membrane potential, disregarding any interactions among substructures within a single neuron. This simplification reduces computational complexity but limits their ability to generate complex neuronal dynamics.

In contrast, real biological neurons can be better modeled by multi-compartment neuron models [18]–[20], which divide a single neuron into several interconnected subunits with interactions among them. Extensive neuroscience experiments have demonstrated the important roles of these interactive dynamics in temporal processing. For example, the nonlinear interactions among ion channels endow neurons with significant computational power for processing temporal signals [21], [22]. Additionally, the interaction among coupled dendritic components acts as temporal filters for signals traveling from dendrites to the soma, thereby facilitating the temporal sequence detection [19], [23]. While detailed multi-compartment models show great potential in temporal processing, they are computationally expensive due to the intricate anatomical structures of dendritic trees and high-dimensional ionic properties. Consequently, they are not well suited for constructing large-scale SNNs to tackle real-world pattern recognition tasks. Therefore, there is a pressing need to develop an efficient spiking neuron model that strikes a balance between computational complexity and the ability to capture the valuable interaction among different substructures of a biological neuron.

Recently, several simplified compartmental spiking neuron models have been specifically designed for deep SNNs. For instance, a two-compartment model was introduced to simulate double-exponential threshold decay [24], and a dendritic neuron model was proposed to endow multiple dendritic compartments with heterogeneous decaying time constants [25]. Furthermore, inspired by the well-known Pinsky-Rinzel (P-R) pyramidal neuron located in the CA3 region of the hippocampus, a two-compartment LIF model was proposed, which divides a single neuron into dendritic and somatic compartments [26].

Nevertheless, there are three significant challenges that remain to be tackled. Firstly, these models overlook the valuable recurrent interactions among neuronal compartments, which play a crucial role in integrating temporal information across multiple timescales in biological neurons [27], [28]. As a result, the expressibility of these models is limited, constraining the complexity of temporal dynamics they can effectively capture. Secondly, these meticulously handcrafted models are inherently constrained by a predetermined number of compartments. While increasing this number allows for richer structural complexity and a broader range of temporal dynamics, reducing it yields a more compact and computationally efficient representation. However, the inability to adjust this abstraction level limits these models’ adaptability to tasks with varying temporal complexity and computational demands. Thirdly, compared to single-compartment models, the increased temporal dynamics in these models result in significantly slower training processes when using the back-propagation through time (BPTT) algorithm. Consequently, these models demonstrate limited scalability to larger networks for tackling challenging real-world temporal processing tasks, particularly those involving long sequences.

To tackle these challenges, we propose a generalized multi-compartment spiking neuron model for SNNs that integrates essential recurrent interactions among interconnected compartments. Moreover, the number of compartments can be flexibly adjusted in the proposed model, allowing for adaptation to the different temporal complexities required in real-world tasks. Furthermore, considering the significant constraint on training speed caused by increasingly complex neuronal dynamics in the proposed model, we introduce two temporal parallelization techniques to accelerate the training process. Firstly, for the linear recurrence in the neuron charging process, we model it as a linear time-invariant (LTI) system, which can be efficiently parallelized over time following existing approaches for linear recurrent models [29]–[31]. Furthermore, for the nonlinear recurrence associated with the neuron firing and reset process, we introduce a novel reset strategy. This mechanism effectively decouples the nonlinear temporal dependency, enabling parallel computation that is mathematically equivalent to a step-by-step, spike-triggered reset mechanism. Altogether, these strategies enable efficient temporal parallelization during training, while preserving the event-driven, online inference on neuromorphic hardware.

The proposed design methodologies have led to the development of a novel Parallel Multi-compartment Spiking Neuron (PMSN) model as illustrated in Fig. 1. The PMSN model effectively captures multiscale temporal information through the interaction among neuronal compartments. Both theoretical analysis and dynamics visualizations are provided to substantiate the efficacy of the PMSN model in establishing temporal dependencies across various timescales. Moreover, the proposed PMSN model offers a significant improvement in training speed, particularly when deployed on GPU-accelerated machine learning (ML) frameworks. The main contributions of this work are summarized as follows:

• We propose a generalized multi-compartment spiking neuron model for deep SNNs that incorporates valuable interactions among different neuronal compartments. The number of compartments in this model can be flexibly adjusted to represent temporal information across diverse timescales.

• We develop two temporal parallelization techniques for the proposed model, achieving more than a \(10\times\) training speedup on GPU-accelerated ML frameworks.

• Our comprehensive experiments demonstrate the superior temporal processing capacity and substantial training acceleration of PMSN, as well as its potential for energy-efficient deployment on neuromorphic platforms.

2 Related Works↩︎

2.1 Advances in Spiking Neuron Models↩︎

Recent enhancements in spiking neuron models for temporal processing can be categorized into single- or multi-compartment models according to the number of membrane potential subunits involved [32]. Several single-compartment models have incorporated adaptive variables to enhance the efficacy of temporal information representation. For instance, [14] and [15] propose to use adaptive firing thresholds, which function as long-term memory to enhance temporal processing; [33] utilizes learnable decay constants to enhance the heterogeneity of neuron populations, enabling them to effectively represent multiscale temporal information; [34] introduces a gating mechanism into the neuron model to explicitly control memory storage and retrieval processes.

Comparatively, multi-compartment models utilize compartmental heterogeneity or recurrent interaction to represent information across different timescales. For example, [24] proposes a two-compartment model where the thresholds undergo a double-exponential decay, facilitating the storage of both short- and long-term information; [25] introduces a heterogeneous model consisting of one somatic compartment and varied dendrite-branch compartments; [26] proposes a two-compartment model that captures the interactive dynamics between the soma and dendrites. Although these handcrafted models demonstrate enhanced performance, there is a pressing need for a principled approach to efficiently scale them and incorporate a variable number of compartments. This is crucial for enhancing their temporal processing capability, which serves as the central focus of this paper.

2.2 Training Techniques for SNNs↩︎

Training SNNs poses fundamental challenges for gradient-based optimization due to the non-differentiable nature of spike generation [35]. Existing training methods for deep SNNs can be broadly categorized into conversion-based and direct training approaches. Conversion-based methods transfer knowledge from pre-trained ANNs to SNNs by matching activation or firing-rate statistics [4], [36]. While effective for static tasks, such approaches are generally unsuitable for temporal processing, as they introduce non-negligible approximation errors in temporal dynamics. In contrast, direct training methods optimize SNNs end-to-end by approximating spike gradients with surrogate functions, enabling effective learning on temporal signals [37], [38]. Recent advances have further demonstrated their potential in efficient language modeling [39], robust visual recognition [40], and event-based optical flow estimation [41]. However, due to the recurrent nature of spiking neuron dynamics, most direct training approaches rely on BPTT, or its online variants [42], resulting in training complexity that scales linearly with sequence length. This severely limits their scalability to long sequences.

Several recent works have explored parallelization strategies to alleviate this limitation. Differentiation on Spike Representation (DSR) reformulates spiking dynamics using differentiable rate-based representations, enabling constant-time gradient computation [6]. However, this approach is restricted to simple Integrate-and-Fire (IF) or LIF neurons and relies on rate-based representations, making it unsuitable for modeling fine-grained temporal dependencies. In parallel, the Parallel Spiking Neuron (PSN) family [43] reformulates membrane potential charging dynamics using a learnable decay matrix and omits the reset mechanism to enable temporal parallelism. Despite their efficiency, PSN models require access to future inputs and are limited to single-compartment structures, making them biologically implausible, incompatible with neuromorphic hardware, and difficult to extend to richer multi-compartment temporal modeling.

3 Revisit: Single-compartment Spiking Neuron Models Struggle to Process Temporal Signals↩︎

In this section, we revisit the LIF model as a representative of single-compartment spiking neuron models. We first introduce basic concepts and then discuss the major challenges associated with processing temporal signals using this model. As illustrated in Fig. 1, the temporal dynamics of the LIF model can be formulated as: \[\begin{align} &\text{Leak \& Charge:} \quad \frac{\mathrm{d} v(t)}{\mathrm{d}t}=-\frac{1}{\tau_m }(v(t)-v_{rest})+I(t),\\ &\text{Fire \& Reset:} \quad \text{if} \;v(t)\geq \theta, s(t)=1, {\;v(t)\rightarrow v(t)-\theta}. \end{align}\]

During the membrane potential leaky and charging phase, the information contained in the input current \(I(t)\) is integrated into the membrane potential \(v\) and further undergoes decay at a rate governed by \(\tau_m\). \(v_{rest}\) is the resting potential. Once \(v(t)\) exceeds the firing threshold \(\theta\), an output spike will be generated and transmitted to subsequent neurons, after which the membrane potential is reset. In practice, the above continuous-time formulation is typically discretized using the Euler method as \[\label{eq:lif} \left\{\begin{array}{lr} V[t]=\alpha V[t-1] +I[t] -\theta S[t-1],& \\ S[t]=H\left(V[t] {-}\theta\right), & \end{array}\right.\tag{1}\] where \(H(\cdot)\) is the Heaviside function and \(\alpha=\mathrm{exp}(-\frac{\mathrm{d}t}{\tau_m})\) is the discrete-time decay factor. Following this discretization, LIF neuron dynamics are simulated via recurrent state updates over time, where the membrane potential \(V[t]\) depends on its previous state \(V[t{-}1]\) as well as the current input. More generally, such recurrent evolution is a fundamental characteristic of SNNs, allowing them to be naturally interpreted as recurrent models with intrinsic temporal dependencies.

Despite its promising results in tasks involving limited temporal context, LIF neurons encounter the following two challenges when dealing with long sequences. Firstly, this model faces challenges in retaining information over an extended time period, primarily due to its inadequate representation of neuronal dynamics. Specifically, the LIF model overlooks the computational role played by interactions among different neuronal substructures. As a result, the membrane potential \(V\) is the only state variable that can integrate and store temporal information. Unfortunately, this state variable is subject to exponential decay and reset, which hinders the establishment of long-term temporal dependencies. Secondly, its training time increases proportionally with the sequence length, making it inefficient for tasks that involve long sequences. This slow training speed arises from the nonlinear recurrent updates in the neuronal dynamics, which must be computed step by step along the temporal dimension as commonly observed in RNNs, leading to underutilization of the full potential of GPU acceleration. Specifically, the membrane potential update of \(V[t]\) depends on the output spike \(S[t-1]\) from the preceding timestep, which is generated by applying a nonlinear Heaviside activation to \(V[t-1]\). As a result, \(V[t]\) is not available until \(V[t-1]\) has been fully processed. This time-coupled relationship prevents the membrane potential dynamics from being unfolded across time, leading to difficulties in temporal parallelization.

4 A Generalized Multi-compartment Spiking Neuron Model with Interactive Dynamics↩︎

In this section, we first briefly introduce the multi-compartment models that have been widely used in neuroscience studies. Drawing inspiration from them, we propose a generalized multi-compartment spiking neuron model that can strike a favorable balance between the richness of neuronal dynamics and computational cost. Notably, in contrast to commonly used multi-compartment neuron models in deep SNNs, our model retains the crucial recurrent interactions among interconnected compartments. It thereby allows more effective representation and processing of temporal signals across different timescales.

Compartmental models divide a single neuron into interconnected subunits or “compartments”, with specific spatial structures. The interactions among these interconnected compartments lead to the rich neuronal dynamics observed in biological neurons that play a crucial role in sensory processing [44]. For instance, Rall’s cable theory [23] suggests each dendrite can be mathematically modeled as a series of interconnected compartments resembling a cable structure. However, given the complex anatomical morphology of dendritic trees, these models pose significant computational challenges and are impractical for large-scale simulations. To reduce the modeling difficulty, quantitative models with reduced compartment numbers and structure complexity have been explored. For instance, the neuronal activities of CA3 pyramidal neurons can be modeled by a 19-compartment cable model, wherein dendritic branching is omitted [45]. This model has been further abstracted to a more computationally efficient two-compartment model [46]. In these quantitative models, compartments are restricted to interacting with their adjacent compartments solely, and the dynamics of each compartment can be described as a differential equation of a resistor-capacitance circuit as \[\label{eq:32real95neuron} \begin{align} C_m \frac{\mathrm{d} v^{(i)}}{\mathrm{d} t}=& -I_{Leak}( v^{(i)}) + g_{i-1,i}(v^{(i-1)}-v^{(i)})\\ &+ g_{i,i+1}(v^{(i+1)}-v^{(i)}) + I. \end{align}\tag{2}\] where \(C_m\) is the capacitance of membrane potential, \(I_{Leak}(\cdot)\) represents the leakage current of membrane potential \(v^{(i)}\), \(g_{i,j}(\cdot)\) is the coupling conductance between compartments \(i\) and \(j\), and \(I\) denotes the summation of other voltage-gated ionic currents.

In order to simplify the ionic compartmental dynamics described in Eq. 2 and ensure compatibility with existing GPU-accelerated ML frameworks, we next introduce a generalized multi-compartment spiking neuron model with a simplified neuronal hyperparameter setup. This model retains the essential interactions among different substructures of a neuron, which are theoretically evidenced in Section 6 as playing a crucial role in capturing and integrating temporal features across multiple timescales. Another notable feature of our model is its flexibility, allowing for a varying number of neuronal compartments \(n\) to suit the complexity of the task at hand. The detailed model formulation is given as \[\begin{align} \label{eq:generalized} &\left\{\begin{array}{lr} \frac{\mathrm{d} v^{(1)}(t)}{\mathrm{d} t}=-\frac{1}{\tau_1}v^{(1)}(t) +\beta_{2,1}v^{(2)}(t) + \gamma_1 I(t), & \\ \frac{\mathrm{d} v^{(2)}(t)}{\mathrm{d} t}=-\frac{1}{\tau_2}v^{(2)}(t)+\beta_{3,2}v^{(3)}(t)+\beta_{1,2}v^{(1)}(t)+\gamma_2 I(t), & \\ \qquad \qquad \qquad ... &\\ \frac{\mathrm{d} v^{(n)}(t)}{\mathrm{d} t}=-\frac{1}{\tau_n}v^{(n)}(t)+\beta_{n-1,n}v^{(n-1)}(t)+\gamma_n I(t), & \end{array}\right.\\ &\begin{array}{ll} \text{if} \;&v^{(n)}(t)\geq \theta, \quad \;s(t)=1, \quad v^{(n)}(t)\leftarrow v^{(n)}(t)-\theta, \end{array} \end{align}\tag{3}\] where \(v^{(i)}\) represents the membrane potential of the compartment \(i\), \(\theta\) is the firing threshold. \(I^l(t)=\mathcal{W}^l S^{l-1}(t)\) denotes the synaptic current induced by the output spikes of the preceding layer, where \(\mathcal{W}^l\) represents the synaptic weight between layer \(l-1\) and \(l\). Once the membrane potential of the final compartment \(v^{(n)}\) exceeds the threshold \(\theta\), it triggers an output spike and simultaneously resets \(v^{(n)}\). The compartmental parameters, including \(\tau_i\), \(\gamma_i\), and \(\beta_{i,j}\), represent the membrane time constant, input attenuation of the compartment \(i\), and coupling strength between interconnected compartments \(i\) and \(j\), respectively.

For the sake of simplicity in representation, these first-order differential equations can be equivalently expressed as an \(n\)-dimensional first-order state-space dynamical system: \[\label{eq:multicompartment} \begin{align} \dot{\mathcal{V}}(t)&= \begin{bmatrix} -\frac{1}{\tau_1} & \!\beta_{2,1}\! & \!0\! & \!\mathclap{\cdots}\! & \!0\! \\ \!\beta_{1,2}\! & \!-\frac{1}{\tau_2}\! & \!\beta_{3,2}\! & \!\mathclap{\cdots}\! & \!0\! \\ \!\mathclap{\vdots}\! & \!\mathclap{\vdots}\! & \!\mathclap{\vdots}\! & \!\ddots\! & \!\mathclap{\vdots}\! \\ \!0\! & \!0\! & \!\mathclap{\cdots}\! & \!-\frac{1}{\tau_{n-1}}\! & \!\beta_{n,n-1}\! \\ \!0\! & \!0\! & \!\mathclap{\cdots}\! & \!\beta_{n-1,n}\! & \!-\frac{1}{\tau_n}\! \\ \end{bmatrix}\mathcal{V}(t) + \boldsymbol{\gamma_n} I(t),\\ S(t)&= H(v^{(n)}(t) - \theta ), \qquad v^{(n)}(t) \leftarrow v^{(n)}(t) - \theta S(t), \end{align}\tag{4}\] where \(\boldsymbol{\gamma_n}=[\gamma_1, \cdots, \;\gamma_n]^T\), \(\mathcal{V}=[v^{(1)} , \cdots, \;v^{(n)}]^T\).

It is important to note that many widely used spiking neuron models can be seen as special cases of our proposed model by adjusting its hyperparameters. For instance, when the number of compartments \(n\) is set to \(1\) and \(\beta=0\), our model is simplified to a standard LIF model. Similarly, when \(n=2\) and \(\gamma_2=0\), our model reduces to a Two Compartment LIF (TC-LIF) model [26]. These examples highlight the generalizability of our model. Moreover, by providing flexibility in increasing the number of compartments, our model surpasses existing spiking neuron models and offers richer neuronal dynamics that are essential for complex temporal processing tasks.

5 PMSN: A Parallel Multi-compartment Spiking Neuron↩︎

The generalized multi-compartment spiking neuron model introduced earlier requires significantly more training time compared to existing spiking neuron models used in deep SNNs, primarily due to its higher computational complexity. This limitation restricts its practical deployment in long sequence temporal processing scenarios. Thus, it becomes crucial to develop strategies to enable parallel training of this model in time. However, it is not straightforward to employ existing parallel computation techniques [29], due to the involvement of the nonlinear function \(H(\cdot)\) for spike generation and reset processes.

In this section, we further propose a PMSN model with two parallelization techniques independently designed for hidden and output compartments. To better explain our idea, we divide the total \(n\) compartments into \(n-1\) hidden compartments \(\mathcal{V}_h\) with linear recurrence and one output compartment \(v^{(n)}=v_s\) with nonlinear firing and reset dynamics. Meanwhile, the feedback from the output to the last hidden compartment \(\beta_{n,n-1}\) is set to \(0\). The resulting neuronal function can be represented as

\[\label{eq:enhance95h} \dot{\mathcal{V}}_h(t)= \begin{bmatrix} -\frac{1}{\tau_1} & \!\beta_{2,1}\! & \!0\! & \!\mathclap{\cdots}\! & \!0\! \\ \!\beta_{1,2}\! & \!-\frac{1}{\tau_2}\! & \!\beta_{3,2}\! & \!\mathclap{\cdots}\! & \!0\! \\ \!\mathclap{\vdots}\! & \!\mathclap{\vdots}\! & \!\mathclap{\vdots}\! & \!\ddots\! & \!\mathclap{\vdots}\! \\ \!0\! & \!0\! & \!0\! & \!\mathclap{\cdots}\! & \!-\frac{1}{\tau_{n-1}}\! \\ \end{bmatrix} \mathcal{V}_h(t) + \boldsymbol{\gamma_{n-1}} I(t),\tag{5}\] \[\label{eq:enhance95s} \dot{v}_s(t)=\beta_{n-1,n} v^{(n-1)}(t) -\frac{1}{\tau_{n}}v_s(t)+\gamma_n I(t)-\theta S(t),\tag{6}\] \[\label{eq:enhance95re} S(t)= H(v_s(t) - \theta ).\tag{7}\]

Below, we will introduce how these two parallelization techniques unfold the linear recurrent states \(\mathcal{V}_h\), and decouple the nonlinear temporal dependency issue associated with \(v_s\) to enable parallel computation.

5.1 Parallel Strategy for Hidden Compartments↩︎

By examining the hidden compartment dynamics of PMSN in Eq. 5 , we observe that the hidden states \(\mathcal{V}_h\) evolve according to an LTI dynamical system. Specifically, the membrane potential of each hidden compartment \(v^{(i)}\) is updated by applying a time-invariant state transition matrix to its previous state. As a result, the hidden compartment dynamics admit a state-space formulation, which shares a common underlying mathematical structure with linear recurrent neural networks (RNNs) [29], [47] and state-space models (SSMs) [30], [31] widely used in modern sequential modeling. This connection bridges biologically grounded neuron dynamics with contemporary deep learning architectures and, importantly, allows PMSN to leverage parallel computation techniques originally developed for parameterized SSMs to efficiently compute the temporal evolution of its hidden compartments.

We first apply the zero-order hold (ZOH) method [48] to discretize the nonlinear continuous system detailed in Eqs. 5 and 6 . Specifically, we utilize a full-rank state transition matrix \(\mathcal{T}\in \mathbb{R}^{ (n-1)\times (n-1) }\) to represent the first matrix in Eq. 5 , and diagonalize this matrix via eigenvalue decomposition \(\mathcal{T}=P\Lambda P^{-1}\), where \(\Lambda, P \in \mathbb{C}^{(n-1)\times (n-1)}\) denote the diagonal eigenvalue matrix and eigenvector matrix, respectively. We then obtain the following discrete-time formulation: \[\label{eq:model95woreset1} {V}_h[t]= \Bar{\mathcal{T}} {V}_h[t-1] + {{\Phi}_c} I[t], \qquad \qquad\tag{8}\] \[\label{eq:I95h} I_h[t]={\Phi}_s {V}_h[t] + \gamma_n I[t], \qquad \qquad \quad \;\tag{9}\] \[\label{eq:model95woreset2} \begin{align} &{v_s}[t]=\alpha v_s[t-1] + I_h[t] -\theta S[t-1], \\ &S[t]= H(v_s[t] - \theta ), \end{align}\tag{10}\] where \(V_h=P^{-1}\mathcal{V}_h\), \(\Bar{\mathcal{T}}=\mathrm{exp}(\Lambda \mathrm{d}t)\), \({\Phi}_c=\Lambda^{-1}(\mathrm{exp}(\Lambda \mathrm{d}t)-I)\phi_c\), and \(\phi_c=P^{-1}\boldsymbol{\gamma_{n-1}}.\) The term \(I_h[t]\) signifies the total input to the output compartment, \({\Phi}_s=[0,..,\beta_{n-1,n}]P\), \(\alpha = \mathrm{exp}(-\frac{\mathrm{d}t}{\tau _n})\). The parameters \(\Lambda \mathrm{d}t,~\phi_c,~\gamma_n\), and \(\Phi_s\) are learnable. It should be noted that these complex number operations are well-supported by existing neuromorphic chips, such as Intel Loihi [49].

The model described in Eq. 8 exhibits a linear recurrence that can be unfolded over time as \[{V}_h[t]=\sum^t_{i=0}{\Bar{\mathcal{T}}^{t-i}{{\Phi}_c}I[i] }. \label{eq:v95h95t}\tag{11}\]

By substituting Eq. 11 into Eq. 9 , we can obtain \[\label{eq:I95h952} I_h[t]=\sum^t_{i=0}\Phi_s\Bar{\mathcal{T}}^{t-i}{{\Phi}_c}I[i] +\gamma_n I[t].\tag{12}\]

Notably, the first term in Eq. 12 can be simplified to a convolution form as \[\label{eq:I95h953} I_h[t] = \sum_{i=0}^t\boldsymbol{I_t}[i]\mathcal{K}[t-i] +\gamma_n I[t] =(\boldsymbol{I_t} * \mathcal{K})[t] +\gamma_n I[t],\tag{13}\] where \(\boldsymbol{I_t}=\{I[0],..,I[t]\}\) denotes the input current set, \(\mathcal{K}=[\Phi_s\Bar{\mathcal{T}}^{0}\Phi_c,\;...,\;\Phi_s\Bar{\mathcal{T}}^{t}\Phi_c]\) is the convolution kernel. Consequently, the computation of set \(\boldsymbol{I_h}=\{I_h[0],...,I_h[t]\}\) can be parallelized over time by applying a convolution operation on \(\boldsymbol{I_t}\) and \(\mathcal{K}\), resulting in \[\label{eq:fft} \boldsymbol{I_h} = \boldsymbol{I_t} * \mathcal{K} + \gamma_n \boldsymbol{I_t}= \mathcal{F}^{-1}(\mathcal{F}(\boldsymbol{I_t})\cdot \mathcal{F}(\mathcal{K}))+\gamma_n \boldsymbol{I_t},\tag{14}\] where \(\mathcal{F},~\mathcal{F}^{-1}\) represent forward and inverse Fourier transforms, respectively. In this way, we could efficiently compute the membrane potential for the first \(n-1\) hidden compartments \(V_h\), and the input for the output compartment \(I_h\) across all time steps in parallel.

5.2 Parallel Strategy for Output Compartment with Reset↩︎

In the previous subsection, we demonstrated that the dynamics of PMSN within the hidden compartments can be efficiently computed in parallel. The remaining question is how to achieve parallel computation in the output compartment.

However, achieving temporal parallelization for the output compartment is non-trivial due to the spike generation and reset mechanism, which introduces a nonlinear dependency across consecutive timesteps. Specifically, in conventional reset formulations in Eq. 10 , the membrane potential \(v_s[t]\) is updated based on the previous spike \(S[t-1]\), where \(S[t-1]\) is generated by applying a non-differentiable Heaviside function to \(v_s[t-1]\). This spike-triggered reset tightly couples \(v_s\) across time, preventing the membrane dynamics from being unfolded into a closed-form expression and fundamentally limiting temporal parallelization.

To address this limitation, we propose a novel reset mechanism termed temporal-decoupled reset (TDR). TDR reformulates the reset operation such that it can be evaluated either in a serial (step-by-step) manner or in a temporally parallel manner, while preserving the same spike-driven reset behavior. Concretely, we introduce a reset term \(v_r[t]\) that represents the discharged membrane voltage after spiking. In the serial form, \(v_r[t]\) is expressed as an integer multiple of the firing threshold \(\theta\), determined by the membrane potential \(v_s[t]\). Its dynamics can be formulated as

\[\label{eq:modified95out} \begin{align} &v_s[t]= {v}_s[t-1]+I_h[t] - v_r[t-1], \\ &S[t]=H(v_s[t] - \theta),\\ &v_r[t]= S[t] \cdot \lfloor{ {v}_s[t]}\rfloor_{\theta}, \end{align}\tag{15}\] where \(\lfloor x \rfloor_{\theta}:=\theta \lfloor x/\theta \rfloor\) signifies the floor division by \(\theta\).

Similar to the commonly used reset mechanisms in conventional spiking neuron models, TDR can be viewed as a biologically inspired reset mechanism that abstracts the depolarization–fire–repolarization process into a simplified form suitable for efficient computation and neuromorphic hardware implementation. Specifically, soft reset removes one fixed threshold after spike generation, while hard reset clamps the membrane potential to a reset value. Following the same abstraction, TDR differs only in the amount of membrane voltage discharged after spike generation: it removes a threshold-quantized voltage determined by the floor operation, which helps keep the post-reset membrane potential within the subthreshold range. Importantly, the same reset dynamics can also be computed in a temporally parallel form. Specifically, the accumulated discharged voltage \(\sum_{i=0}^{t-1} v_r[i]\) can be directly derived from the cumulative input current \(\sum_{i=0}^{t-1} I_h[i]\) across time, thereby bypassing the evolution of the membrane potential \(v_s[t]\) as

\[\begin{align} \label{eq:sumreset} &\sum^{t-1}_{i=0} v_r[i]= \lfloor{\sum^{t-1}_{i=0} I_h[i]}\rfloor_{\theta}, \\ &v_s[t]= \sum^{t}_{i=0} I_h[i] - \lfloor{\sum^{t-1}_{i=0} I_h[i]}\rfloor_{\theta}, \end{align}\tag{16}\]

A detailed derivation is provided in the Supplementary Materials. To ensure equivalence between Eq. 16 and Eq. 15 , the input current \(I_h\) is passed through a rectified linear unit (ReLU) to enforce non-negativity before being injected into the output compartment. As a result, the parallel formulation of TDR can be expressed by directly deriving the output spike train \(\boldsymbol{S_t}=\{S[0],..,S[t]\}\) from the corresponding membrane potential sequence \(\boldsymbol{v_{s,t}}=\{v_s[0], ..., v_s[t]\}\) as \[\boldsymbol{S_t} = H\left( \boldsymbol{v_{s,t}}\right)=H\left(\boldsymbol{I_{h,t}}- \lfloor{\boldsymbol{I_{h,t-1}}}\rfloor_{\theta} \right). \label{eq:32output95parallel}\tag{17}\] where the set \(\boldsymbol{I_{h,t}}=\{ I_h[0],...,\sum^{t}_{i=0} I_h[i]\}\), representing the cumulative input current to the output compartment over time, can be efficiently computed using the parallel prefix sum (Scan) algorithm [50].

Compared with existing reset mechanisms summarized in Table 1, the proposed TDR exhibits several distinctive advantages. Unlike soft reset, hard reset, or fused reset strategies, which inherently rely on spike-dependent, step-by-step updates, TDR enables temporally parallel computation during training, substantially improving training efficiency. In contrast to reset-free or Surrogate Dynamic Network (SDN) approaches [54], TDR preserves spike sparsity and maintains equivalent online, event-driven inference behavior, ensuring full compatibility with neuromorphic hardware deployment. Moreover, the proposed reset mechanism can be generalized to other single-compartment neuron models or multi-compartment models with multiple spike-generating output compartments, thereby enabling parallel training for a broader class of spiking neuron architectures.

6 Effective Temporal Gradient Propagation↩︎

In this section, we provide a theoretical analysis of how the proposed PMSN supports effective multiscale temporal processing. We first derive the gradient flow of PMSN and then discuss how the initialization of compartmental parameters stabilizes the neuronal dynamics and enables multiscale temporal dependency modeling.

To overcome the gradient discontinuity introduced by the Heaviside and floor functions in the proposed reset operation, we employ surrogate gradient methods [66], [67] to facilitate gradient backpropagation in Eq. 17 , yielding a well-defined gradient flow for the PMSN model: \[\begin{small} \begin{align} \label{eq:gradw} \Delta \mathcal{W}^{l}& \propto \frac{\partial \mathcal{L}}{\partial \mathcal{W}^{l}}= \sum_{t=1}^{T}{\frac{\partial \mathcal{L} }{\partial I^{l}[t]}S^{l-1}[t]}, \;\Delta b^{l} \propto \frac{\partial \mathcal{L} }{\partial b^{l}}= \sum_{t=1}^{T}{\frac{\partial \mathcal{L} }{\partial I^{l}[t]}}, \\ \frac{\partial \mathcal{L}}{\partial I^{l}[t]} &=\sum^T_{i=t}\frac{\partial \mathcal{L}}{\partial S^{l}[i]}\frac{\partial S^{l}[i]}{\partial v_s^l[i]}\frac{\partial v_s^l[i]}{\partial I^{l}[t]}+ \sum^{T-1}_{i=t}\frac{\partial \mathcal{L}}{\partial v_s^{l}[i+1]}\frac{\partial v_s^{l}[i+1]}{\partial v_s^{l}[i]}\frac{\partial v_s^l[i]}{\partial I^{l}[t]}\\ &=\underbrace{\frac{\partial \mathcal{L}}{\partial S^{l}[t]}g'[t]\gamma_n}_{\text{Spatial}}+ \underbrace{\sum^T_{i=t}\frac{\partial \mathcal{L}}{\partial S^{l}[i]}g'[i]\Phi_s\Bar{\mathcal{T}}^{i-t}\Phi_c}_{\text{Temporal}}, \end{align}\end{small}\tag{18}\] where \(g'[t]=\frac{\partial S^l[t]}{\partial v_s^l[t]}\) is the surrogate gradient function, \(\mathcal{L}\) is the loss function, and \(\mathcal{W}^l\), \(b^l\) refer to weight and bias terms of layer \(l\), respectively. The first term represents the spatial credit assignment between consecutive network layers, where the gradient propagates from deeper to shallower layers to facilitate hierarchical learning. The second term corresponds to the temporal credit assignment, where gradients flow backward from later to earlier timesteps, facilitating the learning of temporal dependencies across time.

The temporal credit assignment in PMSN is governed by the diagonal state transition matrix \(\Bar{\mathcal{T}}=\mathrm{diag}(\lambda_1,\ldots,\lambda_{n-1})\) after eigenvalue decomposition. Each complex-valued eigenvalue \(\lambda_i\) characterizes the intrinsic temporal dynamics of one compartment. Specifically, the magnitude of \(\lambda_i\) determines the decay rate of temporal information and thus controls the effective memory length, while its phase introduces oscillatory dynamics that arise from inter-compartmental interactions and support temporal integration across different timescales. Together, these eigenvalues determine how PMSN preserves and integrates temporal information across different timescales. Therefore, the initialization of compartmental parameters is important, as it determines the initial distribution of these eigenvalues.

To provide a stable and expressive starting point, we adopt a structured initialization strategy for the compartmental parameters. The membrane leakage terms on the main diagonal of Eq. 5 keep the compartmental dynamics in a stable regime, preventing uncontrolled divergence and supporting long-sequence information retention, as theoretically analyzed in the Supplementary Materials. Meanwhile, the coupling terms on the sub- and super-diagonal entries introduce heterogeneous interactions among neighboring compartments, thereby distributing the dynamics across multiple timescales and enabling PMSN to capture multiscale temporal dependencies, as demonstrated in Fig. 4.

7 Experimental Results↩︎

In this section, we evaluate the proposed PMSN model with a focus on its multiscale temporal processing capacity, static image classification accuracy, simulation acceleration, and computation efficiency. Unless otherwise stated, all tested PMSNs have \(n=5\) compartments to ensure a comparable computational cost to the state-of-the-art (SOTA) parallel spiking neuron model (i.e., 32-receptive-field sliding parallel spiking neuron (SPSN) [43]). We also provide ablation studies in Section 7.7 as well as the detailed experimental settings in Supplementary Materials.

7.1 Temporal Dependency Modeling across Distinct Timescales↩︎

We first compare our PMSN model against other SOTA models on temporal processing tasks involving long-term temporal dependencies. Our experiments begin with three widely-used neuromorphic benchmarks that are commonly adopted to evaluate temporal dependency modeling, including Sequential MNIST (S-MNIST) and Permuted Sequential MNIST (PS-MNIST) datasets with \(784\) time steps [68], and Spiking Heidelberg Digits (SHD) spoken digit classification dataset with \(250\) time steps [69]. As summarized in Table ¿tbl:tab:32mnist?, PMSN model achieves the highest accuracy across all spiking neuron models with fewer or comparable amounts of parameters, demonstrating a superior capacity to capture temporal dependencies. Moreover, PMSN consistently achieves comparable or superior performance to advanced spiking architectures specifically designed for temporal processing, including temporal convolutions [56], [64], spike-based transformers [57], [65], [70]–[72], and SSMs [58], [59]. For the spike-based transformer baselines, we report publicly available results from prior works [56], [62], [72], [73]. These findings suggest that detailed modeling of neuronal structural complexity is an important and complementary dimension to architectural innovations in neuromorphic temporal processing.

We further evaluate our model’s ability to establish spatiotemporal and extended long-term temporal dependencies on the more challenging Sequential CIFAR-10 and CIFAR-100 tasks. For Sequential CIFAR-10, we explore two configurations: column-by-column scanning as per [43] (\(T=32\)), which evaluates the model’s capacity in integrating both spatial and temporal information, and pixel-by-pixel scanning (\(T=1024\)), which poses a greater challenge to learning long-term dependency. For Sequential CIFAR-100, we use the column-by-column configuration. To ensure a fair comparison, we employ the same network architecture for each individual task. As shown in Table 2, our PMSN surpasses the SOTA models by at least 2% accuracy, showcasing its superiority in multiscale temporal processing. As shown in Fig. 3, the learning curves of the PMSN model exhibit faster and more stable training convergence, aligning well with our theoretical analysis in Section 6.

Furthermore, we would like to stress the importance of the neuronal reset mechanism that has been neglected by many previous works [43]. As presented in Table 2, the accuracy consistently improves for both PMSN and LIF models after incorporating the reset mechanism. This is because the reset mechanism prevents the membrane potential from becoming excessively high, yielding a smoother distribution of membrane potentials across time that can facilitate stable information flow.

Beyond these vision-based and audio processing benchmarks, we further assess the scalability and practicality of PMSN on the more challenging UEA multivariate time-series classification benchmark [75] and the PPG-DaLiA time-series regression benchmark [78]. Both benchmarks involve significantly longer sequences, with several tasks containing more than 10,000 timesteps, and capture a broad range of real-world temporal processing scenarios, ranging from scientific time-series analysis to practical applications in healthcare. Following prior work, we adopt the same preprocessing and evaluation protocols to ensure a fair comparison [75], [77]. As reported in Tables 3 and 4, PMSN achieves competitive or superior performance compared with advanced non-spiking architectures that are widely used for temporal processing, including linear RNNs and SSMs, while being the only spike-based model among the compared methods. These results demonstrate that PMSN scales effectively to very long and challenging sequential tasks, underscoring its applicability to real-world temporal processing scenarios and its competitiveness beyond the spiking domain.

7.2 Scalability on Large-scale Image Classification↩︎

While the preceding evaluations primarily focus on temporal processing benchmarks, it is also important to examine whether the proposed spiking neuron model can scale to larger network architectures and datasets. To assess this scalability, we integrate PMSN into ResNet-18 and ResNet-34 architectures and evaluate its performance on the ImageNet-1K benchmark [79], with the results summarized in Table 5. The results show that PMSN achieves SOTA accuracy against other spiking models with identical architectures, reaching 70.98% on ImageNet-1K, a 3% absolute improvement over the IF baseline. This consistent improvement across both sequential benchmarks and large-scale image recognition tasks demonstrates the effectiveness and scalability of the proposed neuron model.

At the same time, it is worth noting that recent transformer-based architectures introduce innovations tailored for static image recognition and can achieve strong ImageNet performance [57], [65], [71]. Under comparable experimental settings, PMSN attains comparable performance to these approaches, despite relying solely on convolutional backbones and being primarily designed for temporal processing rather than static image feature extraction.

7.3 Visualization of Multi-compartment Dynamics↩︎

Having demonstrated the superior performance of PMSN, we next analyze how PMSN integrates multiscale temporal information through compartmental interactions. We first visualize the single-neuron impulse response in Fig.  4. After receiving an impulse at \(t=0\), different compartments exhibit damped oscillations with distinct frequencies and decay rates, characterized by \(\lambda_i=e^{a+ib}\). This reflects the system’s multiscale properties: compartments with higher oscillation frequencies correlate with faster timescales, while lower frequencies indicate slower scales [86]. This compartmental synergy enables the integration of information across various frequency domains and time spans, facilitating multiscale temporal processing.

We further examine the distribution of oscillation frequencies and damping coefficients across neurons in one layer. As shown in the middle and right panels of Fig. 4, PMSNs within the same population exhibit heterogeneous temporal dynamics, allowing the network to cover a broader spectrum of timescales.

7.4 Simulation Acceleration↩︎

To quantitatively evaluate the proposed acceleration techniques, we record the wall-clock forward- and backward-propagation times of PMSN using the GPU-enabled PyTorch library. Besides that, we also compare them against a range of serial and parallel models, including a serial MSN with neuronal dynamics identical to those of PMSN, LIF, PSN, and SPSN. Fig. 5 shows the acceleration ratios across diverse sequence lengths under the same single-hidden layer network architecture consisting of 256 neurons. Our parallel scheme shows a substantial speed-up over serial models, achieving training speed-up ratios up to \(134 \times\) and \(71 \times\) compared to the serial MSN and LIF models, respectively. When compared with SOTA parallel models, our PMSN outperforms SPSN but is slower than PSN. This can be attributed to a larger number of neuronal compartments used (i.e., five in PMSN vs. one in PSN) and less efficient implementation of FFT operations compared to PSN’s matrix multiplication in the current GPU acceleration framework. Nevertheless, this slight slowdown is worthwhile considering the overall effectiveness of our model. Additionally, we observe a positive correlation between the speed-up ratio and sequence length, indicating greater relative acceleration for longer sequences.

7.5 Theoretical Computational Cost Analysis↩︎

The PMSN offers an energy-efficient, hardware-friendly solution for practical applications. To systematically characterize these properties, we first analyze the spike sparsity and theoretical computational cost of PMSN in comparison with the PSN family [43] and conventional LIF neurons. All experiments are conducted on the pixel-by-pixel Sequential CIFAR-10 task using identical backbone architectures to ensure a fair comparison. To compare models across different scales, we vary the compartment number \(n\) for PMSN, the temporal receptive field size \(k\) for masked PSN and SPSN, and the network size for LIF neurons. The total cost is computed by considering both synaptic operations and neuronal operations. As shown in Fig. 6 (a), PMSN exhibits the lowest spike density among all compared models. This efficiency arises from two key mechanisms: firstly, the mutual inhibitory effect resulting from the coupling compartments may contribute to the spike frequency adaptation. Moreover, richer compartmental dynamics reduce the spike-based computational workload and suppress spike redundancy.

To estimate the theoretical computational cost, we count the number of Multiply–Accumulate (MAC) and Accumulate (AC) operations per inference. The analytical cost of each neuron model is provided in Supplementary Materials. Using CMOS energy estimates at 45 nm [87] (\(E_{AC}=0.9\) pJ, \(E_{MAC}=4.6\) pJ), Fig. 6 (b) shows that PMSN consistently achieves higher accuracy than other spiking neuron models at comparable computational costs.

7.6 Hardware Implementation and Energy Efficiency↩︎

To evaluate the practical deployability and efficiency of PMSN, we conduct an FPGA-based neuromorphic hardware implementation and compare its energy efficiency against conventional ANN implementations as well as a baseline SNN with LIF neurons. Details of the hardware platforms, software environments and the measurement pipeline are provided in Supplementary Materials.

Specifically, we develop a dedicated FPGA-based neuromorphic accelerator tailored for PMSN, as illustrated in Fig. 7. This hardware implementation preserves event-driven inter-neuron communication while retaining the rich internal neuronal dynamics of PMSN. Building upon this hardware deployment, we conduct an energy efficiency evaluation of PMSN on the PS-MNIST task and compare it against conventional ANN implementations executed on CPU and GPU platforms, and a baseline feedforward SNN with LIF neurons implemented on the same FPGA platform. The results are summarized in Table 6, where energy efficiency is defined as the ratio between throughput and power consumption. Compared to ANN implementations on CPU and GPU, PMSN achieves two orders of magnitude improvement in energy efficiency, demonstrating its suitability for real-time temporal signal processing in edge scenarios. Moreover, PMSN exhibits a more favorable accuracy-efficiency trade-off than baseline LIF neurons. In particular, a PMSN configuration with network size comparable to the LIF baseline achieves substantially higher accuracy at the cost of moderately reduced energy efficiency, while a more compact variant achieves markedly higher energy efficiency than the LIF baseline while still delivering significantly higher accuracy.

7.7 Ablation Study of Compartment Number↩︎

PMSN provides the flexibility to adjust the number of neuronal compartments according to task complexity and computational resources. To examine this property, we evaluate the effect of compartment number on Sequential CIFAR-10 and CIFAR-100 (\(T=32\)). As summarized in Table 7, increasing the number of compartments consistently improves accuracy on both datasets, indicating that larger compartment numbers provide additional temporal modeling capacity.

However, this performance gain comes at the cost of increased neuronal computation. As shown in Fig. 6 (b), increasing \(n\) from \(2\) to \(17\) substantially improves classification accuracy while incurring a more than 10-fold increase in theoretical computational cost, revealing a trade-off between model expressiveness and computational efficiency.

8 Discussion and Conclusion↩︎

In this work, we proposed a generalized multi-compartment neuron model with superior capacity in multiscale temporal processing. Furthermore, we introduced a parallel implementation for this model, enabling accelerated training on GPU-accelerated ML frameworks. Experimental results demonstrated its strong ability to establish temporal dependencies across multiple timescales, substantially accelerate training, scale to challenging tasks, and achieve a favorable trade-off between accuracy and energy efficiency.

While PMSN has demonstrated accurate and efficient neuromorphic temporal processing, how this neuron-level advance interacts with contemporary Transformer-based architectures remains an interesting question. Transformers provide architecture-level global representation learning, whereas PMSN enhances multiscale temporal integration within individual spiking units. Recent hybrid Transformer–SSM models suggest the benefit of combining global attention with structured temporal dynamics [88], [89]. Therefore, integrating PMSN with attention-based architectures remains a promising direction for future work.

Moreover, our FPGA-based neuromorphic accelerator evaluation provided hardware-oriented evidence for the potential deployability of PMSN on neuromorphic platforms, positioning it as a promising solution for accurate and energy-efficient real-time temporal processing in edge applications. Its compartmental dynamics can also be formulated as programmable state-space neuronal dynamics, suggesting potential compatibility with neuromorphic chips that support programmable neurons [90] and efficient state-space inference [91].

References↩︎