# Heterogeneous Neuronal and Synaptic Dynamics for Spike-Efficient
Unsupervised Learning: Theory and Design Principles

Biswadeep Chakraborty & Saibal Mukhopadhyay
Department of Electrical and Computer Engineering
Georgia Institute of Technology
<EMAIL_ADDRESS>

###### Abstract

This paper shows that the heterogeneity in neuronal and synaptic dynamics
reduces the spiking activity of a Recurrent Spiking Neural Network (RSNN)
while improving prediction performance, enabling spike-efficient
(unsupervised) learning. We analytically show that the diversity in neurons’
integration/relaxation dynamics improves an RSNN’s ability to learn more
distinct input patterns (higher memory capacity), leading to improved
classification and prediction performance. We further prove that heterogeneous
Spike-Timing-Dependent-Plasticity (STDP) dynamics of synapses reduce spiking
activity but preserve memory capacity. The analytical results motivate
Heterogeneous RSNN design using Bayesian optimization to determine
heterogeneity in neurons and synapses to improve $\mathcal{E}$, defined as the
ratio of spiking activity and memory capacity. The empirical results on time
series classification and prediction tasks show that optimized HRSNN increases
performance and reduces spiking activity compared to a homogeneous RSNN.

###  1 Introduction

Spiking neural networks (SNNs) (Ponulak & Kasinski, 2011) use unsupervised
bio-inspired neurons and synaptic connections, trainable with either
biological learning rules such as spike-timing-dependent plasticity (STDP)
(Gerstner & Kistler, 2002; Chakraborty & Mukhopadhyay, 2023a) or supervised
statistical learning algorithms such as surrogate gradient (Neftci et al.,
2019). Empirical results on standard SNNs also show good performance for
various tasks Chakraborty et al. (2023), including spatiotemporal data
classification, (Lee et al., 2017; Khoei et al., 2020), sequence-to-sequence
mapping (Zhang & Li, 2020), object detection (Chakraborty et al., 2021; Kim et
al., 2020), and universal function approximation (Gelenbe et al., 1999;
Iannella & Back, 2001). An important motivation for the application of SNN in
machine learning (ML) is the sparsity in the firing (activation) of the
neurons, which reduces energy dissipation during inference (Wu et al., 2019).
Many prior works have empirically shown that SNN has lower firing activity
than artificial neural networks and can improve energy efficiency (Kim et al.,
2022; Srinivasan & Roy, 2019). However, there are very few analytical studies
on how to reduce the spiking activity of an SNN while maintaining its learning
performance. Understanding and optimizing the relations between spiking
activity and performance will be key to designing energy-efficient SNNs for
complex ML tasks.

In this paper, we derive analytical results and present design principles from
optimizing the spiking activity of a recurrent SNN (RSNN) while maintaining
prediction performance. Most SNN research in ML considers a simplified network
model with a homogeneous population of neurons and synapses (homogeneous RSNN
(MRSNN)) where all neurons have uniform integration/relaxation dynamics, and
all synapses use the same long-term potentiation (LTP) and long-term
depression (LTD) dynamics in STDP learning rules. On the contrary,
neurobiological studies have shown that a brain has a wide variety of neurons
and synapses with varying firing and plasticity dynamics, respectively
(Destexhe & Marder, 2004; Gouwens et al., 2019; Hansel et al., 1995; Prescott
et al., 2008). We show that optimizing neuronal and synaptic heterogeneity
will be key to simultaneously reducing spiking activity while improving
performance.

We define the spike efficiency $\mathcal{E}$ of an RSNN as the ratio of its
memory capacity $\mathcal{C}$ and average spiking activity $\tilde{S}$. Given
a fixed number of neurons and synapses, a higher $\mathcal{C}$ implies a
network can learn more patterns and hence, perform better in classification or
prediction tasks (Aceituno et al., 2020; Goldmann et al., 2020); a lower
spiking rate implies that a network is less active, and hence, will consume
less energy while making inferences (Sorbaro et al., 2020; Rathi et al.,
2021). We analytically show that a Heterogeneous Recurrent SNN (HRSNN) model
leads to a more spike-efficient learning architecture by reducing spiking
activity while improving $\mathcal{C}$ (i.e., performance) of the learning
models. In particular, we make the following contributions to the theoretical
understanding of an HRSNN.

  * •

We prove that for a finite number of neurons, models with heterogeneity among
the neuronal dynamics has higher memory capacity $\mathcal{C}$.

  * •

We prove that heterogeneity in the synaptic dynamics reduces the spiking
activity of neurons while maintaining $\mathcal{C}$. Hence, a model with
heterogeneous synaptic dynamics has a lesser firing rate than a model with
homogeneous synaptic dynamics.

  * •

We connect the preceding results to prove that simultaneously using
heterogeneity in neurons and synapses, as in an HRSNN, improves the spike
efficiency of a network.

We empirically characterize HRSNN considering the tasks of (a) classifying
time series ( Spoken Heidelberg Digits (SHD)) and (b) predicting the evolution
of a dynamical system (a modified chaotic Lorenz system). The theoretical
results are used to develop an HRSNN architecture where a modified Bayesian
Optimization (BO) is used to determine the optimal distribution of neuron and
synaptic parameters to maximize $\mathcal{E}$. HRSNN exhibits a better
performance (higher classification accuracy and lower NRMSE loss) with a
lesser average spike count $\tilde{S}$ than MRSNN.

Related Works Inspired by the biological observations, recent empirical
studies showed potential for improving SNN performance with heterogeneous
neuron dynamics(Perez-Nieves et al., 2021; Chakraborty & Mukhopadhyay, 2023b).
However, there is a lack of theoretical understanding of why heterogeneity
improves SNN performance, which is critical for optimizing SNNs for complex
tasks. She et al. (2022) have analytically studied the universal sequence
approximation capabilities of a feed-forward network of neurons with varying
dynamics. However, they did not consider heterogeneity in plasticity dynamics,
and the results are applicable only for a feed-forward SNN and do not extend
to recurrent SNNs (RSNN). The recurrence is not only a fundamental component
of a biological brain (Soures & Kudithipudi, 2019), but as a machine learning
(ML) model, RSNN also shows good performance in modeling spatiotemporal and
nonlinear dynamics (Pyle & Rosenbaum, 2017; Gilra & Gerstner, 2017). Hence, it
is critical to understand whether heterogeneity can improve learning in an
RSNN. To the best of our knowledge, this is the first work that analytically
studies the impact of heterogeneity in synaptic and neuronal dynamics in an
RSNN. This work shows that only using neuronal heterogeneity improves
performance and does not impact spiking activity. The number of spikes
required for the computation increases exponentially with the number of
neurons. Therefore, simultaneously analyzing and optimizing neuronal and
synaptic heterogeneity, as demonstrated in this work, is critical to design an
energy-efficient recurrent SNN.

###  2 Preliminaries and Definitions

Table 1: Notations Table

Notation | Meaning | Notation | Meaning | Notation | Meaning | Notation | Meaning
---|---|---|---|---|---|---|---
HRSNN |  | Heterogeneous
---
Recurrent Spiking
Neural Network
$S_{i}(t)$ |  | Spike train from
---
neuron i at time t
$z_{j}^{t}$ |  | Spike indicator
---
Function
$\mathbf{r}(t)$ | States of RSNN
MRSNN |  | Homogeneous
---
Recurrent Spiking
Neural Network
$\Delta w$ |  | Synaptic Weight
---
Update
$v_{th}$ | Threshold voltage | $\mathbf{\Sigma}$ |  | Covariance Matrix of
---
$\mathbf{r}(t)$
$\tau_{m}$ |  | Membrane
---
Time Constant
$\tau_{\pm}$ |  | Synaptic Time
---
Constants
$\mathcal{C}$ | Memory Capacity | $x_{i}$ | Model Input
$v_{i}(t)$ |  | Membrane
---
Potential
$\eta_{\pm}$ | Scaling Functions | $y_{\tau}$ |  | Model Output
---
with delay $\tau$
$\mathcal{E}$ | Spike Efficiency
$v_{rest}$ | Resting potential | $\alpha$ |  | Constant
---
Decay Factor
$\mathcal{R}$ |  | Recurrent Layer
---
in RSNN
$S_{R}$, $S_{M}$ |  | Total Spike Count
---
HRSNN (R) , MRSNN (M)
$\tilde{S}_{R},\tilde{S}_{M}$ |  | Avg. Spike Count
---
per neuron of
HRSNN (R), MRSNN (M)
$W_{ji}$ |  | Synapse weight
---
from neuron $i$
to neuron $j$
$\mathcal{H}$ |  | Heterogeneity
---
in Neuronal
Parameters
$\Phi_{R}$, $\Phi_{M}$ |  | Spike frequency
---
of HRSNN (R),
MRSNN (M)

We now define the key terms used in the paper. Table 1 summarizes the key
notations used in this paper. Figure 1 shows the general structure of the
HRSNN model with heterogeneity in both the LIF neurons and the STDP dynamics.
It is to be noted here that there are a few assumptions we use for the rest of
the paper: Firstly, the heterogeneous network hyperparameters are estimated
before the training and inference. The hyperparameters are frozen after
estimation and do not change during the model evaluation. Secondly, this paper
introduces neuronal and synaptic dynamics heterogeneity by using a
distribution of specific parameters. However, other parameters can also be
chosen, which might lead to more interesting/better performance or
characteristics. We assume a mean-field model where the synaptic weights
converge for the analytical proofs. In addition, it must be noted that LIF
neurons have been shown to demonstrate different states Brunel (2000). Hence,
for the analytical study of the network, we use the mean-field theory to
analyze the collective behavior of a dynamical system comprising many
interacting particles.

Figure 1: Concept of HRSNN with variable Neuronal and Synaptic Dynamics

Heterogeneous LIF Neurons We use the Leaky Integrate and Fire (LIF) neuron
model in all our simulations. In this model, the membrane potential of the
$i$-th neuron $u_{i}(t)$ varies over time as:

$\tau_{m}\frac{dv_{i}(t)}{dt}=-\left(v_{i}(t)-v_{rest}\right)+I_{i}(t)$ (1)

where $\tau_{\mathrm{m}}$ is the membrane time constant, $v_{rest}$ is the
resting potential and $I_{\mathrm{i}}$ is the input current. When the membrane
potential reaches the threshold value $v_{\text{th }}$ a spike is emitted,
$v_{i}(t)$ resets to the reset potential $v_{\mathrm{r}}$ and then enters a
refractory period where the neuron cannot spike. Spikes emitted by the $j$ th
neuron at a finite set of times $\left\\{t_{j}\right\\}$ can be formalized as
a spike train $\displaystyle S_{i}(t)=\sum\delta\left(t-t_{i}\right)$.

Let the recurrent layer of an RSNN be $\mathcal{R}$. We incorporate
heterogeneity in the LIF neurons by using different membrane time constants
$\tau_{m,i}$ and threshold voltages $v_{th,i}$ for each LIF neuron $i$ in
$\mathcal{R}$. This gives a distribution of time constants and threshold
voltages of the LIF neurons in $\mathcal{R}$.

Heterogeneous STDP:

The STDP rule for updating a synaptic weight ($\Delta w$) is defined by Pool &
Mato (2011):

$\Delta w(\Delta t)=\left\\{\begin{array}[]{l}A_{+}(w)e^{-\frac{|\Delta
t|}{\tau_{+}}}\text{ if }\Delta t\geq 0\\\ -A_{-}(w)e^{-\frac{|\Delta
t|}{\tau_{-}}}\text{ if }\Delta t<0\end{array}\right.\quad
s.t.,A_{+}(w)=\eta_{+}\left(w_{\max}-w\right),A_{-}(w)=\eta_{-}\left(w-w_{\min}\right)$
(2)

where $\Delta t=t_{\text{post }}-t_{\text{pre }}$ is the time difference
between the post-synaptic spike and the pre-synaptic one, with synaptic time-
constant $\tau_{\pm}$. In heterogeneous STDP, we use an ensemble of values
from a distribution for $\tau_{\pm}$ and the scaling functions $\eta_{\pm}$.

Heterogeneity:

We define heterogeneity as a measure of the variability of the hyperparameters
in an RSNN that gives rise to an ensemble of neuronal dynamics.

Entropy is used to measure population diversity. Assuming that the random
variable for the hyperparameters $X$ follows a multivariate Gaussian
Distribution ($X\sim\mathcal{N}(\mu,\Sigma)$), then the differential entropy
of $x$ on the multivariate Gaussian distribution, is
$\displaystyle\mathcal{H}(x)=\frac{n}{2}\ln(2\pi)+\frac{1}{2}\ln|\Sigma|+\frac{n}{2}$.
Now, if we take any density function $q(\mathrm{x})$ that satisfies
$\displaystyle\int q(\mathrm{x})x_{i}x_{j}d\mathrm{x}=\Sigma_{ij}$ and
$p=\mathcal{N}(0,\Sigma),$ then $\mathcal{H}(q)\leq\mathcal{H}(p).$ (Proof in
Suppl. Sec. A) The Gaussian distribution maximizes the entropy for a given
covariance. Hence, the log-determinant of the covariance matrix bounds
entropy. Thus, for the rest of the paper, we use the determinant of the
covariance matrix to measure the heterogeneity of the network.

Memory Capacity: Given an input signal x(t), the memory capacity $\mathcal{C}$
of a trained RSNN model is defined as a measure for the ability of the model
to store and recall previous inputs fed into the network (Jaeger, 2001; Jaeger
et al., 2001).
In this paper, we use $\mathcal{C}$ as a measure of the performance of the
model, which is based on the network’s ability to retrieve past information
(for various delays) from the reservoir using the linear combinations of
reservoir unit activations observed at the output. Intuitively, HRSNN can be
interpreted as a set of coupled filters that extract features from the input
signal. The final readout selects the right combination of those features for
classification or prediction. First, the $\tau$-delay $\mathcal{C}$ measures
the performance of the $\mathrm{RC}$ for the task of reconstructing the
delayed version of model input $x(t)$ at delay $\tau$ (i.e., $x(t-\tau)$ ) and
is defined as the squared correlation coefficient between the desired output (
$\tau$-time-step delayed input signal, $x(t-\tau))$ and the observed network
output $y_{\tau}(t)$, given as:

$\mathcal{C}=\lim_{\tau_{\max}\rightarrow\infty}\sum_{\tau=1}^{\tau_{\max}}\mathcal{C}(\tau)=\lim_{\tau_{\max}\rightarrow\infty}\sum_{\tau=1}^{\tau_{\max}}\frac{\operatorname{Cov}^{2}\left(x(t-\tau),y_{\tau}(t)\right)}{\operatorname{Var}(x(t))\operatorname{Var}\left(y_{\tau}(t)\right)},\tau\in\mathbb{N},$
(3)

where $\operatorname{Cov}(\cdot)$ and $\operatorname{Var}(\cdot)$ denote the
covariance function and variance function, respectively. The $y_{\tau}(t)$ is
the model output in this reconstruction task. $\mathcal{C}$ measures the
ability of RSNN to reconstruct precisely the past information of the model
input. Thus, increasing $\mathcal{C}$ indicates the network is capable of
learning a greater number of past input patterns, which in turn, helps in
increasing the performance of the model. For the simulations, we use
$\tau_{\max}=100$.

Spike-Efficiency:  Given an input signal x(t), the spike-efficiency
($\mathcal{E}$) of a trained RSNN model is defined as the ratio of the memory
capacity $\mathcal{C}$ to the average total spike count per neuron
$\tilde{S}$.

$\mathcal{E}$ is an analytical measure used to compare how $\mathcal{C}$ and
hence the model’s performance is improved with per unit spike activity in the
model. Ideally, we want to design a system with high $\mathcal{C}$ using fewer
spikes. Hence we define $\mathcal{E}$ as the ratio of the memory capacity
using $N_{\mathcal{R}}$ neurons $\mathcal{C}(N_{\mathcal{R}})$ to the average
number of spike activations per neuron ($\tilde{S}$) and is given as:

${\color[rgb]{0,0,0}\mathcal{E}=\frac{\mathcal{C}(N_{\mathcal{R}})}{\frac{\sum_{i=1}^{N_{\mathcal{R}}}S_{i}}{N_{\mathcal{R}}}}},\quad\quad
S_{i}=\int_{0}^{T}s_{i}(t)dt\approx N_{\text{post
}}\frac{T}{\int_{t_{ref}}^{\infty}t\Phi_{i}dt}$ (4)

where $N_{\text{post }}$ is the number of postsynaptic neurons, $\Phi_{i}$ is
the inter-spike interval spike frequency for neuron $i$, and $T$ is the total
time. It is to be noted here that the total spike count $S$ is obtained by
counting the total number of spikes in all the neurons in the recurrent layer
until the emission of the first spike at the readout layer.

###  3 Heterogeneous RSNN: Analytical Results

We present three main analytical findings. Firstly, neuronal dynamic
heterogeneity increases memory capacity by capturing more principal components
from the input space, leading to better performance and improved
$\mathcal{C}$. Secondly, STDP dynamic heterogeneity decreases spike activation
without affecting $\mathcal{C}$, providing better orthogonalization among the
recurrent network states and a more efficient representation of the input
space, lowering higher-order correlation in spike trains. This makes the model
more spike-efficient since the higher-order correlation progressively
decreases the information available through neural population (Montani et al.,
2009; Abbott & Dayan, 1999). Finally, incorporating heterogeneity in both
neuron and STDP dynamics boosts the $\mathcal{C}$ to spike activity ratio,
i.e., $\mathcal{E}$, which enhances performance while reducing spike counts.

Memory Capacity: The performance of an RSNN depends on its ability to retain
the memory of previous inputs. To quantify the relationship between the
recurrent layer dynamics and $\mathcal{C}$, we note that extracting
information from the recurrent layer is made using a combination of the
neuronal states. Hence, more linearly independent neurons would offer more
variable states and, thus, more extended memory.

Lemma 3.1.1: The state of the neuron can be written as follows: $\displaystyle
r_{i}(t)=\sum_{k=0}^{N_{R}}\sum_{n=1}^{N_{R}}\lambda_{n}^{k}\left\langle
v_{n}^{-1},\mathbf{w}^{\mathrm{in}}\right\rangle\left(v_{n}\right)_{i}x(t-k)$,
where $\mathbf{v}_{n},\mathbf{v}_{n}^{-1}\in\mathbf{V}$ are, respectively, the
left and right eigenvectors of $\mathbf{W}$,  $\mathbf{w}^{\text{in}}$ are the
input weights, and $\lambda_{n}^{k}\in\mathbf{\lambda}$ belongs to the
diagonal matrix containing the eigenvalues of $\mathbf{W}$;
$\mathbf{a}_{i}=\left[a_{i,0},a_{i,1},\ldots\right]$ represents the
coefficients that the previous inputs $\mathbf{x}_{t}=[x(t),x(t-1),\ldots]$
have on $r_{i}(t)$.

Short Proof: (See Suppl. Sec. B for full proof)  As discussed by Aceituno et
al. (2020), the state of the neuron can be represented as
$\mathbf{r}(t)=\mathbf{W}\mathbf{r}(t-1)+\mathbf{w}^{\text{in }}x(t)$, where
$\mathbf{w}^{\text{in }}$ are the input weights. We can simplify this using
the coefficients of the previous inputs and plug this term into the covariance
between two neurons. Hence, writing the input coefficients $\mathbf{a}$ as a
function of the eigenvalues of $\mathbf{W}$,

$\mathbf{r}(t)=\sum_{k=0}^{\infty}\mathbf{W}^{k}\mathbf{w}^{\mathrm{in}}x(t-k)=\sum_{k=0}^{\infty}\left(\mathbf{V}\mathbf{\Lambda}^{k}\mathbf{V}^{-1}\right)\mathbf{w}^{\mathrm{in}}x(t-k)\Rightarrow
r_{i}(t)=\sum_{k=0}^{N_{R}}\sum_{n=1}^{N_{R}}\lambda_{n}^{k}\left\langle
v_{n}^{-1},\mathbf{w}^{\mathrm{in}}\right\rangle\left(v_{n}\right)_{i}x(t-k)\blacksquare$

Theorem 1:  If the memory capacity of the HRSNN and MRSNN networks are denoted
by $\mathcal{C}_{H}$ and $\mathcal{C}_{M}$ respectively, then,
$\displaystyle\mathcal{C}_{H}\geq\mathcal{C}_{M}$, where the heterogeneity in
the neuronal parameters $\mathcal{H}$ varies inversely to the correlation
among the neuronal states measured as
$\sum_{n=1}^{N_{\mathcal{R}}}\sum_{m=1}^{N_{\mathcal{R}}}\operatorname{Cov}^{2}\left(x_{n}(t),x_{m}(t)\right)$
which in turn varies inversely with $\mathcal{C}$.

Intuitive Proof:  (See Suppl. Sec. B for full proof)  Aceituno et al. (2020)
showed that the $\mathcal{C}$ increases when the variance along the
projections of the input into the recurrent layer are uniformly distributed.
We show that this can be achieved efficiently by using heterogeneity in the
LIF dynamics. More formally, let us express the projection in terms of the
state space of the recurrent layer. We show that the raw variance in the
neuronal states $\mathcal{J}$ can be written as
$\displaystyle\mathcal{J}=\frac{\sum_{n=1}^{N_{\mathcal{R}}}\mathbf{\lambda}_{n}^{2}(\mathbf{\mathbf{\Sigma}})}{\left(\sum_{n=1}^{N_{\mathcal{R}}}\mathbf{\lambda}_{n}(\mathbf{\mathbf{\Sigma}})\right)^{2}}$
where $\lambda_{n}(\mathbf{\Sigma})$ is the $n$th eigenvalue of
$\mathbf{\Sigma}$. We further show that with higher $\mathcal{H}$, the
magnitude of the eigenvalues of $\mathbf{W}$ decreases and hence leads to a
higher $\mathcal{J}$. Now, we project the inputs into orthogonal directions of
the network state space and model the system as
$\displaystyle\mathbf{r}(t)=\sum_{\tau=1}^{\infty}\mathbf{a}_{\tau}x(t-\tau)+\varepsilon_{r}(t)$
where the vectors $\mathbf{a}_{\tau}\in\mathbb{R}^{N}$ are correspond to the
linearly extractable effect of $x(t-\tau)$ onto $\mathbf{r}(t)$ and
$\varepsilon_{r}(t)$ is the nonlinear contribution of all the inputs onto the
state of $\mathbf{r}(t)$. First, we show that $\mathcal{C}$ increases when the
variance along the projections of the input into the recurrent layer is more
uniform. Intuitively, the variances at directions $\mathbf{a}_{\tau}$ must fit
into the variances of the state space, and since the projections are
orthogonal, the variances must be along orthogonal directions. Hence, we show
that increasing the correlation among the neuronal states increases the
variance of the eigenvalues, which would decrease our memory bound
$\mathcal{C}^{*}$. We show that heterogeneity is inversely proportional to
$\displaystyle\sum_{n=1}^{N_{\mathcal{R}}}\operatorname{Cov}^{2}\left(x_{n}(t),x_{m}(t)\right)$.
We see that increasing the correlations between neuronal states decreases the
heterogeneity of the eigenvalues, which reduces $\mathcal{C}$. We show that
the variance in the neuronal states is bounded by the determinant of the
covariance between the states; hence, covariance increases when the neurons
become correlated. As $\mathcal{H}$ increases, neuronal correlation decreases.
Aceituno et al. (2020) proved that the neuronal state correlation is inversely
related to $\mathcal{C}$. Hence, for HRSNN, with $\mathcal{H}>0$,
$\displaystyle\mathcal{C}_{H}\geq\mathcal{C}_{M}$. $\blacksquare$

Spiking Efficiency We analytically prove that the average firing rate of HRSNN
is lesser than the average firing rate of the MRSNN model by considering a
subnetwork of the HRSNN network and modeling the pre-and post-synaptic spike
trains using a nonlinear interactive Hawkes process with inhibition, as
discussed by Duval et al. (2022). The details of the model are discussed in
Suppl. Sec. B.

Lemma 3.2.1:  If the neuronal firing rate of the HRSNN network with only
heterogeneity in LTP/LTD dynamics of STDP is represented as $\Phi_{R}$ and
that of MRSNN represented as $\Phi_{M}$, then the HRSNN model promotes
sparsity in the neural firing which can be represented as
$\displaystyle\Phi_{R}<\Phi_{M}$.

Short Proof: (See Suppl. Sec. B for full proof) We show that the average
firing rate of the model with heterogeneous STDP (LTP/LTD) dynamics (averaged
over the population of neurons) is lesser than the corresponding average
neuronal activation rate for a model with homogeneous STDP dynamics. We prove
this by taking a sub-network of the HRSNN model. Now, we model the input spike
trains of the pre-synaptic neurons using a multivariate interactive, nonlinear
Hawkes process with multiplicative inhibition. Let us consider a population of
neurons of size $N$ that is divided into population $A$ (excitatory) and
population $B$ (inhibitory). We use a particular instance of the model given
in terms of a family of counting processes
$\left(Z_{t}^{1},\ldots,Z_{t}^{N_{A}}\right.$ ) (population $\left.A\right)$
and $\left(Z_{t}^{N_{A}+1},\ldots,Z_{t}^{N}\right)$ (population $B$ ) with
coupled conditional stochastic intensities given respectively by $\lambda^{A}$
and $\lambda^{B}$ as follows:

$\displaystyle\lambda_{t}^{A,N}:=\Phi_{A}\left(\frac{1}{N}\sum_{j\in
A}\int_{0}^{t^{-}}h_{1}(t-u)\mathrm{d}Z_{u}^{j}\right)\Phi_{B\rightarrow
A}\left(\frac{1}{N}\sum_{j\in
B}\int_{0}^{t^{-}}h_{2}(t-u)\mathrm{d}Z_{u}^{j}\right)$
$\displaystyle\lambda_{t}^{B,N}:=\Phi_{B}\left(\frac{1}{N}\sum_{j\in
B}\int_{0}^{t^{-}}h_{3}(t-u)\mathrm{d}Z_{u}^{j}\right)+\Phi_{A\rightarrow
B}\left(\frac{1}{N}\sum_{j\in
A}\int_{0}^{t^{-}}h_{4}(t-u)\mathrm{d}Z_{u}^{j}\right)$ (5)

where $A,B$ are the populations of the excitatory and inhibitory neurons,
respectively, $\lambda_{t}^{i}$ is the intensity of neuron $i,\Phi_{i}$ a
positive function denoting the firing rate, and $h_{j\rightarrow i}(t)$ is the
synaptic kernel associated with the synapse between neurons $j$ and $i$.
Hence, we show that the heterogeneous STDP dynamics increase the synaptic
noise due to the heavy tail behavior of the system. This increased synaptic
noise leads to a reduction in the number of spikes of the post-synaptic
neuron. Intuitively, a heterogeneous STDP leads to a non-uniform scaling of
correlated spike-trains leading to de-correlation. Hence, we can say that
heterogeneous STDP models have learned a better-orthogonalized subspace
representation, leading to better encoding of the input space with fewer
spikes. $\blacksquare$

Theorem 2: For a given number of neurons $N_{\mathcal{R}}$, the spike
efficiency of the model
$\mathcal{E}=\frac{\mathcal{C}(N_{\mathcal{R}})}{\tilde{S}}$ for HRSNN
($\mathcal{E}_{R}$) is greater than MRSNN ($\mathcal{E}_{M}$) i.e.,
$\mathcal{E}_{R}\geq\mathcal{E}_{M}$

Short Proof:  (See Suppl. Sec. B for full proof) First, using Lemma 3.2.1, we
show that the number of spikes decreases when we use heterogeneity in the
LTP/LTD Dynamics. Hence, we compare the efficiencies of HRSNN with that of
MRSNN as follows:

$\displaystyle{\color[rgb]{0,0,0}\frac{\mathcal{E}_{R}}{\mathcal{E}_{M}}=\frac{\mathcal{C}_{R}(N_{\mathcal{R}})\times\tilde{S}_{M}}{\tilde{S}_{R}\times\mathcal{C}_{M}(N_{\mathcal{R}})}=\frac{\sum_{\tau=1}^{N_{\mathcal{R}}}\frac{\operatorname{Cov}^{2}\left(x(t-\tau),\mathbf{a}^{R}_{\tau}\mathbf{r}_{R}(t)\right)}{\operatorname{Var}\left(\mathbf{a}^{R}_{\tau}\mathbf{r}_{R}(t)\right)}\times\int\limits_{t_{ref}}^{\infty}t\Phi_{R}dt}{\sum_{\tau=1}^{N_{\mathcal{R}}}\frac{\operatorname{Cov}^{2}\left(x(t-\tau),\mathbf{a}^{M}_{\tau}\mathbf{r}_{M}(t)\right)}{\operatorname{Var}\left(\mathbf{a}^{M}_{\tau}\mathbf{r}_{M}(t)\right)}\times\int\limits_{t_{ref}}^{\infty}t\Phi_{M}dt}}$
(6)

Since $S_{R}\leq S_{M}$ and also,the covariance increases when the neurons
become correlated, and as neuronal correlation decreases, $\mathcal{H}_{X}$
increases (Theorem 1), we see that $\mathcal{E}_{R}/\mathcal{E}_{M}\geq
1\Rightarrow\mathcal{E}_{R}\geq\mathcal{E}_{M}\quad\blacksquare$

Optimal Heterogeneity using Bayesian Optimization for Distributions To get
optimal heterogeneity in the neuron and STDP dynamics, we use a modified
Bayesian Optimization (BO) technique. However, using BO for high-dimensional
problems remains a significant challenge. In our case, optimizing HRSNN model
parameters for 5000 neurons requires the optimization of two parameters per
neuron and four parameters per STDP synapse, where standard BO fails to
converge to an optimal solution. However, the parameters to be optimized are
correlated and can be drawn from a probability distribution as shown by Perez-
Nieves et al. (2021). Thus, we design a modified BO to estimate parameter
distributions instead of individual parameters for the LIF neurons and the
STDP synapses, for which we modify the BO’s surrogate model and acquisition
function. This makes our modified BO highly scalable over all the variables
(dimensions) used. The loss for the surrogate model’s update is calculated
using the Wasserstein distance between the parameter distributions. We use the
modified Matern function on the Wasserstein metric space as a kernel function
for the BO problem. The detailed BO methods are discussed in Suppl. Sec. A. BO
uses a Gaussian process to model the distribution of an objective function and
an acquisition function to decide points to evaluate. For data points $x\in X$
and the corresponding output $y\in Y$, an SNN with network structure
$\mathcal{V}$ and neuron parameters $\mathcal{W}$ acts as a function
$f_{\mathcal{V},\mathcal{W}}(x)$ that maps input data $x$ to $y$. The
optimization problem can be defined as:
$\min_{\mathcal{V},\mathcal{W}}\sum_{x\in X,y\in
Y}\mathcal{L}\left(y,f_{\mathcal{V},\mathcal{W}}(x)\right)$ where
$\mathcal{V}$ is the set of hyperparameters of the neurons in $\mathcal{R}$
and $\mathcal{W}$ is the multi-variate distribution constituting the
distributions of: (i) the membrane time constants $\tau_{m-E},\tau_{m-I}$ of
LIF neurons, (ii) the scaling function constants $(A_{+},A_{-})$ and (iii) the
decay time constants $\tau_{+},\tau_{-}$ for the STDP learning rule in
$\mathcal{S}_{\mathcal{RR}}$.

Figure 2: Block Diagram showing the methodology using HRSNN for prediction

###  4 Experimental Results

Model and Architecture We empirically verify our analytical results using
HRSNN for classification and prediction tasks. Fig. 2 shows the overall
architecture of the prediction model. Using a rate-encoding methodology, the
time-series data is encoded to a series of spike trains. This high-dimensional
spike train acts as the input to HRSNN. The output spike trains from HRSNN act
as the input to a decoder and a readout layer that finally gives the
prediction results. For the classification task, we use a similar method.
However, we do not use the decoding layer for the signal but directly feed the
output spike signals from HRSNN into the fully connected layer. The complete
details of the models used and description of the different modules used in
Fig. 2 is discussed in Suppl. Sec. A.

Figure 3: Figure showing the results of ablation studies of BO for the
following three cases: (a) for prediction problem - the radius indicates the
normalized NRMSE loss with a smaller radius indicating lower (better) NRMSE
(b) for classification problem - the radius indicates the normalized accuracy
with larger radius indicating higher (better) accuracy. The numbers represent
the number of neurons used in each model, and the line joins the three
corresponding models with the same model size. The detailed results are shown
in Suppl. Sec. A

Datasets:  Classification:  We use the Spoken Heidelberg Digits (SHD) spiking
dataset to benchmark the HRSNN model with other standard spiking neural
networks (Cramer et al., 2020).
Prediction:  We use a multiscale Lorenz 96 system (Lorenz, 1996) which is a
set of coupled nonlinear ODEs and an extension of Lorenz’s original model for
multiscale chaotic variability of weather and climate systems which we use as
a testbed for the prediction capabilities of the HRSNN model (Thornes et al.,
2017). Further details on both datasets are provided in Suppl. Sec. A.

Bayesian Optimization Ablation Studies: First, we perform an ablation study of
BO for the following three cases: (i) Using Memory Capacity $\mathcal{C}$ as
the objective function (ii) Using Average Spike Count $\tilde{S}$ as the
objective function and (iii) Using $\mathcal{E}$ as the objective function. We
optimize both LIF neuron parameter distribution and STDP dynamics
distributions for each. We plot $\mathcal{C}$, $\tilde{S}$, the empirical
spike efficiency $\hat{\mathcal{E}}$, and the observed RMSE of the model
obtained from BO with different numbers of neurons. The results for
classification and prediction problems are shown in Fig. 3(a) and (b),
respectively. Ideally, we want to design networks with high $\mathcal{C}$ and
low spike count, i.e., models in the upper right corner of the graph. The
observed results show that BO using $\mathcal{E}$ as the objective gives the
best accuracy with the fewest spikes. Thus, we can say that this model has
learned a better-orthogonalized subspace representation, leading to better
encoding of the input space with fewer spikes. Hence, for the remainder of
this paper, we focus on this BO model, keeping the $\mathcal{E}$ as the
objective function. This Bayesian Optimization process to search for the
optimal hyperparameters of the model is performed before training and
inference using the model and is generally equivalent to the network
architecture search process used in deep learning. Once we have these optimal
hyper-parameters, we freeze these hyperparameters, learn (unsupervised) the
network parameters (i.e., synaptic weights) of the HRSNN while using the
frozen hyperparameters, and generate the final HRSNN model for inference. In
other words, the hyperparameters, like the distribution of membrane time
constants or the distribution of synaptic time constants for STDP, are fixed
during the learning and inference. Further details of the Bayesian
Optimization procedure, including the parameterized variables and the final
searched distributions of the hyperparameters, are shown in Suppl. Sec. A,
where we also discuss the convergence analysis of the three different BOs
discussed above.

Heterogeneity Parameter Importance:  We use SAGE (Shapley Additive Global
importancE) (Covert et al., 2020), a game-theoretic approach to understand
black-box models to calculate the significance of adding heterogeneity to each
parameter for improving $\mathcal{C}$ and $\tilde{S}$. SAGE summarizes the
importance of each feature based on the predictive power it contributes and
considers complex feature interactions using the principles of Shapley value,
with a higher SAGE value signifying a more important feature. We tested the
HRSNN model using SAGE on the Lorenz96 and the SHD datasets. The results are
shown in Fig. 4. We see that $\tau_{m}$ has the greatest SAGE values for
$\mathcal{C}$, signifying that it has the greatest impact on improving
$\mathcal{C}$ when heterogeneity is added. Conversely, we see that
heterogeneous STDP parameters (viz., $\tau_{\pm},\eta_{\pm}$) play a more
critical role in determining the average neuronal spike activation. Hence, we
confirm the notions proved in Sec. 3 that heterogeneity in neuronal dynamics
improves the $\mathcal{C}$ while heterogeneity in STDP dynamics improves the
spike count. Thus, we need to optimize the heterogeneity of both to achieve
maximum $\mathcal{E}$.

Figure 4: Bar chart showing the global importance of different heterogeneous
parameters using HRSNN on the dataset. The experiments were repeated five
times with different parameters from the same distribution (a) Classification
(b) Prediction

Results:  We perform an ablation study to evaluate the performance of the
HRSNN model and compare it to standard BP-based spiking models. We study the
performances of both the SHD dataset for classification and the Lorenz system
for prediction. The results are shown in Table 2. We compare the Normalized
Root Mean Squared Error (NRMSE) loss (prediction), Accuracy (classification),
Average Spike Count $\tilde{S}$ and the application level empirical spiking
efficiency $\hat{\mathcal{E}}$ calculated as
$\displaystyle\frac{1}{\text{NRMSE}\times\tilde{S}}$ (prediction) and
$\displaystyle\frac{\text{Accuracy}}{\tilde{S}}$ (classification).  We perform
the experiments using 5000 neurons in $\mathcal{R}$ on both classification and
prediction datasets. We see that the HRSNN model with heterogeneous LIF and
heterogeneous STDP outperforms other HRSNN and MRSNN models in terms of NRMSE
scores while keeping the $\tilde{S}$ much lower than HRSNN with heterogeneous
LIF and homogeneous STDP. From the experiments, we can conclude that the
heterogeneous LIF neurons have the greatest contribution to improving the
model’s performance. In contrast, heterogeneity in STDP has the most
significant impact on a spike-efficient representation of the data. HRSNN with
heterogeneous LIF and STDP leverages the best of both worlds by achieving the
best RMSE with low spike activations, as seen from Table 2. Further detailed
results on limited training data are added in Suppl. Sec. A. We also compare
the generalizability of the HRSNN vs. MRSNN models, where we empirically show
that the heterogeneity in STDP dynamics helps improve the overall model’s
generalizability. In addition, we discuss how HRSNN reduces the effect of
higher-order correlations, thereby giving rise to a more efficient
representation of the state space.

Table 2: Table showing the comparison of the Accuracy and NRMSE losses for the
SHD Classification and Lorenz System Prediction tasks, respectively. We show
the average spike rate, calculated as the ratio of the moving average of the
number of spikes in a time interval $T$. For this experiment, we choose
$T=4ms$ and a rolling time span of $2ms$, which is repeated until the first
spike appears in the final layer. Following the works of Paul et al. (2022),
we show that the normalized average spike rate is the total number of spikes
generated by all neurons in an RSNN averaged over the time interval $T$. The
results marked with * denotes we implemented the open-source code for the
model and evaluated the given results.

| Method |  | SHD
---
(Classification)
| Chaotic Lorenz System
---
(Prediction)
| Accuracy
---
$(A)$
| Normalized
---
Avg. Firing
Rate ($\frac{\tilde{S}}{T}$)
| Efficiency
---
($\hat{\mathcal{E}}={\color[rgb]{0,0,0}\frac{A}{\tilde{S}}}$)
NRMSE |  | Normalized
---
Avg. Firing
Rate ($\frac{\tilde{S}}{T}$)
| Efficiency
---
$\hat{\mathcal{E}}=\frac{1}{\text{NRMSE}\times\tilde{S}}$
Unsupervised RSNN |  | MRSNN
---
(Homogeneous LIF,
Homogeneous STDP)
73.58 | -0.508 | ${\color[rgb]{0,0,0}18.44\times 10^{-3}}$ | 0.395 | -0.768 | ${\color[rgb]{0,0,0}0.787\times 10^{-3}}$
| HRSNN
---
(Heterogeneous LIF,
Homogeneous STDP)
78.87 | 0.277 | ${\color[rgb]{0,0,0}17.19\times 10^{-3}}$ | 0.203 | -0.143 | ${\color[rgb]{0,0,0}1.302\times 10^{-3}}$
| HRSNN
---
(Homogeneous LIF,
Heterogeneous STDP)
74.03 | -1.292 | ${\color[rgb]{0,0,0}22.47\times 10^{-3}}$ | 0.372 | -1.102 | ${\color[rgb]{0,0,0}0.932\times 10^{-3}}$
| HRSNN
---
(Heterogeneous LIF,
Heterogeneous STDP)
80.49 | -1.154 | ${\color[rgb]{0,0,0}24.35\times 10^{-3}}$ | 0.195 | -1.018 | ${\color[rgb]{0,0,0}1.725\times 10^{-3}}$
RSNN with BP |  | MRSNN-BP
---
(Homogeneous LIF,
BP)
81.42 | 0.554 | ${\color[rgb]{0,0,0}16.9\times 10^{-3}}$ | 0.182 | 0.857 | ${\color[rgb]{0,0,0}1.16\times 10^{-3}}$
| HRSNN-BP
---
(Heterogeneous LIF,
BP)
83.54 | 1.292 | ${\color[rgb]{0,0,0}15.42\times 10^{-3}}$ | 0.178 | 1.233 | ${\color[rgb]{0,0,0}1.09\times 10^{-3}}$
| Adaptive SRNN
---
(Yin et al., 2020)
84.46 | $0.831^{*}$ | $17.21^{*}\times 10^{-3}$ | $0.174^{*}$ | $0.941^{*}$ | $1.19^{*}\times 10^{-3}$

###  5 Conclusion

This paper analytically and empirically proved that heterogeneity in neuronal
(LIF) and synaptic (STDP) dynamics leads to an unsupervised RSNN with more
memory capacity, reduced spiking count, and hence, better spiking efficiency.
We show that HRSNN can achieve similar performance as an MRSNN but with sparse
spiking leading to the improved energy efficiency of the network. In
conclusion, this work establishes important mathematical properties of an RSNN
for neuromorphic machine learning applications like time series classification
and prediction. However, it is interesting to note that the mathematical
results from this paper also conform to the recent neurobiological research
that suggests that the brain has large variability between the types of
neurons and learning methods. For example, intrinsic biophysical properties of
neurons, like densities and properties of ionic channels, vary significantly
between neurons where the variance in synaptic learning rules invokes reliable
and efficient signal processing in several animals (Marder & Taylor, 2011;
Douglass et al., 1993). Experiments in different brain regions and diverse
neuronal types have revealed a wide range of STDP forms with varying neuronal
dynamics that vary in plasticity direction, temporal dependence, and the
involvement of signaling pathways (Sjostrom et al., 2008; Korte & Schmitz,
2016). Thus, heterogeneity is essential in encoding and decoding stimuli in
biological systems. In conclusion, this work establishes connections between
the mathematical properties of an RSNN for neuromorphic machine learning
applications like time series classification and prediction with
neurobiological observations.  There are some key limitations to the analyses
in this paper. First, the properties discussed are derived independently. An
important extension will be to consider all factors simultaneously. Second, we
assumed an idealized spiking network where the memory capacity is used to
measure its performance, and the spike count measures the energy. Also, we
mainly focused on the properties of RSNN trained using STDP. An interesting
connection between synchronization and heterogeneous STDP remains a topic that
needs to be studied further - whether we can optimally engineer the
synchronization properties to improve the model’s performance. Finally, the
empirical evaluations were presented for the prediction task on a single
dataset. More experimental evaluations, including other tasks and datasets,
will strengthen the empirical validations.

### Acknowledgement

This work is supported by the Army Research Office and was accomplished under
Grant Number W911NF-19-1-0447. The views and conclusions contained in this
document are those of the authors and should not be interpreted as
representing the official policies, either expressed or implied, of the Army
Research Office or the U.S. Government.

### References

  * Abbott & Dayan (1999) Larry F Abbott and Peter Dayan.  The effect of correlated variability on the accuracy of a population code.  _Neural computation_ , 11(1):91–101, 1999.
  * Aceituno et al. (2020) Pau Vilimelis Aceituno, Gang Yan, and Yang-Yu Liu.  Tailoring echo state networks for optimal learning.  _iscience_ , 23(9):101440, 2020.
  * Blumenthal & Getoor (1960) Robert M Blumenthal and Ronald K Getoor.  Some theorems on stable processes.  _Transactions of the American Mathematical Society_ , 95(2):263–273, 1960.
  * Boerlin & Denève (2011) Martin Boerlin and Sophie Denève.  Spike-based population coding and working memory.  _PLoS computational biology_ , 7(2):e1001080, 2011\.
  * Boerlin et al. (2013) Martin Boerlin, Christian K Machens, and Sophie Denève.  Predictive coding of dynamical variables in balanced spiking networks.  _PLoS computational biology_ , 9(11):e1003258, 2013.
  * Bourdoukan et al. (2012) Ralph Bourdoukan, David Barrett, Sophie Deneve, and Christian K Machens.  Learning optimal spike-based representations.  _Advances in neural information processing systems_ , 25, 2012.
  * Brunel (2000) Nicolas Brunel.  Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons.  _Journal of computational neuroscience_ , 8(3):183–208, 2000.
  * Buzsáki & Mizuseki (2014) György Buzsáki and Kenji Mizuseki.  The log-dynamic brain: how skewed distributions affect network operations.  _Nature Reviews Neuroscience_ , 15(4):264–278, 2014.
  * Chakraborty & Mukhopadhyay (2021) Biswadeep Chakraborty and Saibal Mukhopadhyay.  Characterization of generalizability of spike time dependent plasticity trained spiking neural networks.  _arXiv preprint arXiv:2105.14677_ , 2021.
  * Chakraborty & Mukhopadhyay (2023a) Biswadeep Chakraborty and Saibal Mukhopadhyay.  Brain-inspired spiking neural network for online unsupervised time series prediction.  _arXiv preprint arXiv:2304.04697_ , 2023a.
  * Chakraborty & Mukhopadhyay (2023b) Biswadeep Chakraborty and Saibal Mukhopadhyay.  Heterogeneous recurrent spiking neural network for spatio-temporal classification.  _Frontiers in Neuroscience_ , 17, 2023b.  ISSN 1662-453X.  doi: 10.3389/fnins.2023.994517.  URL https://www.frontiersin.org/articles/10.3389/fnins.2023.994517.
  * Chakraborty et al. (2021) Biswadeep Chakraborty, Xueyuan She, and Saibal Mukhopadhyay.  A fully spiking hybrid neural network for energy-efficient object detection.  _arXiv preprint arXiv:2104.10719_ , 2021.
  * Chakraborty et al. (2023) Biswadeep Chakraborty, Uday Kamal, Xueyuan She, Saurabh Dash, and Saibal Mukhopadhyay.  Brain-inspired spatiotemporal processing algorithms for efficient event-based perception.  In _2023 Design, Automation & Test in Europe Conference & Exhibition (DATE)_, pp. 1–6. IEEE, 2023.
  * Chattopadhyay et al. (2020) Ashesh Chattopadhyay, Pedram Hassanzadeh, and Devika Subramanian.  Data-driven predictions of a multiscale lorenz 96 chaotic system using machine-learning methods: reservoir computing, artificial neural network, and long short-term memory network.  _Nonlinear Processes in Geophysics_ , 27(3):373–389, 2020.
  * Chevallier et al. (2015) Julien Chevallier, María José Cáceres, Marie Doumic, and Patricia Reynaud-Bouret.  Microscopic approach of a time elapsed neural model.  _Mathematical Models and Methods in Applied Sciences_ , 25(14):2669–2719, 2015.
  * Chornoboy et al. (1988) ES Chornoboy, LP Schramm, and AF Karr.  Maximum likelihood identification of neural point process systems.  _Biological cybernetics_ , 59(4-5):265–275, 1988\.
  * Cossell et al. (2015) Lee Cossell, Maria Florencia Iacaruso, Dylan R Muir, Rachael Houlton, Elie N Sader, Ho Ko, Sonja B Hofer, and Thomas D Mrsic-Flogel.  Functional organization of excitatory synaptic strength in primary visual cortex.  _Nature_ , 518(7539):399–403, 2015.
  * Covert et al. (2020) Ian Covert, Scott M Lundberg, and Su-In Lee.  Understanding global feature contributions with additive importance measures.  _Advances in Neural Information Processing Systems_ , 33:17212–17223, 2020.
  * Cramer et al. (2020) Benjamin Cramer, Yannik Stradmann, Johannes Schemmel, and Friedemann Zenke.  The heidelberg spiking data sets for the systematic evaluation of spiking neural networks.  _IEEE Transactions on Neural Networks and Learning Systems_ , 2020\.
  * Denève & Machens (2016) Sophie Denève and Christian K Machens.  Efficient codes and balanced networks.  _Nature neuroscience_ , 19(3):375–382, 2016.
  * Destexhe & Marder (2004) Alain Destexhe and Eve Marder.  Plasticity in single neuron and circuit computations.  _Nature_ , 431(7010):789–795, 2004.
  * Douglass et al. (1993) John K Douglass, Lon Wilkens, Eleni Pantazelou, and Frank Moss.  Noise enhancement of information transfer in crayfish mechanoreceptors by stochastic resonance.  _Nature_ , 365(6444):337–340, 1993.
  * Duval et al. (2022) Céline Duval, Eric Luçon, and Christophe Pouzat.  Interacting hawkes processes with multiplicative inhibition.  _Stochastic Processes and their Applications_ , 148:180–226, 2022.
  * Farkaš & Gergel’ (2017) Igor Farkaš and Peter Gergel’.  Maximizing memory capacity of echo state networks with orthogonalized reservoirs.  In _2017 International Joint Conference on Neural Networks (IJCNN)_ , pp. 2437–2442. IEEE, 2017.
  * Farkaš et al. (2016) Igor Farkaš, Radomír Bosák, and Peter Gergel’.  Computational analysis of memory capacity in echo state networks.  _Neural Networks_ , 83:109–120, 2016.
  * Feydy et al. (2019) Jean Feydy, Thibault Séjourné, François-Xavier Vialard, Shun-ichi Amari, Alain Trouvé, and Gabriel Peyré.  Interpolating between optimal transport and mmd using sinkhorn divergences.  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , pp. 2681–2690. PMLR, 2019.
  * Frazier (2018) Peter I Frazier.  A tutorial on bayesian optimization.  _arXiv preprint arXiv:1807.02811_ , 2018.
  * Galves & Löcherbach (2016) Antonio Galves and Eva Löcherbach.  Modeling networks of spiking neurons as interacting processes with memory of variable length.  _Journal de la Société Française de Statistique_ , 157(1):17–32, 2016.
  * Gelenbe et al. (1999) Erol Gelenbe, Zhi-Hong Mao, and Yan-Da Li.  Function approximation with spiked random networks.  _IEEE Transactions on Neural Networks_ , 10(1):3–9, 1999.
  * Gerhard et al. (2017) Felipe Gerhard, Moritz Deger, and Wilson Truccolo.  On the stability and dynamics of stochastic spiking neuron models: Nonlinear hawkes process and point process glms.  _PLoS computational biology_ , 13(2):e1005390, 2017.
  * Gerstner & Kistler (2002) Wulfram Gerstner and Werner M Kistler.  Mathematical formulations of hebbian learning.  _Biological cybernetics_ , 87(5):404–415, 2002\.
  * Gilra & Gerstner (2017) Aditya Gilra and Wulfram Gerstner.  Predicting non-linear dynamics by stable local learning in a recurrent spiking neural network.  _Elife_ , 6:e28295, 2017.
  * Goldmann et al. (2020) Mirko Goldmann, Felix Köster, Kathy Lüdge, and Serhiy Yanchuk.  Deep time-delay reservoir computing: Dynamics and memory capacity.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , 30(9):093124, 2020.
  * Gouwens et al. (2019) Nathan W Gouwens, Staci A Sorensen, Jim Berg, Changkyu Lee, Tim Jarsky, Jonathan Ting, Susan M Sunkin, David Feng, Costas A Anastassiou, Eliza Barkan, et al.  Classification of electrophysiological and morphological neuron types in the mouse visual cortex.  _Nature neuroscience_ , 22(7):1182–1195, 2019\.
  * Hansel et al. (1995) David Hansel, Germán Mato, and Claude Meunier.  Synchrony in excitatory neural networks.  _Neural computation_ , 7(2):307–337, 1995.
  * Hansen et al. (2015) Niels Richard Hansen, Patricia Reynaud-Bouret, and Vincent Rivoirard.  Lasso and probabilistic inequalities for multivariate point processes.  2015\.
  * Hayashi & Igarashi (2009) Hatsuo Hayashi and Jun Igarashi.  Ltd windows of the stdp learning rule and synaptic connections having a large transmission delay enable robust sequence learning amid background noise.  _Cognitive neurodynamics_ , 3(2):119–130, 2009\.
  * Iannella & Back (2001) Nicolangelo Iannella and Andrew D. Back.  A spiking neural network architecture for nonlinear function approximation.  _Neural Networks_ , 14(6):933–939, 2001.  ISSN 0893-6080.  doi: https://doi.org/10.1016/S0893-6080(01)00080-6.  URL https://www.sciencedirect.com/science/article/pii/S0893608001000806.
  * Jaeger (2001) Herbert Jaeger.  The “echo state” approach to analysing and training recurrent neural networks-with an erratum note.  _Bonn, Germany: German National Research Center for Information Technology GMD Technical Report_ , 148(34):13, 2001.
  * Jaeger (2002) Herbert Jaeger.  Short term memory in echo state networks. gmd-report 152.  In _GMD-German National Research Institute for Computer Science (2002), http://www. faculty. jacobs-university. de/hjaeger/pubs/STMEchoStatesTechRep. pdf_. Citeseer, 2002.
  * Jaeger et al. (2001) Herbert Jaeger et al.  _Short term memory in echo state networks_ , volume 5.  GMD-Forschungszentrum Informationstechnik Bremen, Germany, 2001.
  * Khoei et al. (2020) Mina A Khoei, Amirreza Yousefzadeh, Arash Pourtaherian, Orlando Moreira, and Jonathan Tapson.  Sparnet: Sparse asynchronous neural network execution for energy efficient inference.  In _2020 2nd IEEE International Conference on Artificial Intelligence Circuits and Systems (AICAS)_ , pp. 256–260. IEEE, 2020.
  * Kim et al. (2022) Daehyun Kim, Biswadeep Chakraborty, Xueyuan She, Edward Lee, Beomseok Kang, and Saibal Mukhopadhyay.  Moneta: A processing-in-memory-based hardware platform for the hybrid convolutional spiking neural network with online learning.  _Frontiers in Neuroscience_ , 16, 2022.
  * Kim et al. (2020) Seijoon Kim, Seongsik Park, Byunggook Na, and Sungroh Yoon.  Spiking-yolo: Spiking neural network for energy-efficient object detection.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 34, pp. 11270–11277, 2020.
  * Kingma & Ba (2014) Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Korte & Schmitz (2016) Martin Korte and Dietmar Schmitz.  Cellular and system biology of memory: timing, molecules, and beyond.  _Physiological reviews_ , 96(2):647–693, 2016\.
  * Kuśmierz et al. (2020) Łukasz Kuśmierz, Shun Ogawa, and Taro Toyoizumi.  Edge of chaos and avalanches in neural networks with heavy-tailed synaptic weight distribution.  _Physical Review Letters_ , 125(2):028101, 2020\.
  * Lee et al. (2017) Jaehoon Lee, Yasaman Bahri, Roman Novak, Samuel S Schoenholz, Jeffrey Pennington, and Jascha Sohl-Dickstein.  Deep neural networks as gaussian processes.  _arXiv preprint arXiv:1711.00165_ , 2017.
  * Li et al. (2017) Hongmin Li, Hanchao Liu, Xiangyang Ji, Guoqi Li, and Luping Shi.  Cifar10-dvs: an event-stream dataset for object classification.  _Frontiers in neuroscience_ , 11:309, 2017.
  * Löcherbach (2017) Eva Löcherbach.  Spiking neurons: interacting hawkes processes, mean field limits and oscillations.  _ESAIM: Proceedings and Surveys_ , 60:90–103, 2017.
  * Lorenz (1963) Edward N Lorenz.  Deterministic nonperiodic flow.  _Journal of atmospheric sciences_ , 20(2):130–141, 1963.
  * Lorenz (1996) Edward N Lorenz.  Predictability: A problem partly solved.  In _Proc. Seminar on predictability_ , volume 1, 1996.
  * Ly (2015) Cheng Ly.  Firing rate dynamics in recurrent spiking neural networks with intrinsic and network heterogeneity.  _Journal of computational neuroscience_ , 39(3):311–327, 2015.
  * Mantegna & Stanley (1995) Rosario N Mantegna and H Eugene Stanley.  Scaling behaviour in the dynamics of an economic index.  _Nature_ , 376(6535):46–49, 1995.
  * Marder & Taylor (2011) Eve Marder and Adam L Taylor.  Multiple models to capture the variability in biological neurons and networks.  _Nature neuroscience_ , 14(2):133–138, 2011.
  * Mascart (2021) Cyrille Mascart.  _Efficient simulation of point processes with applications to neurosciences_.  PhD thesis, Université Côte d’Azur, Nice, France, 2021.
  * McKenzie et al. (2021) Sam McKenzie, Roman Huszár, Daniel F English, Kanghwan Kim, Fletcher Christensen, Euisik Yoon, and György Buzsáki.  Preexisting hippocampal network dynamics constrain optogenetically induced place fields.  _Neuron_ , 109(6):1040–1054, 2021.
  * Montani et al. (2009) Fernando Montani, Robin AA Ince, Riccardo Senatore, Ehsan Arabzadeh, Mathew E Diamond, and Stefano Panzeri.  The impact of high-order interactions on the rate of synchronous discharge and information transmission in somatosensory cortex.  _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 367(1901):3297–3310, 2009.
  * Neftci et al. (2019) Emre O Neftci, Hesham Mostafa, and Friedemann Zenke.  Surrogate gradient learning in spiking neural networks.  _IEEE Signal Processing Magazine_ , 36:61–63, 2019.
  * Oja (1982) Erkki Oja.  Simplified neuron model as a principal component analyzer.  _Journal of mathematical biology_ , 15(3):267–273, 1982.
  * Oja (1989) Erkki Oja.  Neural networks, principal components, and subspaces.  _International journal of neural systems_ , 1(01):61–68, 1989.
  * Paul et al. (2022) Ankita Paul, Stefan Wagner, and Anup Das.  Learning in feedback-driven recurrent spiking neural networks using full-force training.  In _2022 International Joint Conference on Neural Networks (IJCNN)_ , pp. 1–10. IEEE, 2022.
  * Perez-Nieves et al. (2021) Nicolas Perez-Nieves, Vincent CH Leung, Pier Luigi Dragotti, and Dan FM Goodman.  Neural heterogeneity promotes robust learning.  _Nature communications_ , 12(1):1–9, 2021.
  * Petro et al. (2019) Balint Petro, Nikola Kasabov, and Rita M Kiss.  Selection and optimization of temporal spike encoding methods for spiking neural networks.  _IEEE transactions on neural networks and learning systems_ , 31(2):358–370, 2019.
  * Pfaffelhuber et al. (2022) Peter Pfaffelhuber, Stefan Rotter, and Jakob Stiefel.  Mean-field limits for non-linear hawkes processes with excitation and inhibition.  _Stochastic Processes and their Applications_ , 153:57–78, 2022.
  * Ponulak & Kasinski (2011) Filip Ponulak and Andrzej Kasinski.  Introduction to spiking neural networks: Information processing, learning and applications.  _Acta neurobiologiae experimentalis_ , 71(4):409–433, 2011.
  * Pool & Mato (2011) R Rossi Pool and Germán Mato.  Spike-timing-dependent plasticity and reliability optimization: the role of neuron dynamics.  _Neural computation_ , 23(7):1768–1789, 2011\.
  * Prescott et al. (2008) Steven A Prescott, Yves De Koninck, and Terrence J Sejnowski.  Biophysical basis for three distinct dynamical mechanisms of action potential initiation.  _PLoS computational biology_ , 4(10):e1000198, 2008.
  * Pyle & Rosenbaum (2017) Ryan Pyle and Robert Rosenbaum.  Spatiotemporal dynamics and reliable computations in recurrent spiking neural networks.  _Physical review letters_ , 118(1):018103, 2017\.
  * Rathi et al. (2021) Nitin Rathi, Amogh Agrawal, Chankyu Lee, Adarsh Kumar Kosta, and Kaushik Roy.  Exploring spike-based learning for neuromorphic computing: Prospects and perspectives.  In _2021 Design, Automation & Test in Europe Conference & Exhibition (DATE)_, pp. 902–907. IEEE, 2021.
  * Reynaud-Bouret et al. (2014) Patricia Reynaud-Bouret, Vincent Rivoirard, Franck Grammont, and Christine Tuleau-Malot.  Goodness-of-fit tests and nonparametric adaptive estimation for spike train analysis.  _The Journal of Mathematical Neuroscience_ , 4:1–41, 2014\.
  * She et al. (2022) Xueyuan She, Saurabh Dash, and Saibal Mukhopadhyay.  Sequence approximation using feedforward spiking neural network for spatiotemporal learning: Theory and optimization methods.  In _Proceedings of 10th International Conference on Learning Representations, ICLR (forthcoming)_ , 2022.  URL https://openreview.net/forum?id=bp-LJ4y_XC.
  * Shlesinger et al. (1987) Michael F Shlesinger, BJ West, and Joseph Klafter.  Lévy dynamics of enhanced diffusion: Application to turbulence.  _Physical Review Letters_ , 58(11):1100, 1987\.
  * Simsekli et al. (2019) Umut Simsekli, Levent Sagun, and Mert Gurbuzbalaban.  A tail-index analysis of stochastic gradient noise in deep neural networks.  _arXiv preprint arXiv:1901.06053_ , 2019.
  * Simsekli et al. (2020) Umut Simsekli, Ozan Sener, George Deligiannidis, and Murat A Erdogdu.  Hausdorff dimension, heavy tails, and generalization in neural networks.  _Advances in Neural Information Processing Systems_ , 33, 2020.
  * Sjostrom et al. (2008) P Jesper Sjostrom, Ede A Rancz, Arnd Roth, and Michael Hausser.  Dendritic excitability and synaptic plasticity.  _Physiological reviews_ , 88(2):769–840, 2008\.
  * Sorbaro et al. (2020) Martino Sorbaro, Qian Liu, Massimo Bortone, and Sadique Sheik.  Optimizing the energy consumption of spiking neural networks for neuromorphic applications.  _Frontiers in neuroscience_ , 14:662, 2020.
  * Soures & Kudithipudi (2019) Nicholas Soures and Dhireesha Kudithipudi.  Spiking reservoir networks: Brain-inspired recurrent algorithms that use random, fixed synaptic strengths.  _IEEE Signal Processing Magazine_ , 36(6):78–87, 2019.
  * Srinivasan & Roy (2019) Gopalakrishnan Srinivasan and Kaushik Roy.  Restocnet: Residual stochastic binary convolutional spiking neural network for memory-efficient neuromorphic computing.  _Frontiers in neuroscience_ , 13:189, 2019.
  * Staude et al. (2010) Benjamin Staude, Stefan Rotter, and Sonja Grün.  Cubic: cumulant based inference of higher-order correlations in massively parallel spike trains.  _Journal of computational neuroscience_ , 29(1):327–350, 2010.
  * Thornes et al. (2017) Tobias Thornes, Peter Düben, and Tim Palmer.  On the use of scale-dependent precision in earth system modelling.  _Quarterly Journal of the Royal Meteorological Society_ , 143(703):897–908, 2017.
  * Wu et al. (2019) Yujie Wu, Lei Deng, Guoqi Li, Jun Zhu, Yuan Xie, and Luping Shi.  Direct training for spiking neural networks: Faster, larger, better.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 33, pp. 1311–1318, 2019.
  * Yin et al. (2020) Bojian Yin, Federico Corradi, and Sander M Bohté.  Effective and efficient computation with multiple-timescale spiking recurrent neural networks.  In _International Conference on Neuromorphic Systems 2020_ , pp. 1–8, 2020.
  * Zhang & Li (2020) Wenrui Zhang and Peng Li.  Temporal spike sequence learning via backpropagation for deep spiking neural networks.  _Advances in Neural Information Processing Systems_ , 33:12022–12033, 2020.
  * Zhou et al. (2020) Yan Zhou, Yaochu Jin, and Jinliang Ding.  Surrogate-assisted evolutionary search of spiking neural architectures in liquid state machines.  _Neurocomputing_ , 406:12–23, 2020.

##  Supplementary Section

###  Appendix A Supplementary Section A

####  A.1 Experimental Setup

#####  A.1.1 HRSNN and MRSNN

Table 3: Table showing the description of the models described in the paper

| LIF Neurons | STDP parameters
---|---|---
HRSNN | Heterogeneous | Heterogeneous
MRSNN | Homogeneous | Homogeneous

The models used in this paper were the Heterogeneous Recurrent SNN (HRSNN) and
the Homogeneous Recurrent SNN (MRSNN). Both models use STDP as the learning
method. For MRSNN, we use STDP with uniform parameters for all the synapses.
However, for HRSNN, we use a distribution for each parameter to get a rich
class of diverse LTP/LTD dynamics. But, at the core, all the training is done
using STDP.

#####  A.1.2 LIF Neuron Numerical Implementation

To implement the LIF model, we discretize time into multiples of a small-time
step $\Delta t$ so that spikes can only happen at multiples of $\Delta t$.
(Cramer et al., 2020; Perez-Nieves et al., 2021) Thus, we can approximately
solve Eq. 1 as

$v_{i}(t+\Delta
t)=v_{i}[t+1]=\beta\left(v_{i}[t]-v_{0}\right)+v_{0}+(1-\beta)I_{i}[t]-\left(v_{th}-v_{r}\right)S_{i}[t]\quad\text{s.t.
}\beta=\exp^{(-\Delta t/\tau_{\mathrm{m}})}$ (7)

It is to be noted here that we use this approximation for numerically solving
the LIF neurons. Hence, although we use continuous notations for the remainder
of the paper, it is to be noted that we use the discrete form discussed here
for numerical solutions.

#####  A.1.3 Spike Coding

Encoding:  For the RSNN to process our time series, the signal must be
represented as spikes. We use a temporal encoding technique for representing
signals in this paper. The spikes are only generated whenever the signal
changes in value. The implementation of the temporal encoding used in this
research is based on the Step-Forward (SF) algorithm (Petro et al., 2019). The
percentage of neurons to input the spikes to ($\alpha$) is also chosen to
provide good recurrent layer dynamics.

Decoding:  To represent the recurrent state, we use an exponentially
decreasing rate decoding strategy by taking the sum of all the spikes $s$ over
the last $\tau$ timesteps into account as follows:

$x_{i}^{X}(t)=\sum_{n=0}^{\tau}\gamma^{n}s_{i}(t-n)\quad\forall i\in E$

where X denotes the model representation. The parameters $\tau$ and $\gamma$
are balanced to optimize the memory size of the stored data (e.g., $\tau\leq
50$ ) and its containment of information, which includes adjusting $\tau$ to
the pace at which the temporal data is presented and processed. The state of
the recurrent layer will be only based on the output of excitatory neurons.
Thus, it is crucial for the discount $\gamma$ not to be too small, as it
possibly flattens older values in the window to 0, making part of the sliding
window unusable. Recent spikes hardly affect the recurrent layer state when
setting $\gamma$ too high in combination with a large window size. This causes
the decoder to react too late to recent information provided by the recurrent
layer and complicates the learning process of the readout layer.

#####  A.1.4 Readout

After the initialization of the recurrent layer, the readout is the only
component of the LSM with trainable parameters. It consists of a single fully
connected layer for regression or classification. The readout does not have to
be any deeper, as the output of the recurrent layer is already a high-
dimensional representation of the processed input. $x$ and $y$ present the
continuous signals of the time series $T$ and model representation $X$.

$x^{T}(t+k)=y^{T}(t)\approx\hat{y}^{X}(t)=f_{\theta}\left(x^{X}(t)\right)$

The mean squared error (MSE) is used as the loss function to train the
readout, and the network was trained using the stochastic optimizer Adam
(Kingma & Ba, 2014).

$\mathcal{L}\left(y^{T},\hat{y}^{X}\right)=\frac{1}{n}\sum_{i=0}^{n}\left(y_{i}^{T}-\hat{y}_{i}^{X}\right)^{2}$

#####  A.1.5 Datasets

Lorenz96: (Lorenz, 1996)Our objective is more clearly demonstrated using the
canonical chaotic system we will use as a test bed for the prediction
capabilities of the HRSNN model. We use a multiscale Lorenz 96 system which is
a set of coupled nonlinear ODEs and an extension of Lorenz’s original model
(Thornes et al., 2017), (Chattopadhyay et al., 2020).

$\displaystyle\frac{\mathrm{d}X_{k}}{\mathrm{~{}d}t}=\quad$ $\displaystyle
X_{k-1}\left(X_{k+1}-X_{k-2}\right)+F-\frac{hc}{b}\Sigma_{j}Y_{j,k}$
$\displaystyle\frac{\mathrm{d}Y_{j,k}}{\mathrm{~{}d}t}=$ $\displaystyle-
cbY_{j+1,k}\left(Y_{j+2,k}-Y_{j-1,k}\right)-cY_{j,k}+\frac{hc}{b}X_{k}-\frac{he}{d}\Sigma_{i}Z_{i,j,k}$
$\displaystyle\frac{\mathrm{d}Z_{i,j,k}}{\mathrm{~{}d}t}=$ $\displaystyle
edZ_{i-1,j,k}\left(Z_{i+1,j,k}-Z_{i-2,j,k}\right)-geZ_{i,j,k}+\frac{he}{d}Y_{j,k}$
(8)

This set of coupled nonlinear ordinary differential equations (ODEs) is a
three-tier extension of Lorenz’s original model (Lorenz, 1963) and has been
proposed by Thornes et al.Thornes et al. (2017) as a fitting prototype for
multiscale chaotic variability of the weather and climate system and a useful
test bed for novel methods. In these equations, $F=20$ is a large-scale
forcing that makes the system highly chaotic, and $b=c=e=d=g=10$ and $h=1$ are
tuned to produce appropriate spatiotemporal variability. For this paper, we
focus on predicting $Y$ axes, which have relatively moderate amplitudes
compared to $X,Z$ and demonstrate high-frequency variability and
intermittency, which makes the prediction problem difficult. It is to be noted
here that the Lorenz 96 is a complex, difficult dataset for climate modeling.
A snippet of the time series is shown in Fig. 5.

Figure 5:  Figure showing a snippet for the Y dimension of the Lorenz96 time
series used for the prediction problem

SHD dataset: We use the Spoken Heidelberg Digits spiking dataset to benchmark
the HRSNN model with other standard spiking neural networks (Cramer et al.,
2020). It was created based on the Heidelberg Digits (HD) audio dataset which
comprises 20 classes of spoken digits from zero to nine in English and German,
spoken by 12 individuals. For training and evaluation, the dataset (10420
samples) is split into a training set (8156 samples) and test set (2264
samples). To apply our RSNNs, we converted all audio samples into 250- by-700
binary matrices. For this, all samples fit within a 1 the second window,
shorter samples were padded with zeros, and longer samples were cut by
removing the tail. Spikes were then binned in time bins, both of sizes 10ms
and 4ms; for the RSNNs, the presence or non-presence of any spikes in the time
bin is noted as a single binary event.

#####  A.1.6 Hyperparameters

The hyperparameters used in this paper are summarized in Table 4

Table 4: Table showing the hyperparameters used in the experiments and their
values

Parameter | Value | Description
---|---|---
|E|/|N| | 80% | Excitatory/inhibitory ratio
$\lambda$ | 2 | Leak Exponent
$\tau$ | 50 | Sliding window size
$\gamma^{\tau-1}$ | 0.02 | Sliding Window Leak
SHD Parameters | Lorenz-63 Parameters
Parameter | Value | Description | Parameter | Value | Description
$\tau$ | 17 | Time delay | $\rho$ | 28.0 |  $\rho$-parameter
$a$ | 0.2 | a parameter | $\sigma$ | 10.0 |  $\sigma$-parameter
$b$ | 0.1 | b parameter | $\beta$ | 8/3 |  $\beta$-parameter
n | 10 | n parameter | $x_{0}$ | [1.0,1.0,1.0] | Initial Condition
$x_{0}$ | 1.2 | Initial Condition | h | 0.03 |  | Time delta between
---
two discrete timesteps
h | 1.0 |  | Time delta between
---
two discrete timesteps
|  |

####  A.2 Maximum entropy distribution

This subsection proves that the maximum entropy distribution with a fixed
covariance matrix is Gaussian.

Lemma: Let $q(\mathbf{r})$ be any density satisfying $\int
q(\mathbf{r})x_{i}x_{j}d\mathbf{r}=\Sigma_{ij}$. Let
$p=\mathcal{N}(\mathbf{0},\bm{\Sigma})$. Then $h(q)\leq h(p)$
Proof.

$\displaystyle 0$ $\displaystyle\leq\mathbb{K}\mathbb{L}(q\|p)=\int
q(\mathbf{r})\log\frac{q(\mathbf{r})}{p(\mathbf{r})}d\mathbf{r}$
$\displaystyle=-h(q)-\int q(\mathbf{r})\log p(\mathbf{r})d\mathbf{r}$
$\displaystyle=\quad-h(q)-\int p(\mathbf{r})\log p(\mathbf{r})d\mathbf{r}$
$\displaystyle=-h(q)+h(p)$

since $q$ and $p$ yield the same moments for the quadratic form encoded by
$\log p(\mathbf{r})$.

####  A.3 Optimal Hyperparameter Selection using Bayesian Optimization

Most recent research in Bayesian Optimization (BO) applications is limited to
low-dimensional problems, as BO fails catastrophically when generalizing to
high-dimensional problems (Frazier, 2018). However, in this paper, we aim to
use BO to optimize the neuronal and synaptic parameters of a heterogeneous
RSNN model. This BO problem thus entails a huge number of hyperparameters to
be optimized; hence, using standard BO algorithms remains a significant
challenge. Hence, to overcome this issue, we used a novel BO algorithm based
on the assumption that our hyperparameters to be optimized are not completely
random and uncorrelated but can be thought of as being drawn from a
probability distribution as shown by Perez et al.Perez-Nieves et al. (2021).
Thus, instead of searching for the individual parameters themselves, we use a
modified BO to estimate parameter distributions for the LIF neurons and the
STDP dynamics. After learning the optimal distributions, we simply sample from
the distribution to get the distribution of hyperparameters used in the model.
To learn the probability distribution of the data, we modify BO’s surrogate
model and acquisition function to treat the parameter distributions instead of
individual variables. This makes our modified BO highly scalable over all the
variables (dimensions) used. The loss for the surrogate model’s update is
calculated using the Wasserstein distance between the parameter distributions.

BO uses a Gaussian process to model the distribution of an objective function
and an acquisition function to decide on points to evaluate. For data points
in a target dataset $x\in X$ and the corresponding label $y\in Y$, an SNN with
network structure $\mathcal{V}$ and neuron parameters $\mathcal{W}$ acts as a
function $f_{\mathcal{V},\mathcal{W}}(x)$ that maps input data $x$ to
predicted label $\tilde{y}$. The optimization problem in this work is defined
as

$\min_{\mathcal{V},\mathcal{W}}\sum_{x\in X,y\in
Y}\mathcal{L}\left(y,f_{\mathcal{V},\mathcal{W}}(x)\right)$ (9)

where $\mathcal{V}$ is the set of hyperparameters of the neurons in
$\mathcal{R}$ (Details of hyperparameters given in the Supplementary) and
$\mathcal{W}$ is the multi-variate distribution constituting the distributions
of (i) the membrane time constants $\tau_{m-E},\tau_{m-I}$ of the LIF neurons,
(ii) the scaling function constants $(A_{+},A_{-})$ and (iii) the decay time
constants $\tau_{+},\tau_{-}$ for the STDP learning rule in
$\mathcal{S}_{\mathcal{RR}}$.

Again, BO needs a prior distribution of the objective function $f(\vec{x})$ on
the given data
$\mathcal{D}_{1:k}=\left\\{\vec{x}_{1:k},f\left(\vec{x}_{1:k}\right)\right\\}.$
In the Gaussian Process (GP)-based BO, we assume that the prior distribution
of $f\left(\vec{x}_{1:k}\right)$ follows the multivariate Gaussian
distribution, which follows a GP with mean $\vec{\mu}_{\mathcal{D}_{1:k}}$ and
covariance $\vec{\Sigma}_{\mathcal{D}_{1:k}}$. Thus, we estimate
$\vec{\Sigma}_{\mathcal{D}_{1:k}}$ using the modified Matern kernel function.
We use the loss function as $d(x,x^{\prime})$, which is the Wasserstein
distance between the multivariate distributions of the different parameters.
That is, given two distributions of hyperparameters $x_{1},x_{2}$, the
distance between these two distributions (given as $d(x_{1},x_{2})$ is used as
the loss function in the Matern kernel for the modified BO. We want to learn
the optimal distribution of hyperparameters $x^{\prime}$, which maximizes the
performance. It is to be noted here that for higher-dimensional metric spaces,
we use the Sinkhorn distance as a regularized version of the Wasserstein
distance to approximate the Wasserstein distance (Feydy et al., 2019).

$\mathcal{D}_{1:k}$ are the points evaluated by the objective function. The GP
will estimate the mean $\vec{\mu}_{\mathcal{D}_{k:n}}$ and variance
$\vec{\sigma}_{\mathcal{D}_{k:n}}$ for the rest unevaluated data
$\mathcal{D}_{k:n}$. The acquisition function used in this work is the
expected improvement (EI) of the prediction fitness as:

$EI\left(\vec{x}_{k:n}\right)=\left(\vec{\mu}_{\mathcal{D}_{k:n}}-f\left(x_{\text{best
}}\right)\right)\Phi(\vec{Z})+\vec{\sigma}_{\mathcal{D}_{k:n}}\phi(\vec{Z})$
(10)

where $\Phi(\cdot)$ and $\phi(\cdot)$ denote the probability distribution
function and the cumulative distribution function of the prior distributions,
respectively. $f\left(x_{\text{best }}\right)=\max
f\left(\vec{x}_{1:k}\right)$ is the maximum value that has been evaluated by
the original function $f$ in all evaluated data $\mathcal{D}_{1:k}$ and
$\vec{Z}=\frac{\vec{\mu}_{\mathcal{D}_{k:n}}-f\left(x_{\text{best
}}\right)}{\vec{\sigma}_{\mathcal{D}_{k:n}}}$. BO will choose the data
$x_{j}=$
$\operatorname{argmax}\left\\{EI\left(\vec{x}_{k:n}\right);x_{j}\subseteq\vec{x}_{k:n}\right\\}$
as the next point to be evaluated using the original objective function.

#####  A.3.1 Optimized Hyperparameters

The list of the hyperparameters optimized using the Bayesian Optimization
technique is shown in Table 5. We also show the range of the hyperparameters
used and the initial values. In addition to this, Table 6 enlist the final
optimized distributions of the STDP and the LIF parameters obtained using BO.

Table 5: The list of parameter settings for the Bayesian Optimization-based
hyperparameter search

Parameter | Initial Value | Range
---|---|---
$\eta$ | 10 | (0,50)
$\gamma$ | 5 | (0,10)
$\zeta$ | 2.5 | (0,10)
$\eta^{*}$ | 1 | (0,3)
$g$ | 2 | (0,10)
$\omega$ | 0.5 | (0,1)
$k$ | 50 | (0,100)
$\lambda$ (SHD) | 1 | (0,2)
$\lambda$ (Lorenz) | 1.5 | (0,4)
$P_{IR}$ | 0.05 | (0,0.1)
$\tau_{n-E},\tau_{n-I}$ (SHD) | $50ms$ | $(0ms,100ms)$
$\tau_{n-E},\tau_{n-I}$ (Lorenz) | $100ms$ | $(0ms,300ms)$
$A_{en-R},A_{EE},A_{EI},A_{IE},A_{II}$ | 30 | (0,60)

Table 6:  Table showing the average final distributions of the hyperparameters

| Parameter | Distribution |
---|---|---|---
STDP Parameter | $\tau_{+}$ | Normal |  | $\bar{\mu}=18.235$
---
$\bar{\sigma}=1.522$
$\tau_{-}$ | Normal |  | $\bar{\mu}=22.382$
---
$\bar{\sigma}=1.768$
$\eta_{+}$ | Normal |  | $\bar{\mu}=0.516$
---
$\bar{\sigma}=0.0055$
$\eta_{-}$ | Normal |  | $\bar{\mu}=0.448$
---
$\bar{\sigma}=0.0057$
LIF Parameter | $\tau_{m}^{(e)}$ | Gamma |  | $\bar{\alpha}=2.89$
---
$1/\bar{\beta}=0.248$
$\tau_{m}^{(i)}$ | Gamma |  | $\bar{\alpha}=5.14$
---
$1/\bar{\beta}=0.313$

Figure 6: Figure Showing the convergence behaviors of the three types of BO
described in the paper (a) BO optimizing the memory capacity $\mathcal{C}$ (b)
BO optimizing the average spike count $\tilde{S}$ and (c) BO optimizing the
spike efficiency $\mathcal{E}$

#####  A.3.2 Convergence Analysis

We compare the convergence analysis of the three Bayesian Optimization
techniques and the results are shown in Fig. 6. Each of the experiments was
repeated five times and the mean and variance of the observations are shown in
the Figure. It is to be noted here that since we define the BO as a
minimization principle, we minimize $\frac{1}{\mathcal{C}}$, $\tilde{S}$ and
$\frac{1}{\mathcal{E}}$.

####  A.4 Comparing Bayesian Optimization Objective Functions

We show the results of Bayesian Optimization results for the three cases we
are considering in this paper for both the classification and prediction
problems. The results for the classification problem are shown in Table 7. We
tabulate the memory capacity, the average spike count and the observed
accuracy for the three BO cases. Similarly, the results for the prediction
problem are shown in Table 8. In that case, we tabulate the memory capacity,
the average spike count and the observed NRMSE for the three BO cases. We
rerun each of the experiments 5 times and report the mean and standard
deviation of the results obtained.

Table 7: Table showing the performance of the Bayesian Optimization on the SHD
Classification dataset for the three different cases where BO 1 optimizes
$\mathcal{C}$ , BO 2 optimizes $\tilde{S}$ and BO 3 optimizes $\mathcal{E}$

N_R |  | BO 1
---
Memory
Capacity
| BO 1
---
Average Spike
Count
| BO 1
---
Accuracy
| BO 2
---
Memory
Capacity
| BO 2
---
Average Spike
Count
| BO 2
---
Accuracy
| BO 3
---
Memory
Capacity
| BO 3
---
Average Spike
Count
| BO 3
---
Accuracy
100 | $5.31\pm 0.36$ | $1399.14\pm 113.66$ | $67.41\pm 4.86$ | $5.22\pm 0.58$ | $1189.37\pm 86.95$ | $66.23\pm 5.74$ | $5.12\pm 0.31$ | $1268.86\pm 91.15$ | $68.67\pm 4.97$
200 | $5.45\pm 0.39$ | $1479.86\pm 133.72$ | $67.85\pm 4.16$ | $5.38\pm 0.53$ | $1233.19\pm 122.56$ | $67.21\pm 5.49$ | $5.28\pm 0.34$ | $1328.18\pm 131.45$ | $69.54\pm 4.27$
300 | $6.72\pm 0.4$ | $1541.14\pm 148.83$ | $68.41\pm 4.87$ | $6.42\pm 0.57$ | $1325.25\pm 128.47$ | $67.95\pm 5.11$ | $6.54\pm 0.36$ | $1391.68\pm 140.37$ | $70.1\pm 4.69$
400 | $7.21\pm 0.44$ | $1682.68\pm 239.95$ | $70.11\pm 4.33$ | $6.91\pm 0.51$ | $1425.69\pm 129.57$ | $68.23\pm 5.27$ | $7.01\pm 0.35$ | $1511.79\pm 200.96$ | $71.54\pm 4.06$
500 | $8.59\pm 0.48$ | $1768.24\pm 287.94$ | $71.05\pm 3.98$ | $7.69\pm 0.46$ | $1555.29\pm 139.17$ | $69.31\pm 5.03$ | $8.63\pm 0.38$ | $1621.8\pm 250.46$ | $73.05\pm 3.87$
1000 | $11.22\pm 0.46$ | $2251.17\pm 319.75$ | $72.93\pm 3.38$ | $9.03\pm 0.44$ | $2015.24\pm 147.44$ | $70.89\pm 4.91$ | $12.25\pm 0.39$ | $2102.59\pm 279.86$ | $75.32\pm 3.44$
2000 | $13.3\pm 0.51$ | $2566.21\pm 348.68$ | $75.36\pm 3.29$ | $9.89\pm 0.46$ | $2314.59\pm 151.18$ | $72.33\pm 4.88$ | $13.95\pm 0.37$ | $2410.08\pm 301.57$ | $77.25\pm 3.17$
3000 | $14.47\pm 0.53$ | $2825.47\pm 355.87$ | $77.14\pm 3.38$ | $10.48\pm 0.41$ | $2623.41\pm 177.94$ | $74.63\pm 4.93$ | $14.88\pm 0.42$ | $2708.52\pm 315.34$ | $78.21\pm 3.24$
4000 | $15.17\pm 0.52$ | $3551.07\pm 366.19$ | $78.05\pm 3.25$ | $11.57\pm 0.45$ | $3045.28\pm 225.53$ | $75.15\pm 4.85$ | $15.87\pm 0.38$ | $3218.42\pm 328.19$ | $79.36\pm 3.13$
5000 | $15.64\pm 0.57$ | $4186.49\pm 383.09$ | $78.92\pm 3.31$ | $11.68\pm 0.48$ | $3294.62\pm 241.14$ | $75.87\pm 4.81$ | $16.03\pm 0.41$ | $3573.51\pm 331.18$ | $80.49\pm 3.15$

Table 8: Table showing the performance of the Bayesian Optimization on the
Lorenz System Prediction dataset for the three different cases where BO 1
optimizes $\mathcal{C}$ , BO 2 optimizes $\tilde{S}$ and BO 3 optimizes
$\mathcal{E}$

N_R |  | BO 1
---
Memory
Capacity
| BO 1
---
Average Spike
Count
| BO 1
---
RMSE
| BO 2
---
Memory
Capacity
| BO 2
---
Average Spike
Count
| BO 2
---
RMSE
| BO 3
---
Memory
Capacity
| BO 3
---
Average Spike
Count
| BO 3
---
RMSE
100 | $3.56\pm 0.36$ | $1354.35\pm 108.96$ | $0.617\pm 0.019$ | $3.02\pm 0.52$ | $1101.25\pm 75.93$ | $0.684\pm 0.026$ | $3.45\pm 0.3$ | $1207.52\pm 67.26$ | $0.654\pm 0.0207$
200 | $4.12\pm 0.39$ | $1443.59\pm 127.23$ | $0.587\pm 0.02$ | $3.89\pm 0.48$ | $1207.35\pm 118.57$ | $0.639\pm 0.027$ | $4.01\pm 0.33$ | $1302.47\pm 96.05$ | $0.613\pm 0.0218$
300 | $5.26\pm 0.37$ | $1499.62\pm 141.73$ | $0.503\pm 0.027$ | $4.44\pm 0.55$ | $1257.26\pm 1257.26$ | $0.558\pm 0.034$ | $5.15\pm 0.34$ | $1335.81\pm 168.02$ | $0.531\pm 0.0287$
400 | $6.37\pm 0.41$ | $1528.73\pm 228.87$ | $0.459\pm 0.033$ | $5.87\pm 0.49$ | $1304.35\pm 126.18$ | $0.467\pm 0.04$ | $6.05\pm 0.35$ | $1415.32\pm 213.49$ | $0.482\pm 0.0347$
500 | $7.25\pm 0.45$ | $1601.27\pm 277.97$ | $0.389\pm 0.036$ | $6.25\pm 0.45$ | $1365.35\pm 137.04$ | $0.421\pm 0.043$ | $6.87\pm 0.36$ | $1507.29\pm 219.58$ | $0.411\pm 0.0377$
1000 | $10.12\pm 0.43$ | $1868.14\pm 301.17$ | $0.316\pm 0.04$ | $7.41\pm 0.51$ | $1563.25\pm 146.76$ | $0.396\pm 0.047$ | $9.03\pm 0.37$ | $1699.27\pm 275.79$ | $0.332\pm 0.0417$
2000 | $11.84\pm 0.48$ | $2105.95\pm 331.54$ | $0.293\pm 0.045$ | $8.02\pm 0.5$ | $1854.35\pm 150.28$ | $0.351\pm 0.052$ | $11.44\pm 0.39$ | $2014.12\pm 280.03$ | $0.301\pm 0.0467$
3000 | $13.71\pm 0.51$ | $2408.35\pm 348.26$ | $0.258\pm 0.042$ | $8.94\pm 0.53$ | $2195.82\pm 179.75$ | $0.335\pm 0.058$ | $13.87\pm 0.4$ | $2236.59\pm 281.05$ | $0.241\pm 0.0482$
4000 | $14.15\pm 0.52$ | $2951.56\pm 352.66$ | $0.242\pm 0.063$ | $9.55\pm 0.55$ | $2445.31\pm 217.73$ | $0.326\pm 0.07$ | $14.63\pm 0.41$ | $2546.25\pm 289.81$ | $0.227\pm 0.0649$
5000 | $14.45\pm 0.54$ | $3784.44\pm 353.51$ | $0.203\pm 0.064$ | $9.96\pm 0.58$ | $2684.59\pm 234.63$ | $0.302\pm 0.071$ | $15.12\pm 0.42$ | $2898.27\pm 307.14$ | $0.195\pm 0.0655$

####  A.5 Comparing the Generalizability

We observed that increasing the neuronal heterogeneity increases the memory
capacity of the network. However, this increment in the memory capacity might
lead to a model which overfits the training data. However, the heterogeneous
STDP model with varying synaptic dynamics gives rise to a heavy-tailed Feller
process. Recent works Simsekli et al. (2020), Chakraborty & Mukhopadhyay
(2021) show that the Hausdorff dimension of the trajectories of the sample
paths of the learning algorithm can control the generalization error. This is
intimately linked to the tail behavior of the driving process. The authors
showed that heavier-tailed processes achieve better generalization. Thus, the
tail index of the process can be used as a notion of capacity metric that
estimates the generalization error, which does not necessarily grow with the
number of parameters. The authors discuss that the stochastic process for the
synaptic weights behaves like a Lévy motion around a local point. Because of
this locally regular behavior, the Hausdorff dimension can be bounded by the
Blumenthal-Getoor (BG) index (Blumenthal & Getoor, 1960), which depends on the
tail behavior of the Lévy process. Thus, we can use the BG index as a bound
for the Hausdorff dimension of the trajectories from the STDP learning
process. Now, as the Hausdorff dimension is a measure of the generalization
error and is also controlled by the tail behavior of the process, heavier
tails imply less generalization error. In this paper, we empirically study the
generalization ability of the HRSNN network using the BG index as a metric. We
did the experiments on the 4 ablation study models for the classification task
on the SHD dataset, and the results are reported in Table 9. From the table,
we see that the heterogeneity in STDP improves the generalization error the
most, while the heterogeneity in the LIF neurons increases the training and
testing accuracies.

Table 9:  Table showing the Ablation Study for the comparison of the
Generalizability of heterogeneous networks

| BG Index |  | Training
---
Accuracy (A)
| Testing
---
Accuracy (B)
| Generalization Error
---
(|A-B|)
| Hom LIF
---
Hom STDP
1.522 | 87.33 | 73.58 | 13.75
| Hom LIF
---
Het STDP
1.438 | 85.31 | 74.03 | 11.28
| Het LIF
---
Hom STDP
1.835 | 95.29 | 78.87 | 16.42
| Het LIF
---
Het STDP
1.711 | 94.32 | 80.49 | 13.83

####  A.6 Results on Limited Training Data

Table 10: Table showing results with limited training data

| Percentange
---
Training Data
| Train
---
Accuracy
(A)
| Test
---
Accuracy
(B)
| Generalization
---
Error
|A-B|
| Train
---
Accuracy
(A)
| Test
---
Accuracy
(B)
| Generalization
---
Error
|A-B|
|  | Heterogeneous LIF,
---
Heterogeneous STDP
| Homogeneous LIF,
---
Homogeneous STDP
100 | 94.32 | 80.49 | 13.83 | 87.33 | 73.58 | 13.75
90 | 94.89 | 78.34 | 16.55 | 87.83 | 67.84 | 19.99
80 | 95.47 | 76.72 | 18.75 | 88.86 | 65.86 | 23
70 | 96.15 | 74.92 | 21.23 | 89.95 | 62.19 | 27.76
60 | 96.92 | 70.34 | 26.58 | 91.58 | 61.25 | 30.33
50 | 97.69 | 69.44 | 28.25 | 94.38 | 59.51 | 34.87
40 | 98.21 | 63.76 | 34.45 | 96.85 | 55.93 | 40.92
30 | 98.43 | 54.01 | 44.42 | 98.43 | 45.86 | 52.57
20 | 99.43 | 43.87 | 55.56 | 99.49 | 42.68 | 56.81
10 | 100 | 31.43 | 68.57 | 100 | 30.18 | 69.82
5 | 100 | 15.32 | 84.68 | 100 | 14.38 | 85.62
|  | Heterogeneous LIF,
---
Homogeneous STDP
| Homogeneous LIF,
---
Heterogeneous STDP
100 | 97.29 | 78.87 | 18.42 | 86.31 | 74.03 | 12.28
90 | 97.41 | 77.48 | 19.93 | 86.94 | 68.59 | 18.35
80 | 97.65 | 76.32 | 21.33 | 87.75 | 67.58 | 20.17
70 | 97.95 | 74.03 | 23.92 | 88.17 | 65.25 | 22.92
60 | 98.03 | 71.16 | 26.87 | 89.52 | 63.11 | 26.41
50 | 98.43 | 68.48 | 29.95 | 90.48 | 60.86 | 29.62
40 | 98.79 | 61.93 | 36.86 | 93.15 | 57.31 | 35.84
30 | 99.56 | 51.68 | 47.88 | 96.34 | 48.41 | 47.93
20 | 100 | 44.52 | 55.48 | 98.43 | 43.59 | 54.84
10 | 100 | 30.68 | 69.32 | 99.56 | 32.57 | 66.99
5 | 100 | 14.15 | 85.85 | 100 | 18.48 | 81.52

We have trained the models with limited training data. We observe that the
HRSNN model with heterogeneous LIF and STDP dynamics not only has better
testing accuracy but also shows better generalization behavior when compared
to other homogeneous RSNN or the other ablation heterogeneous models (with
heterogeneity in only one of them). Also, we see that the HRSNN model with
heterogeneous STDP shows distinctly better generalization ability than the
generalization ability of HRSNN with heterogeneous LIF neurons. On the other
hand, the latter showcases significantly higher training and testing accuracy
compared to the former model. This can be interpreted as follows: since
heterogeneous LIF dynamics increase the memory capacity, it leads to an
overfitting of the data. Heterogeneous STDP dynamics help in obtaining more
generalizable solutions from this. Each has its own downsides; however, using
HRSNN with both heterogeneous LIF and STDP dynamics shows better performance
and generalization abilities, as seen from Table 10.

####  A.7 Further Evaluations

In Section B, we argued that as the heterogeneity in the neuronal parameters
increases, the covariance decreases; hence the neurons become less correlated.
In this section, we give empirical results to support the theory. We tested
the model on more complex datasets - (i) The Spiking Heidelberg Digits (SHD)
dataset (ii) the Spiking Speech Command (SSC) dataset are both audio-based
classification datasets for which input spikes and output labels are provided
Cramer et al. (2020) and (iii) CIFAR10 DVS dataset Li et al. (2017).

  * •

Impact of Heterogeneity on Covariance: We plot the covariance matrices for
different levels of heterogeneity $\mathcal{J}$ (Eq. 46) for a small network
with 50 neurons. The covariance matrix is calculated by taking the average
neuronal states before the appearance of the first spike in the final layer.
We see that as the heterogeneity in the neuronal parameters increases, the
correlation between the neurons decreases. The results are shown in Fig. 7

  * •

Impact of Heterogeneity on Principal Components: From the covariance plots, we
see that increasing $\mathcal{J}$. reduces the correlation between neurons. We
also plot the probability density functions of the eigenvalues of the
covariance matrix of the neurons with increasing heterogeneity in the neuronal
parameters. We see that with higher heterogeneity in the neuronal parameters
$\mathcal{J}$, the distribution of the eigenvalues of the covariance becomes
flatter. This signifies that the covariance matrix has a lower variance for
higher $\mathcal{J}$. A flatter distribution also indicates that a larger
number of principal components are active. This supports our hypothesis that
heterogeneity in the neuronal parameters increases the number of principal
components and helps increase the model’s memory capacity. The result is shown
in this Fig. 8

  * •

Impact of Heterogeneity in STDP on Firing Rate: We plot the mean firing rate
of the neurons for the four types of HRSNNs and MRSNN with homogeneous LIF and
STDP dynamics. We plot the results for a smaller network with 100 neurons and
a Poisson input process. The MRSNN model shows a much higher firing rate,
especially at a higher frequency, demonstrating that MRSNN requires
significantly more spikes than the HRSNN model. The result is shown in Fig. 9

  * •

Coupling Strength: We note here that in this paper, we use (homogeneous or
heterogeneous) STDP to learn the synaptic conductance connecting various
neurons in the SNN. Therefore, we do not control the synaptic coupling
strength as independent variables and hence, cannot perform control
experiments with various extents of coupling strength. An interesting future
extension of the results will be quantifying the coupling strength for HRSNN
with heterogeneity in LIF and STDP dynamics. We can leverage McKenzie et
al.McKenzie et al. (2021), where the authors proposed statistical tools to
estimate synaptic coupling dynamics from spike-spike correlations.

Figure 7: Figure showing as the heterogeneity in the neuronal parameters
increases, the covariance between the neurons decreases Figure 8: With higher
heterogeneity in the neuronal parameters $\mathcal{H}$, the distribution of
the eigenvalues of the covariance becomes flatter. This signifies that the
covariance matrix has a lower variance for higher $\mathcal{H}$. A flatter
distribution also signifies a greater number of principal components are
active, which supports our hypothesis that heterogeneity in the neuronal
parameters increases the number of principal components and helps in
increasing the memory capacity of the model. Figure 9: Figure showing the
histograms of the firing rates for the four kinds of heterogeneous RSNN with
homogeneous /heterogeneous neurons and synapses

| SHD | SSC | CIFAR10 DVS
---|---|---|---
| Training
---
Accuracy
($A$)
| Testing
---
Accuracy
($B$)
| Generalization
---
Error
$|A-B|$
| Training
---
Accuracy
($A$)
| Testing
---
Accuracy
($B$)
| Generalization
---
Error
$|A-B|$
| Training
---
Accuracy
($A$)
| Testing
---
Accuracy
($B$)
| Generalization
---
Error
$|A-B|$
Hom LIF Hom STDP | $86.92\pm 1.35$ | $72.89\pm 1.85$ | $14.03\pm 1.67$ | $74.69\pm 1.72$ | $47.94\pm 1.94$ | $26.75\pm 1.42$ | $82.41\pm 1.8$ | $65.33\pm 3.41$ | $17.08\pm 1.35$
Hom LIF Het STDP | $85.76\pm 1.27$ | $73.91\pm 1.49$ | $11.85\pm 1.25$ | $76.79\pm 1.58$ | $52.96\pm 1.73$ | $23.86\pm 1.29$ | $83.48\pm 1.52$ | $67.06\pm 2.97$ | $16.42\pm 1.24$
Het LIF Hom STDP | $95.29\pm 1.16$ | $78.36\pm 1.42$ | $16.93\pm 1.13$ | $84.26\pm 1.33$ | $55.11\pm 1.65$ | $29.15\pm 1.12$ | $86.93\pm 1.79$ | $68.37\pm 3.05$ | $18.56\pm 1.42$
Het LIF Het STDP | $94.07\pm 1.03$ | $80.01\pm 1.13$ | $14.06\pm 1.02$ | $86.41\pm 1.49$ | $59.28\pm 1.35$ | $27.13\pm 0.97$ | $87.49\pm 1.76$ | $70.54\pm 1.82$ | $16.95\pm 1.38$

###  Appendix B Supplementary Section B

####  B.1 Approximations and assumptions

We make several approximations and assumptions for this section’s theoretical
analysis of the heterogeneous RSNN networks. Firstly, it must be noted that in
this paper, the analytical relations are derived by taking the heterogeneity
individually. i.e., when we consider heterogeneity in the neuronal parameters,
we consider homogeneous STDP dynamics and vice-versa. In addition to this, we
assume diffusion approximation. That is, if a neuron receives Poissonian
uncorrelated input spike trains and the contribution of a single synaptic
connection is small compared to the distance between reset and threshold
$w\ll\left(V_{\Theta}-V_{0}\right)$, the random input can be approximated by
Gaussian white noise with mean $\mu$ and noise intensity $\sigma^{2}$. This
approximation does not hold if the network features highly correlated activity
or receives strong external input common to many neurons. Also, we assume a
fast/slow synaptic regime in which the synaptic time constant
$\tau_{\mathrm{s}}$ is much shorter/longer than the membrane time constant
$\tau_{\mathrm{m}}$. In this work, we consider a mean-field approximation of
the HRSNN network with heterogeneity in the parameters of the LIF neurons and
the STDP dynamics independently.

####  B.2 Mean-Field Reduction Model of HRSNN

In this section, we model the HRSNN network using heterogeneity in only the
LIF neuron parameters. Following the works of Ly et al.Ly (2015), we can write
the equations for the excitatory neurons indexed by $j\in$
$\left\\{1,2,\ldots,N_{e}\right\\}$ are:

$\displaystyle\tau_{m}\frac{dv_{j}}{dt}$
$\displaystyle=-v_{j}-g_{ie}(t)\left(v_{j}-\mathcal{E}_{I}\right)-g_{ee}(t)\left(v_{j}-\mathcal{E}_{E}\right)+\sigma_{E}\eta_{j}(t)$
(11) $\displaystyle v_{j}\left(t^{*}\right)$
$\displaystyle\geq\theta_{j}(\text{ refractory period })\Rightarrow
v_{j}\left(t^{*}+\tau_{ref}\right)=0$ (12)
$\displaystyle\tau_{n}\frac{d\eta_{j}}{dt}$
$\displaystyle=-\eta_{j}+\sqrt{\tau_{n}\xi_{j}(t)}$ (13) $\displaystyle
g_{ee}(t)$
$\displaystyle=q_{j}\frac{\gamma_{ee}}{p_{ee}N_{e}}\sum_{j^{\prime}\in\\{\text{
presyn E cells\\} }}G_{j^{\prime}}(t)$ (14) $\displaystyle g_{ei}(t)$
$\displaystyle=\frac{\gamma_{ei}}{p_{ei}N_{i}}\sum_{k^{\prime}\in\\{\text{
presyn I cells\\} }}G_{k^{\prime}}(t)$ (15)
$\displaystyle\tau_{d}\frac{dG_{j}}{dt}$ $\displaystyle=-G_{j}+A_{j}$ (16)
$\displaystyle\tau_{r}\frac{dA_{j}}{dt}$
$\displaystyle=-A_{j}+\tau_{r}\alpha\sum_{l}\delta\left(t-t_{l}\right)$ (17)

where the inhibitory and excitatory reversal potentials are $\mathcal{E}_{I}$,
and $\mathcal{E}_{E}$, respectively, with $\mathcal{E}_{I}<0<\mathcal{E}_{E}$.
$\xi_{j}(t)$ are uncorrelated white noise processes, $p_{xy}$ is the
proportion of neuron type $y$ (randomly chosen) that provides presynaptic
input to neuron type $x$ $(x,y\in\\{e,i\\})$. The second line in the equations
describes the refractory period at spike time $t^{*}$. When the neuron’s
voltage crosses threshold $\theta_{j}$, the neuron goes into a refractory
period for $\tau_{ref}$ where the voltage is undefined, after which we set the
neuron’s voltage to 0. In the last equation, $t_{l}$ denotes the spike times
of the $j$ th excitatory neuron. Now, for the mean-field analysis, we use
$q_{ji}$ to model the synaptic heterogeneity between the pre-and post-synaptic
neurons by modulating the synaptic conductance for both the excitatory and
inhibitory neurons.

We note here the numerical assumptions for the mean-field analysis:

  1. 1.

finite size effects are negligible ($\mathrm{N}$ e/i $\gg 1$ )

  2. 2.

the firing rate of presynaptic neurons is governed by a Poisson process

  3. 3.

the population firing rate averaged over $q$ and $\tau_{m}$ is a good
approximation to the average presynaptic input rate and

  4. 4.

a single p.d.f. function is sufficient to describe the population behavior)
(finite $N$)

Similarly, for the inhibitory neurons indexed by $k\in$
$\left\\{1,2,\ldots,N_{i}\right\\}$, the equations are:

$\displaystyle\tau_{m}\frac{dv_{k}}{dt}$
$\displaystyle=-v_{k}-g_{ii}(t)\left(v_{k}-\mathcal{E}_{I}\right)-g_{ei}(t)\left(v_{k}-\mathcal{E}_{E}\right)+\sigma_{I}\eta_{k}(t)$
(18) $\displaystyle v_{k}\left(t^{*}\right)$ $\displaystyle\geq 1(\text{
refractory period })\Rightarrow v_{j}\left(t^{*}+\tau_{ref}\right)=0$ (19)
$\displaystyle\tau_{n}\frac{d\eta_{k}}{dt}$
$\displaystyle=-\eta_{k}+\sqrt{\tau_{n}\xi_{k}(t)}$ (20) $\displaystyle
g_{ie}(t)$
$\displaystyle=q_{j}\frac{\gamma_{ie}}{p_{ie}N_{e}}\sum_{k^{\prime}\in\\{\text{
presyn I cells\\} }}G_{k^{\prime}}(t)$ (21) $\displaystyle g_{ii}(t)$
$\displaystyle=\frac{\gamma_{ii}}{p_{ii}N_{i}}\sum_{k^{\prime}\in\text{
\\{presyn I cells\\} }}G_{k^{\prime}}(t)$ (22)
$\displaystyle\tau_{d}\frac{dG_{k}}{dt}$ $\displaystyle=-G_{k}+A_{k}$ (23)
$\displaystyle\tau_{r}\frac{dA_{k}}{dt}$
$\displaystyle=-A_{k}+\tau_{r}\alpha\sum_{l}\delta\left(t-t_{l}\right)$ (24)

Please refer to the paper by Ly et al. Ly (2015) for details regarding the
equations. Since the recurrent coupled stochastic network is difficult to
describe theoretically, we use population density methods, where an equation
determines the probability of a neuron being in a particular state. The
variables in the populations are determined using distribution functions. The
two forms of heterogeneity introduce a large number of dimensions. For
simplicity, one can track a family of probability density functions for each
$\left(q_{j},\tau_{j}\right)$ pair for each neuron. The subsequent equations
are a good approximation to the HRSNN network with the following assumptions:
(i) finite size effects are negligible $\left(N_{e/i}\gg 1\right)$ (ii) the
firing rate of presynaptic neurons is governed by a Poisson process (iii) the
population firing rate averaged over $\bm{q}$ and $\bm{\tau}_{m}$ is a good
approximation to the average presynaptic input rate (iv) a single p.d.f.
function is sufficient to describe the population behavior, and the
heterogeneity is driven by $\left(q_{j},\tau_{m},j\right)$ For each pair of
values $\left(q_{j},\tau_{m},j\right)$, the probability density function
$\rho$ is defined by:

$\displaystyle\int_{\Omega}\rho\left(v_{E},\mathbf{w}_{E},v_{I},\mathbf{w}_{I},t\right)dv_{E}d\mathbf{w}_{E}dv_{I}d\mathbf{w}_{I}$
$\displaystyle=\operatorname{Pr}((v_{E}(t),\mathbf{w}_{E}(t),v_{I}(t),\mathbf{w}_{I}(t))$
$\displaystyle\in\Omega)$ (25)

where $\mathbf{w}_{X}$ denotes the other states variables of the corresponding
neuron type $X\in\\{E,I\\}$, consisting of conductance, colored noise:
$\mathbf{w}_{X}=\left(g_{X},a_{X},\eta_{X}\right)$. The evolution of the
p.d.f.’s is governed by a continuity equation and boundary conditions:

$\displaystyle\frac{\partial\rho}{\partial t}=-\nabla\cdot\mathbf{J}$ (26)
$\displaystyle\mathbf{J}:=\left(J_{v_{E}},J_{g_{E}},J_{a_{E}},J_{\eta_{E}},J_{v_{I}},J_{g_{I}},J_{a_{I}},J_{\eta_{I}}\right)$
(27) $\displaystyle
J_{v_{E}}:=-\frac{1}{\tau_{m}}\left[v_{E}+q\gamma_{ei}g_{I}\left(v_{E}-\mathcal{E}_{I}\right)+q\gamma_{ee}g_{E}\left(v_{E}-\mathcal{E}_{E}\right)+\sigma_{E}\eta_{E}\right]\rho$
(28) $\displaystyle
J_{v_{I}}:=-\frac{1}{\tau_{m}}\left[v_{I}+\gamma_{ii}g_{I}\left(v_{I}-\mathcal{E}_{I}\right)+\gamma_{ie}g_{E}\left(v_{I}-\mathcal{E}_{E}\right)+\sigma_{I}\eta_{I}\right]\rho$
(29) $\displaystyle
J_{g_{X}}:=-\frac{1}{\tau_{d}}\left[g_{X}-a_{X}\right]\rho$ (30)
$\displaystyle
J_{a_{X}}:=-\frac{a_{X}}{\tau_{r}}+v_{X}(t)\int_{a_{X}-\alpha_{X}}^{a_{X}}\rho\left(\ldots,a_{X}^{\prime},\ldots\right)da_{X}^{\prime}$
(31) $\displaystyle
J_{\eta_{X}}:=-\frac{1}{\tau_{n}}\eta_{X}\rho+\frac{1}{\tau_{n}}\frac{\partial^{2}\rho}{\partial\eta_{X}^{2}}$
(32) $\displaystyle
v_{X}(t):=\iiint\frac{1}{\tau_{m}}J_{v_{X}}d\mathbf{w}_{X}dqd\tau_{m}$ (33)
$\displaystyle J_{\mathbf{w}_{X}}\mid\partial\mathbf{w}_{X}=0$ (34)

The definitions of $g_{XY}$ in the LIF neuron equations defined above result
in a total conductance of $\gamma_{XY}g_{Y}$ on average.

We describe an insightful analytic reduction that captures how the range of
excitatory firing rates changes in different regimes. We focus on only the
excitatory neurons, which have fewer state variables if the inhibitory
population is ignored or assumed to be known.

Let us denote the approximate excitatory firing rate(s) $v_{E}$ as $r$. The
deterministic firing rate of the equation

$\tau_{m}\frac{dv_{E}}{dt}=-v_{E}-q\tilde{g_{I}}\left(v_{E}-\mathcal{E}_{I}\right)-q\tilde{g_{E}}\left(v_{E}-\mathcal{E}_{E}\right)+\tilde{\eta_{E}}$
(35)

is given by

$r_{0}(q,\tau_{m};\tilde{\mathbf{w}}_{E})=\begin{cases}0,&\text{ if
}\frac{q\left(\tilde{g_{E}}\mathcal{E}_{E}+\tilde{g_{I}}\mathcal{E}_{I}\right)+\tilde{\eta_{E}}}{1+q\left(\tilde{g_{E}}+\tilde{g_{I}}\right)}\leq\theta\\\
\frac{1+q\left(\tilde{g_{E}}+\tilde{g_{I}}\right)}{\tau_{m}\left(\tilde{g_{E}^{*}}\mathcal{E}_{E}+\tilde{g_{I}}\mathcal{E}_{I}\right)+\tilde{\eta_{E}}}&\text{
if
}\frac{q\left(\tilde{g_{E}}\mathcal{E}_{E}+\tilde{g_{I}}\mathcal{E}_{I}\right)+\tilde{\eta_{E}}}{1+q\left(\tilde{g_{E}}+\tilde{g_{I}}\right)}>\theta\end{cases}$
(36)

We define:
$\tilde{g_{E}}:=\gamma_{ee}g_{E},\tilde{g_{I}}:=\gamma_{ei}g_{I},\tilde{\eta_{E}}:=\sigma_{E}\eta_{E}$.
Finally, the given state variables are integrated against their marginal
density to get:

$r(q,\theta)=\mathbb{E}\left[\frac{r_{0}}{1+r_{0}\tau_{ref}}\right]=\int\frac{r_{0}}{1+r_{0}\tau_{ref}}\tilde{\rho}\left(\tilde{g_{E}},\tilde{g_{I}},\tilde{\eta_{E}}\right)d\tilde{\mathbf{w}_{E}}$
(37)

There is a slight abuse of notation because the auxiliary variables $a_{X}$
effect the conductances but are not written in the previous equation; the
emphasis is on how $\left(\tilde{g_{E}},\tilde{g_{I}},\tilde{\eta_{E}}\right)$
directly effects $r$. Since the external noise is applied indiscriminately,
$\tilde{\eta_{E}}$ is independent of the other variables and the marginal
density factors into:

$\tilde{\rho}\left(\tilde{g_{E}},\tilde{g_{I}},\tilde{\eta_{E}}\right)=\tilde{\rho}\left(\tilde{g_{E}},\tilde{g_{I}}\right)\frac{e^{-\left(\tilde{\eta_{E}}/\sigma_{E}\right)^{2}}}{\sigma_{E}\sqrt{\pi}}$
(38)

However, $\tilde{\rho}\left(\tilde{g_{E}},\tilde{g_{I}}\right)$ is still not
analytically tractable, leading us to rely on Monte Carlo simulations to
numerically estimate $\tilde{\rho}\left(\tilde{g_{E}},\tilde{g_{I}}\right)$.

It must be noted here that this is a reduction model for the HRSNN network
with many simplifying assumptions. It is not a complete mean-field derivation
of the HRSNN model with heterogeneous LIF neurons, and heterogeneous STDP
dynamics is a fascinating research question but beyond the scope of this
paper.

####  B.3 Analytical Results

Figure 10:  (a)Figure showing the variation of memory capacity with neuronal
heterogeneity (b) Figure showing the variation of efficiency $\mathcal{E}$
with number of neurons $N_{\mathcal{R}}$

Analytical Results of Memory Capacity Neuroscience networks of spiking neurons
are increasingly used to understand mechanisms underlying phenomena observed
in electrophysiological recordings. There are two complementary strategies for
studying such a recurrent network of spiking neurons - (a) numerical
simulations and (b) analytical methods using mean field models. With numerical
simulations, we can simulate any network model without any approximation.
However, this method typically works in high dimensional parameter space and
is, thus, hard to interpret. Also, it is generally hard to characterize
parameter regions where specific behaviors are found using numerical
simulations. On the other hand, with analytical calculations, we obtain deeper
insights into mechanisms underlying specific behaviors and can obtain critical
parameters that control specific behaviors. So, now, we analytically study the
variance of the estimated memory capacity with the change in the heterogeneity
of neuronal parameters. We plot the change in the estimated memory capacity
$\mathcal{C}$, calculated using Eq. 3. We plot this with respect to the
neuronal heterogeneity $\mathcal{H}$, measured using the entropy of the
neuronal parameters for the HRSNN model. The result is plotted in Fig. 10(a).
We use a HRSNN model with $N_{\mathcal{R}}=1000$ and sequences of 4,000 random
inputs chosen from $\mathcal{U}[-1;1]$. We see that, as predicted, the memory
capacity of the model increases linearly with the increase in heterogeneity
within the limits of the application, as proved in Theorem 1. The error bars
in Fig. 10(a) represent the standard deviation of the observations.

Analytical Study of Spike Efficiency We calculate the average firing rate of
the heterogeneous spiking neural network for the prediction task during
inference, and the results are shown in Fig 10(b). Using heterogeneity in the
STDP parameters reduces the average number of spiking activations while
keeping the memory capacity almost equal. This result shows that Heterogeneous
STDP leads to sparse activation of neurons, as proved in Theorem 2.

Comparison with Neuroscience Works:  We compare the analytical results
obtained with some of the standard recurrent LIF network models in the
literature. Brunel et al.Brunel (2000) analytically study the dynamics of
sparsely connected a network of sparsely connected excitatory and inhibitory
integrate-and-fire neurons. The authors showed the existence of a diverse set
of states, including synchronous states in which neurons fire regularly;
asynchronous states with stationary global activity and very irregular
individual cell activity; and states in which the global activity oscillates
but individual cells fire irregularly, typically at rates lower than the
global oscillation frequency. In this paper, we use heterogeneity in the LIF
neurons. This leads to a diverse set of states for the neurons, which
consequently helps orthogonalize the state space dynamics to increase the
information stored in the memory of the network.

Deneve et al. Denève & Machens (2016) discussed the inefficiency of irregular
Poisson rate encoding in the brain. The authors argue that the Poisson point
process, which we use to model the spike firing rate, is extremely inefficient
as it exponentially increases the number of spikes required to convey
information. The authors further discuss that a continuum exists between
loosely balanced and tightly balanced spike-coding networks in neuroscience.
Though loosely balanced networks are inefficient, they are cheap in terms of
the number of connections per neuron and structure (Boerlin et al., 2013;
Boerlin & Denève, 2011; Bourdoukan et al., 2012). On the other hand, tightly-
balanced spike-coding networks are highly efficient but extremely structured,
dense connections that STDP rules must constantly maintain. For the HRSNN
model, since we are engineering an artificial spiking neural network model,
our network is highly structured and constantly updated using the
heterogeneous STDP rules. Thus, we might say that the HRSNN model is a
tightly-coupled network that helps in an efficient transfer of information.
This hypothesis is supported by the results shown in Table 2, where the HRSNN
model shows a higher performance using a lesser number of spikes.

####  B.4 Memory Capacity

Let $x(t)\in U$ (where $-\infty<t<+\infty$ and $U\subset\mathbb{R}$ is a
compact interval) be a single-channel stationary input signal. Assume that we
have an RSNN, specified by its internal weight matrix $\mathbf{W}$, its input
weight vector $\mathbf{w}^{\text{in }}$ and the unit output functions
$\mathbf{f},\mathbf{f}^{\text{out }}$. The network receives $x(t)$ at its
input unit. For a given delay $\tau$ and an output unit $y_{\tau}$ with
connection weight vector $\mathbf{w}_{\tau}^{\text{out }}$ we consider the
determination coefficient

$\displaystyle d\left[\mathbf{w}_{\tau}^{\text{out
}}\right]\left(x(t-\tau),y_{\tau}(t)\right)=$
$\displaystyle=d\left(x(t-\tau),\mathbf{w}_{\tau}^{\text{out
}}\left(\begin{array}[]{c}x(t)\\\ \mathbf{r}(t)\end{array}\right)\right.$
$\displaystyle=\frac{\operatorname{Cov}^{2}\left(x(t-\tau),y_{\tau}(t)\right)}{\sigma^{2}(x(t))\sigma^{2}\left(y_{\tau}(t)\right)}$

where Cov denotes covariance and $\sigma^{2}$ variance. The $\tau$-delay
Memory capacity of the network is defined by
$\displaystyle\mathcal{C}_{\tau}=\max_{\mathbf{w}_{\tau}^{\text{out
}}}d\left[\mathbf{w}_{\tau}^{\text{out
}}\right]\left(x(t-\tau),y_{\tau}(t)\right)$. The Memory capacity of the
network is
$\displaystyle\mathcal{C}=\sum_{\tau=1}^{\infty}\mathcal{C}_{\tau}.$ The
determination coefficient of two signals is the squared correlation
coefficient. It ranges between 0 and 1 and represents the fraction of variance
explainable in one signal by the other. Thus, the Memory capacity measures how
much variance of the delayed input signal can be recovered from optimally
trained output units, summed over all delays. Note that the output units do
not interfere; arbitrarily, many output units $y_{\tau}$ can be attached to
the same network.

The performance of the heterogeneous network model derives from its ability to
retain the memory of previous inputs. To quantify the relationship between the
recurrent layer dynamics and the memory capacity, we note that the extraction
of information from the recurrent layer is made through a linear combination
of the neurons’ states. Hence, more linearly independent neurons would offer
more variable states and, thus, more extended memory.

For reservoir computing (RC), Jaeger et al. Jaeger (2002) shows that
$\mathcal{C}$ is bounded by the reservoir network size of the linear RC with
the identity activation function and the independent and identically
distributed (i.i.d.) model input. Memory capacity ( $\mathcal{C}$) is used to
quantify the memory of RSNN. Such memory capacity measures the ability of RC
to reconstruct precisely the past information of the model input. Also, the
network’s structural properties can greatly impact the $\mathcal{C}$ of the
linear RC. Now, the question arises what is the need to maximize the memory
capacity of the network? The $\mathcal{C}$ normally serves as a global index
to quantify the memory property of the network. To comprehensively examine the
memory property deeply, the local measurement of its memory property is
indispensable. Thus, maximizing the $\mathcal{C}$ acts as an estimator for
better prediction results of the trained network.

Since the first-order approximation of the model is linear, the heterogeneity
between state variables depends on all the eigenvalues of the adjacency
matrix, with a larger mean eigenvalue meaning higher heterogeneity. Hence we
can use the eigenvalues $\left\\{\lambda_{i}\right\\}$ of the weight matrix
$W$ to quantify approximately how fast the input decays in the recurrent
layer. In other words, the eigenvalues of $W$ should be related to the memory
capacity of the heterogeneous neural network model. Indeed, we find that the
average eigenvalue modulus:
$\langle|\lambda|\rangle=1/N_{R}\sum_{i=1}^{N_{R}}\left|\lambda_{i}\right|$
strongly correlates with $\mathcal{H}$ and therefore with $\mathcal{C}$ as
well. Note that, instead of $\mathcal{C}$ and
$\mathcal{H},\langle|\lambda|\rangle$ is much easier to compute and is solely
determined by the recurrent layer network.

The memory capacity reflects the precision with which previous inputs can be
recovered. The nonlinearity of the recurrent layer and other far-in-the-past
inputs induce noise that complicates recovery. Thus, similar to the analysis
done by Aceituno et al.Aceituno et al. (2020) for Echo state networks, the
variance of the linear part of the recurrent layer is placed to maximize the
recoverable information. Thus, the inputs are projected into orthogonal
directions of the recurrent layer state space to not add noise to each other.
The variance spread across the different dimensions should be evenly
distributed within those orthogonal directions, quantified by the neurons’
covariance.

We start by noticing that the linear nature of the projection vector
$\mathbf{w}_{\text{out }}$ implies that we are treating the system as

$\mathbf{r}(t)=\sum_{\tau=0}^{\infty}\mathbf{a}_{\tau}x(t-\tau)+\varepsilon(t)$
(39)

where the vectors $\mathbf{a}_{\tau}\in\mathbb{R}^{N_{R}}$ correspond to the
linearly extractable effect of $x(t-\tau)$ onto $\mathbf{r}(t)$ and
$\varepsilon(t)$ is the nonlinear contribution of all the inputs onto the
state of $\mathbf{r}(t)$.

Previous works have shown that linear recurrent layers have more extended
memory, but nonlinearity is needed to perform interesting computations. Here
we show that for a fixed ratio of the nonlinearity, greater heterogeneity
leads to a lesser neuronal correlation, leading to a higher memory capacity.

To maintain this trade-off between linear and non-linear behavior, we will
assume that linear and non-linear strengths distribution is fixed. This can be
achieved if we impose that the probabilities of the neuron states do not
change, meaning that the mean, variance, and other moments of the neuron
outputs are unchanged; hence, the strength of the non-linear effects is
unchanged. A first constraint can also be obtained from the maintained
strength of the linear side of Eq.39

$\operatorname{Var}\left(\sum_{\tau=1}^{\infty}\mathbf{a}_{\tau}x(t-\tau)\right)=c$
(40)

where $c$ is a constant.

Lemma 3.1.1:  The state of the neuron can be written as follows:

$\displaystyle r_{i}(t)$
$\displaystyle=\sum_{k=0}^{N_{R}}\sum_{n=1}^{N_{R}}\lambda_{n}^{k}\left\langle
v_{n}^{-1},\mathbf{w}^{\mathrm{in}}\right\rangle\left(v_{n}\right)_{i}x(t-k)$
(41)

where $\mathbf{v}_{n},\mathbf{v}_{n}^{-1}\in\mathbf{V}$ are, respectively, the
left and right eigenvectors of $\mathbf{W}$, and
$\lambda_{n}^{k}\in\mathbf{\lambda}$ belongs to the diagonal matrix containing
the eigenvalues of $\mathbf{W}$;
$\mathbf{a}_{i}=\left[a_{i,0},a_{i,1},\ldots\right]$ represents the
coefficients that the previous inputs $\mathbf{x}_{t}=[x(t),x(t-1),\ldots]$
have on $r_{i}(t)$.

Proof:  We build on the work of Aceituno et al.Aceituno et al. (2020) where
they showed that higher heterogeneity among the neuronal states implies higher
memory capacity. Here we aim to show that as the number of neurons $N_{R}$ in
the recurrent layer decreases, heterogeneity increases the spectral radius.
More formally, the spectral radius $|\lambda_{n}|$ is directly proportional to
$\mathcal{H}$ as $N_{R}$ decreases. We express the state of a neuron
$r_{i}(t)$ as

$r_{i}(t)=\sum_{k=0}^{\infty}\left(W^{k}\mathbf{w}^{\mathrm{in}}\right)_{i}x(t-k)=\sum_{k=0}^{\infty}a_{i,k}x(t-k)=\left\langle\mathbf{a}_{i},\mathbf{x}_{t}\right\rangle$
(42)

where the vector $\mathbf{a}_{i}=\left[a_{i,0},a_{i,1},\ldots\right]$
represents the coefficients that the previous inputs
$\mathbf{x}_{t}=[x(t),x(t-1),\ldots]$ have on $r_{i}(t).$ We can then plug
this into the covariance between two neurons,

$\displaystyle\operatorname{Cov}\left(r_{i},r_{j}\right)$
$\displaystyle=\lim_{T\rightarrow\infty}\frac{1}{T}\sum_{q=t}^{t+T}\left\langle\mathbf{a}_{i},\mathbf{x}_{q}\right\rangle\left\langle\mathbf{a}_{j},\mathbf{x}_{q}\right\rangle$
$\displaystyle=\left\langle\mathbf{a}_{i},\mathbf{a}_{j}\right\rangle\lim_{T\rightarrow\infty}\frac{1}{T}\sum_{q_{i}=0}^{T}\sum_{q_{j}=0}^{T}\left\langle\mathbf{x}_{q_{i}},\mathbf{x}_{q_{j}}\right\rangle$
$\displaystyle=\left\langle\mathbf{a}_{i},\mathbf{a}_{j}\right\rangle\lim_{T\rightarrow\infty}\frac{1}{T}\sum_{q=0}^{T}\left\langle\mathbf{x}_{q},\mathbf{x}_{q}\right\rangle$
$\displaystyle=\left\langle\mathbf{a}_{i},\mathbf{a}_{j}\right\rangle\times\mathbb{E}\left[x^{2}(t)\right]$
$\displaystyle=\left\langle\mathbf{a}_{i},\mathbf{a}_{j}\right\rangle$ (43)

Now we write $\mathbf{a}_{i}$ as a function of the eigenvalues of
$\mathbf{W}$. Using the eigenvalue decomposition of the weight matrix
$\mathbf{W}$, we rewrite the state of the neuron as follows:

$\displaystyle r_{i}(t)$
$\displaystyle=\sum_{k=0}^{N_{R}}\sum_{n=1}^{N_{R}}\lambda_{n}^{k}\left\langle
v_{n}^{-1},\mathbf{w}^{\mathrm{in}}\right\rangle\left(v_{n}\right)_{i}x(t-k)$
(44)

where $\mathbf{v}_{n},\mathbf{v}_{n}^{-1}\in\mathbf{V}$ are, respectively, the
left and right eigenvectors of $\mathbf{W}$, and
$\lambda_{n}^{k}\in\mathbf{\lambda}$ belongs to the diagonal matrix containing
the eigenvalues of $\mathbf{W}$;
$\mathbf{a}_{i}=\left[a_{i,0},a_{i,1},\ldots\right]$ represents the
coefficients that the previous inputs $\mathbf{x}_{t}=[x(t),x(t-1),\ldots]$
have on $r_{i}(t)$. $\blacksquare$

Theorem 1:  If the memory capacity of the HRSNN and MRSNN networks are denoted
by $\mathcal{C}_{H}$ and $\mathcal{C}_{M}$ respectively, then,
$\displaystyle\mathcal{C}_{H}\geq\mathcal{C}_{M}$, where the heterogeneity in
the neuronal parameters $\mathcal{H}$ varies inversely to the correlation
among the neuronal states measured as
$\sum_{n=1}^{N_{\mathcal{R}}}\sum_{m=1}^{N_{\mathcal{R}}}\operatorname{Cov}^{2}\left(x_{n}(t),x_{m}(t)\right)$
which in turn varies inversely with $\mathcal{C}$.

Proof:  As shown by Aceituno et al.Aceituno et al. (2020), the memory capacity
increases when the variance along the projections of the input into the
recurrent layer state has higher heterogeneity. This can be expressed in terms
of the state space of the recurrent layer. Now, we aim to project the inputs
into orthogonal directions of the network state space. Thus, we model the
system as

$\mathbf{r}(t)=\sum_{\tau=1}^{\infty}\mathbf{a}_{\tau}x(t-\tau)+\varepsilon(t)$
(45)

where the vectors $\mathbf{a}_{\tau}\in\mathbb{R}^{N}$ correspond to the
linearly extractable effect of $x(t-\tau)$ onto $\mathbf{r}(t)$ and
$\varepsilon(t)$ is the nonlinear contribution of all the inputs onto the
state of $\mathbf{r}(t)$.

Since our goal is to have a variance as homogeneous as possible along with the
directions of $a_{\tau}$, we need a variance that is as homogeneous along with
orthogonal directions, where the vectors $\mathbf{a}_{\tau}\in\mathbb{R}^{N}$
correspond to the linearly extractable effect of the input variable $x(t)$
onto the states of the neurons ($\mathbf{r}(t)$). Since the eigenvectors of
$\mathbf{\Sigma}$ preserve orthogonality across the covariance matrix
$\mathbf{\Sigma}$, the new variances are given by the eigenvalues of the
covariance matrix, $\lambda_{n}(\mathbf{\mathbf{\Sigma}})$. Thus, we work on
the distribution of the eigenvalues of the covariance matrix. Specifically, we
want to show that increasing the heterogeneity in the neuronal membrane time
constants decreases the correlation between the neuron states, which decreases
the variance of the neuronal states of the eigenvalues, which would increase
the memory capacity $\mathcal{C}$. We quantify the heterogeneity using the
mean with respect to the square root of the raw variance of the eigenvalues of
the covariance matrix given by

$\mathcal{J}=\frac{\sum_{n=1}^{N_{\mathcal{R}}}\mathbf{\lambda}_{n}^{2}(\mathbf{\Sigma})}{\left(\sum_{n=1}^{N_{\mathcal{R}}}\mathbf{\lambda}_{n}(\mathbf{\mathbf{\Sigma}})\right)^{2}}$
(46)

where $\lambda_{n}(\mathbf{\Sigma})$ is the $n$th eigenvalue of
$\mathbf{\Sigma}$. To get an intuition of how this metric reflects the
heterogeneity in the neuronal parameters, consider the case of two eigenvalues
$\lambda_{1},\lambda_{2}$; when $\lambda_{1}=\lambda_{2}$-very homogeneous $-$
then $\mathcal{J}=\frac{1}{2}$, but when $\lambda_{1}>0,\lambda_{2}=0-$
heterogeneity is more and hence, $\mathcal{J}=1$. The membrane time constant
is given by the product of the membrane resistance $R_{m}$ and membrane
capacitance $C_{m}$, such that $\tau_{m}=R_{m}C_{m}$. $R_{m}$ is the inverse
of the permeability; the higher the permeability, the lower the resistance,
and vice versa. Thus, the lower the time constant, the faster or more rapidly
a membrane will respond to a stimulus. The effects of the time constant on
propagation velocity will become clear below. Hence, variability in the
membrane time constants will lead to variability in the propagation velocity
of action potentials.

Now,

$\left(\sum_{n=1}^{N_{\mathcal{R}}}\mathbf{\lambda}_{n}(\mathbf{\Sigma})\right)^{2}=(\operatorname{tr}[\mathbf{\Sigma}])^{2}=(\sum_{n=1}^{N_{\mathcal{R}}}\operatorname{Var}\left(r_{n}(t)\right))^{2}$
(47)

which is constant by the assumption that the probability distributions of the
neuron activities are fixed. Hence we can focus on the value of
$\sum_{n=1}^{N_{\mathcal{R}}}\mathbf{\lambda}_{n}^{2}(\mathbf{\Sigma})$ which
is true since

$\mathbf{\Sigma}^{k}\mathbf{e}_{n}(\mathbf{\Sigma})=\mathbf{\lambda}_{n}(\mathbf{\Sigma})\mathbf{\Sigma}^{k-1}\mathbf{e}_{n}(\mathbf{\Sigma})=\mathbf{\lambda}_{n}^{k}(\mathbf{\Sigma})\mathbf{e}_{n}(\mathbf{\Sigma})\Rightarrow\sum_{n=1}^{N_{\mathcal{R}}}\mathbf{\lambda}_{n}^{2}(\mathbf{\Sigma})=\operatorname{tr}\left[\mathbf{\Sigma}^{2}\right]$
(48)

where $\mathbf{e}_{n}(\mathbf{\Sigma})$ and
$\mathbf{\lambda}_{n}(\mathbf{\Sigma})$ are, resp. the $n$th eigenvector and
eigenvalue of $\mathbf{\Sigma}$. Hence, we can compute this by decomposing the
square of the covariance matrix as follows:

$\sum_{n=1}^{N_{\mathcal{R}}}\mathbf{\lambda}_{n}^{2}(\mathbf{\Sigma})=\sum_{n=1}^{N_{\mathcal{R}}}\sum_{m=1}^{N_{\mathcal{R}}}\mathbf{\Sigma}_{nm}\mathbf{\Sigma}_{mn}=\sum_{n=1}^{N_{\mathcal{R}}}\sum_{m=1}^{N_{\mathcal{R}}}\operatorname{Cov}^{2}\left(x_{n}(t),x_{m}(t)\right)$
(49)

where $\mathbf{\Sigma}_{ij}$ are the factor matrices obtained using Cholesky
decomposition of $\mathbf{\Sigma}$. Thus,
$\sum_{n=1}^{N_{\mathcal{R}}}\mathbf{\lambda}_{n}^{2}(\mathbf{\Sigma})$
increases as the neurons become more correlated; hence heterogeneity
decreases.

Thus, from Eqs. 46, 49 we can write the heterogeneity as inversely
proportional to
$\sum_{n=1}^{N_{\mathcal{R}}}\operatorname{Cov}^{2}\left(x_{n}(t),x_{m}(t)\right)$.
We see that increasing the correlations between neuronal states decreases the
heterogeneity of the eigenvalues, which would reduce the memory capacity of
the model. We show that the determinant of the covariance between neuronal
parameters bounds the heterogeneity. Thus, as $\mathcal{H}$ increases
$\rightarrow$ covariance decreases $\rightarrow$ neurons become less
correlated. Aceituno et al.Aceituno et al. (2020) proved that the neuronal
state correlation is inversely related to the memory capacity of the network.
Hence, we claim that as $\mathcal{H}$ increases, the memory capacity
$\mathcal{C}$ also increases. Hence, for HRSNN, with $\mathcal{H}>0$,
$\displaystyle\mathcal{C}_{H}\geq\mathcal{C}_{M}$. $\blacksquare$

####  B.5 Spiking Efficiency

In this section, we model the spiking activity using a point process called
the multivariate Point process model. A point process is a collection of
random points on some underlying mathematical space, such as the real line,
the Cartesian plane, or more abstract spaces.

The notion of using point process models, especially the interactive Hawkes
processes, to model the spiking dynamics of LIF network dynamics has been
studied in the literature previously (Löcherbach, 2017; Galves & Löcherbach,
2016; Mascart, 2021; Pfaffelhuber et al., 2022). We leverage these results to
prove that heterogeneity in the synaptic dynamics can help reduce the spike
count, as already discussed in the paper. We highlighted the key assumptions
used in deriving the results in the Suppl. Sec. C. We apologize if there is
still confusion, and we will add more in-depth discussion in the final
manuscript as discussed below. In their paper, Locherbach et al.Löcherbach
(2017) survey some aspects of the study of Hawkes processes in high dimensions
to model biological neural systems and study their long-term behavior. Galves
et al.Galves & Löcherbach (2016) provided an overview of point processes used
as stochastic models for interacting neurons in discrete and continuous time.
Similarly, Hawkes processes have met a recent interest in the mathematical
neuroscience literature for their ability to model the dependence of a
neuron’s activity in the network’s history (Mascart, 2021; Pfaffelhuber et
al., 2022; Galves & Löcherbach, 2016; Gerhard et al., 2017; Zhou et al., 2020;
Duval et al., 2022). Other works have also used a nonlinear interactive Hawkes
process to model spiking neural networks with excitatory and inhibitory
neurons (Chevallier et al., 2015; Chornoboy et al., 1988; Hansen et al., 2015;
Reynaud-Bouret et al., 2014). Drawing from these works, we use a microscopic
model describing a large network of interacting neurons that can generate
oscillations in a macroscopic frame. In the model, the activity of each neuron
is represented by a point process indicating the successive times at which the
neuron emits a spike, where each realization of this point process is the
spike train. We take the spiking intensity of a neuron as the probability of
emitting a spike during the next instant, depending on the history of the
neuron and the activity of other neurons in the network. The neurons interact
through their synapses. This means that a spike of a pre-synaptic neuron leads
to an increase of the membrane potential of the post-synaptic neuron if the
synapse is excitatory or a decrease if the synapse is inhibitory, possibly
after some delay, like the process of synaptic integration. The neuron fires a
spike when the membrane potential reaches a certain upper threshold. Thus,
excitatory inputs from the neurons in the network increase the firing
intensity, and inhibitory inputs decrease it. Hawkes processes provide good
models of this synaptic integration phenomenon by the structure of their
intensity processes. This paper uses a general class of mean-field interacting
Hawkes processes, modeling the reciprocal interactions between a population of
excitatory neurons and a population of inhibitory neurons.

Let us consider a subsection of the HRSNN network as shown in Fig. 11 denoted
by $N_{x}$. We use the multivariate Point process model to create a
probabilistic model that relates the inner structure of the sub-network and
its spiking activity. In this model, each neuron $i$ has a background spiking
intensity $\nu_{i}$ caused by neurons outside the network. We know that when a
neuron spikes, it impacts its spiking activity and the spiking activity of its
output neurons. The impact of a neuron $j$ on neuron $i$ is modeled by a real
function $h_{j\rightarrow i}(t)$. This impact can be excitatory or inhibitory
depending on whether the pre-synaptic neuron is excitatory or inhibitory, as
shown in Fig. 11. While the spikes from excitatory neurons try to excite
another spike, spikes originating from inhibitory neurons try to inhibit the
spiking of the cascading neuron.

A Hawkes process is a point process in which each point is commonly associated
with event occurrences in time, where every event time impacts the probability
that other events will take place subsequently. These processes are
characterized by the conditional intensity function, seen as an instantaneous
measure of the probability of event occurrences. A Hawkes process is a point
process in which each point is commonly associated with event occurrences. In
this past-dependent model, every event time impacts the probability that other
events take place subsequently. These processes are characterized by the
conditional intensity function, seen as an instantaneous measure of the
probability of event occurrences. Although the self-exciting Hawkes process
remains widely studied, there has been a growing interest in modeling the
opposite effect, known as inhibition, in which the apparition of certain
events lowers the probability of observing an event. In practice, this amounts
to considering negative kernel functions. To maintain the positivity of the
intensity function, a non-linear operator is added to the expression, which in
turn entails the loss of the cluster representation. This model is known as
the non-linear Hawkes process, where the existence of such processes was
proved via construction using bi-dimensional marked Poisson processes. The
general Hawkes framework can be written as:

$\lambda_{t}^{i}=\Phi_{i}\left(\sum_{j\in\mathcal{S}_{i,E}}\int_{0}^{t}h_{j\rightarrow
i}(t-u)\mathrm{d}Z_{u}^{j}\right),$ (50)

where $\lambda_{t}^{i}$ is the intensity of neuron $i,\Phi_{i}$ a positive
function, $Z_{j,t}$ is the counting process associated with neuron
$j,h_{j\rightarrow i}(t)$ is the synaptic kernel associated with the synapse