Neuronal populations in the brain exhibit remarkable heterogeneity, even among neurons of the same physiological class. Emerging evidence suggests that this diversity is not merely noise but a fundamental feature of neural computation and processing (Markram et al., 1997; Marder and Goaillard, 2006; Brémaud et al., 2007; Gast et al., 2024). This heterogeneity spans structural, genetic, environmental, and electrophysiological dimensions, including spiking thresholds (Gast et al., 2024; Huang et al., 2016), membrane time constants (Perez-Nieves et al., 2021; Wu et al., 2025), external currents (Chen and Campbell, 2022; Montbrió and Pazó, 2020; Montbrió et al., 2015), electrical coupling strengths (Parker, 2003) and reset mechanisms (Leng and Aihara, 2020). These variations reflect the heterogeneity in cellular composition and network organization across brain regions (Billeh et al., 2020; Hawrylycz et al., 2012).

Although this neural heterogeneity might appear detrimental to the reliability of neural networks, both theoretical and empirical studies have demonstrated that it can, in fact, enhance information encoding, learning robustness, and task-specific computation (Chelaru and Dragoi, 2008; Padmanabhan and Urban, 2010; Mejias and Longtin, 2012). For example, diversity in membrane and synaptic time constants improves generalization in learning models (Perez-Nieves et al., 2021), while differences in spiking thresholds allow SNNs to flexibly gate, encode, and transform signals (Gast et al., 2024). However, focusing on a single source offers an incomplete picture of how heterogeneity shapes neural computation. Given its multifaceted nature, understanding the interplay between different types of heterogeneity remains an important and largely uncharted direction.

Building on prior theoretical frameworks (Ly, 2015), we classify neural heterogeneity into three main categories: external, network, and intrinsic. External heterogeneity stems from variations in input currents, which are experimenter-controlled and reflect properties of sensory stimuli rather than network structure (Chen and Campbell, 2022; Montbrió and Pazó, 2020; Montbrió et al., 2015; Goldobin et al., 2021; Luccioli et al., 2019). Network heterogeneity refers to structural diversity within the network, such as variability in synaptic connectivity or electrical coupling strengths (Marder and Goaillard, 2006; Oswald et al., 2009). Intrinsic heterogeneity, in contrast, originates from neuron-specific properties that are independent of network interactions (Padmanabhan and Urban, 2010). This classification provides a useful lens to systematically evaluate their computational contributions. While individual forms of heterogeneity have been explored in isolation, it remains unclear whether their computational effects are consistent across different learning methods and task domains. Addressing this gap could yield a more unified understanding of how heterogeneity shapes computation in both biological and artificial neural systems.

Previous studies have mainly explored neural heterogeneity in the special task and learning methods. For example, Ridge Least Square (RLS) learning method has demonstrated performance gains in basic curve-fitting tasks when neuron spiking thresholds or neural time constants are varied (Gast et al., 2024; Wu et al., 2025). To extend beyond such constrained scenarios, recent work has incorporated two broader learning methods (Nicola and Clopath, 2017). First, originating in the field of reservoir computing, First-Order, Reduced and Controlled Error (FORCE) learning method leverages the dynamics of high-dimensional recurrent systems to perform computations (Sussillo and Abbott, 2009; Schliebs et al., 2011; Buonomano and Merzenich, 1995). Learning is driven by an external supervisory signal that provides error feedback. This method expands the potential applications for more complicated target functions, such as chaotic system prediction, songbird generation and memory recall. Second, Surrogate Gradient Descent (SGD) (Neftci et al., 2019) enables SNNs to classify vision (Orchard et al., 2015; Lichtsteiner et al., 2008) and auditory (Cramer et al., 2020) stimuli, and this approach is well-suited for integration into modern SNNs architectures. While both methods have expanded the learning capacity of SNNs, it remains unclear whether neural heterogeneity offers consistent benefits across this range of learning methods from simple regression to high-dimensional classification.

In this study, we investigate the general computational role of neural heterogeneity across three representative forms: external current, synaptic coupling strength, and partial reset. These forms have been widely used in modeling studies, yet their impact remains poorly understood. We evaluate their effects across three complementary learning methods: RLS, FORCE, and SGD, and apply them to a diverse set of tasks ranging from curve fitting and network reconstruction to real-world speech and image classification. Our results demonstrate that neural heterogeneity, regardless of its origin, consistently improves both performance and robustness in SNNs. These findings underscore its fundamental computational advantages and offer a unified perspective on the functional role of heterogeneity in neural systems. We argue that heterogeneity is not merely a biological artifact but a key design principle for efficient and adaptable learning in next-generation SNNs. This work thus highlights a potentially unifying computational role of neural heterogeneity across multiple SNNs learning methods.

The Izhikevich (IK) model is a reduced form of the Hodgkin-Huxley-type neuron model, derived via bifurcation analysis (Izhikevich, 2007). Despite its reduced complexity, it preserves essential features of neuronal excitability through an adaptive quadratic mechanism. In this study, we incorporate external, network, and intrinsic heterogeneity by allowing three model parameters to vary across neurons.

External heterogeneity is introduced by assigning each neuron a unique external current, a method commonly employed in previous studies (Chen and Campbell, 2022; Montbrió and Pazó, 2020; Montbrió et al., 2015). This design mimics the inherent variability in sensory input processing in biological systems: sensory neurons receive non-uniform external drives due to differences in receptor sensitivity. For example, variations in the density of light-sensitive opsins in retinal cells or mechanosensitive channels in auditory hair cells lead to divergent responses to the same stimulus, which is well-documented in neuroscience (Chen and Campbell, 2022; Montbrió et al., 2015)

Network heterogeneity is captured by scaling the electrical coupling strengths, such that neurons receive varying degrees of network influence (Xue et al., 2014). This approach reflects extensive experimental evidence of heterogeneity in synaptic and electrical coupling (Marder and Goaillard, 2006; Oswald et al., 2009; Parker, 2003). By incorporating g_i heterogeneity our model replicates this structural diversity ensuring that network interactions align with the non-uniform connectivity patterns of real brains.

Intrinsic heterogeneity is implemented via the partial reset mechanism following spike generation, representing incomplete membrane repolarization (Leng and Aihara, 2020; Rospars and Lánskỳ, 1993). This mechanism allows neuron-specific post-spike dynamics, enhancing variability that arises intrinsically from the individual neuron's state. Similar approaches have been used to account for stochasticity and diversity in spiking behaviors (Bugmann et al., 1997; Kirst et al., 2009).

The network model for a population of coupled IK neurons is described by the following equations (Chen and Campbell, 2022)

where v_i denotes the membrane potential of the ith neuron and w_i its associated recovery variable, which governs spike-frequency adaptation. The input current is modeled as I_i = η_i + I_ext(t), where η_i represents neuron-specific external input and I_ext(t) is a global time-varying external drive. This formulation introduces external heterogeneity through the distributed input η_i. Network heterogeneity is incorporated by assigning electrical coupling strengths g_i, which scales the strength of recurrent input each neuron receives. The core biophysical parameters of each neuron include membrane capacitance C, the leakage parameter k, the resting potential v_r, the spike threshold potential v_t, the synaptic reversal potential E, the recovery variable time constant τ_w, and the scaling factor β. When v_i reaches the spike cutoff v_peak, it is reset according to a partial reset rule:

where 0 < θ ≤ 1 is the reset coefficient for neuron i, introducing intrinsic heterogeneity. A value of θ_i = 0 corresponds to a full reset, while θ_i > 0 retains partial voltage history. Simultaneously, the adaptation variable undergoes a discrete jump by an amount w_jump. The synaptic gating variable s lies between 0 and 1, representing the synaptic activation of each neuron. For analytical simplicity, we assume an all-to-all connectivity structure, which is a standard approximation widely used in prior studies (Montbrió et al., 2015; Byrne et al., 2020). The temporal dynamics of synaptic transmission are formally described by a first-order differential equation (Ermentrout and Terman, 2010)

where δ(t) is the Dirac delta function, and tjk represents the time of the kth spike of the jth neuron. τ_s is a decay time constant, and s_jump is the coupling strength. We consider the network system of Equations 1–3 as the basic components of the following SNNs.

To introduce different sources of heterogeneity, we assume that the neuron-specific parameters {η_i, g_i, θ_i} are distributed according to a Lorentzian probability distribution parameterized by a center and half-width-at-half-maximum, as described by the density function

where * = η, g, θ, denoting the three sources of heterogeneity respectively.

Here, Δ^* denotes the half-width-at-half-maximum of the Lorentzian distribution, which controls the spread of heterogeneity. Larger Δ^* values correspond to greater variability in the parameter. This means that Δ_η controls the range of neuron-specific external currents, Δ_g controls the spread of coupling strengths, and Δ_θ controls the variability of partial reset coefficients).

The Lorentzian distribution is chosen not only for mathematical tractability (facilitating mean-field analysis in Appendix A) but also for biological relevance: it captures the heavy-tailed variability of sensory input drives observed in empirical studies unlike Gaussian distributions that underestimate the range of real-world sensory neuron responses (Gast et al., 2024; Perez-Nieves et al., 2021).

We investigate the role of heterogeneity in SNNs trained using methods ranging from traditional machine learning to state-of-the-art deep learning, considering three perspectives: external input, network structure, and intrinsic neuronal variability. Our findings consistently show that heterogeneity facilitates learning and enhances computational performance. Furthermore, the generalizability of its benefits across tasks and models supports its role as a fundamental mechanism in neural computation. These results offer promising directions for advancing biologically inspired neural system design.

A key mechanistic underpinning of heterogeneity's performance gains lies in its ability to enhance the dynamical richness of the network, as implied by the mean-field bifurcation analysis in Appendix A. The mean-field model reveals how the population firing rate transitions across distinct dynamical regimes, as the mean external current (η̄) varies (Figure 1a). By introducing heterogeneity in external current (η_i), coupling strength (g_i), and partial reset (θ_i), we expand the network's capacity to explore these regimes. For instance, in reservoir computing (Section 3.1), heterogeneity enables the network to move beyond the quiescent regime (where low firing rates hinder information encoding) and stably operate in oscillatory or persistently active regimes, stating that support robust temporal feature encoding. This dynamical expansion is critical for tasks requiring precise temporal processing, such as fitting complex sinusoidal signals (Section 3.1) and classifying time-varying stimuli in the DVS128, SHD, and SSC datasets (Section 3.3). Without heterogeneity, homogeneous networks remain confined to narrow dynamical ranges, limiting their ability to adapt to diverse input patterns.

Another theoretical basis for heterogeneity's benefits is its role in reducing sensitivity to hyperparameters, particularly evident in FORCE learning (Section 3.2). Homogeneous FORCE networks exhibit fragile performance, relying on tightly tuned hyperparameters G (gain of static recurrent weights) and Q (coefficient of learned weights) to maintain functional dynamics (Figure 4a). In contrast, introducing any single form of heterogeneity, whether in external current (Δ_η), coupling strength (Δ_g), or partial reset (θ), broadens the viable (G, Q) parameter space. For example, external current heterogeneity extends the lower boundary of G from 3,000 to 1,000 in certain configurations, while combining all three heterogeneities fully restores performance across the entire tested range of G and Q in Figures 4b–f. This hyperparameter robustness can be attributed to heterogeneity's ability to break symmetry in homogeneous networks: by introducing variability in neuron-specific properties, it reduces redundant neural activity and creates more distinct input-output mapping pathways. This symmetry breaking allows the network to efficiently utilize its high-dimensional dynamics, minimizing the risk of performance collapse due to slight hyperparameter deviations.

Our model's accuracy on real-world datasets such as DVS128 and SHD (e.g., a maximum accuracy of 73.6 ± 0.7% on DVS128, as shown in Table 1) has not yet reached that of state-of-the-art (SOTA) implementations. A major contributing factor is that our framework is built upon the Izhikevich neuron model, which, while providing a more biologically realistic description of spiking dynamics through the recovery variable (w_i) that captures phenomena such as spike-frequency adaptation, has not been widely adopted in current deep SNN pipelines due to its higher computational cost and less optimized parameterization. As a result, our experiments focus primarily on establishing the fundamental computational role of neural heterogeneity, rather than achieving peak task-specific performance. Future work will incorporate our proposed principles into more advanced, SOTA architectures (e.g., hybrid convolutional SNNs or transformer-like event encoders) to further validate the generality of heterogeneity under high-performance settings.

Our findings align with and extend prior hypotheses about neural heterogeneity in biological systems. On one hand, we confirm the long-standing hypothesis that heterogeneity enhances learning robustness, which is a conclusion previously restricted to LIF networks and single heterogeneity types (Perez-Nieves et al., 2021; Mejias and Longtin, 2012). Our work broadens this insight by demonstrating that such robustness is universal across three distinct heterogeneity forms, three learning methods, and a variety of tasks from simple curve fitting to complex real-world classification. On the other hand, we also reveal that the relationship between the number and type of heterogeneity sources and performance is nonlinear: combined heterogeneity (e.g., joint variation in η, g, and θ) does not always outperform single-type heterogeneity. For instance, on the SHD dataset, independent heterogeneity in η alone yielded higher accuracy than the combined use of η and g. This finding suggests that biological systems may employ selective rather than accumulative heterogeneity, prioritizing specific forms according to task demands.

Nevertheless, our current results, obtained on moderately scaled networks, are sufficient to establish heterogeneity as a general and biologically plausible computational mechanism. We acknowledge that further Validation on more complex architectures and neuron models is necessary to fully assess the scalability and task-specific impact of heterogeneity. In future work, we aim to extend our framework toward SOTA network designs and larger datasets to bridge the gap between conceptual generality and high-performance realization.

In conclusion, our study provides systematic evidence that external, network, and intrinsic heterogeneity collectively enhance the learning capacity and robustness of SNNs. By linking these empirical findings to dynamical systems theory and hyperparameter sensitivity analysis, we strengthen the theoretical basis for heterogeneity as a core design principle, which bridges biological observations of neural diversity and engineering goals of optimizing SNN performance. This unifying perspective advances both our understanding of biological computation and our ability to construct more adaptive, efficient artificial neural systems.