Multi-scale Evolutionary Neural Architecture Search for Deep Spiking
  Neural Networks
Spiking Neural Networks (SNNs) have received considerable attention not only
for their superiority in energy efficiency with discrete signal processing but
also for their natural suitability to integrate multi-scale biological
plasticity. However, most SNNs directly adopt the structure of the
well-established Deep Neural Networks (DNNs), and rarely automatically design
Neural Architecture Search (NAS) for SNNs. The neural motifs topology, modular
regional structure and global cross-brain region connection of the human brain
are the product of natural evolution and can serve as a perfect reference for
designing brain-inspired SNN architecture. In this paper, we propose a
Multi-Scale Evolutionary Neural Architecture Search (MSE-NAS) for SNN,
simultaneously considering micro-, meso- and macro-scale brain topologies as
the evolutionary search space. MSE-NAS evolves individual neuron operation,
self-organized integration of multiple circuit motifs, and global connectivity
across motifs through a brain-inspired indirect evaluation function,
Representational Dissimilarity Matrices (RDMs). This training-free fitness
function could greatly reduce computational consumption and NAS's time, and its
task-independent property enables the searched SNNs to exhibit excellent
transferability on multiple datasets. Furthermore, MSE-NAS show robustness
against the training method and noise. Extensive experiments demonstrate that
the proposed algorithm achieves state-of-the-art (SOTA) performance with
shorter simulation steps on static datasets (CIFAR10, CIFAR100) and
neuromorphic datasets (CIFAR10-DVS and DVS128-Gesture). The thorough analysis
also illustrates the significant performance improvement and consistent
bio-interpretability deriving from the topological evolution at different
scales and the RDMs fitness function.
Wenxuan Pan, Feifei Zhao, Guobin Shen, Yi Zeng
# Emergence of Brain-inspired Small-world Spiking Neural Network through Neuroevolution

###### Abstract

Human brain is the product of evolution during hundreds over millions of years and can engage in multiple advanced cognitive functions with low energy consumption. Brain-inspired artificial intelligence serves as a computational continuation of this natural evolutionary process, is imperative to take inspiration from the evolutionary mechanisms of brain structure and function. Studies suggest that the human brain's high efficiency and low energy consumption may be closely related to its small-world topology and critical dynamics. However, existing efforts on the performance-oriented structural evolution of spiking neural networks (SNNs) are time-consuming and ignore the core structural properties of the brain. In this paper, we propose a multi-objective Evolutionary Liquid State Machine (ELSM) with the combination of small-world coefficient and criticality as evolution goals and simultaneously integrate the topological properties of spiking neural networks from static and dynamic perspectives to guide the emergence of brain-inspired efficient structures. Extensive experiments show a consistent and comparable performance of the proposed model compared to LSM-based and hierarchical SNNs algorithms: it achieves 97.23%
on NMNIST, and reaches the state-of-art performance compared to all LSM models on MNIST and Fashion-MNIST (98.05% and 88.81%, respectively). A thorough analysis reveals the spontaneous emergence of hub nodes, short paths, long-tailed degree distributions, and numerous community structures in evolutionary models. This work evolves recurrent spiking neural networks into brain-inspired efficient structures and dynamics, providing the potential to achieve adaptive general artificial intelligence.

keywords: Spiking Neural Networks, Neuroevolution, Small-world Topologies, Critical Dynamics, Liquid State Machines +
Footnote †: journal: Journal of Computational and Graphical Statistics

## 1 Introduction

How can the human brain perform many complex advanced cognitive functions while running on less power than a light bulb? Its mysterious wiring rules and firing patterns have attracted much research interest. There is a degree of commonality in brain anatomy across the human species: different regions are often thought to be responsible for specific cognitive functions [1]. It is worth mentioning that the densely connected community structure and hub nodes existing in these specific regions help efficient information processing and integration in the brain [2; 3].

From a static topology perspective, researchers have demonstrated that the mammalian cortical (including the human brain) is a complex network whose topological properties are neither random nor regular, but somewhere in between [4; 5; 6; 7], with small-world properties of dense local clustering and short path length [4; 8; 9]. From a network dynamics perspective, when dealing with complex and changeable environments, the human brain can exhibit powerful and flexible adaptive processing capabilities and achieve a delicate balance between efficiency and robustness due to its well-evolved internal wiring rules. In this case, biological neural network dynamics achieve optimal computational and processing capabilities near a certain point called a critical state where networks oscillate between order and disorder, synchronous and asynchronous [10; 11; 12],
and can reach optimal value of information transmission [13].

Efficient transmission topology and optimal dynamics enable the human brain to exhibit powerful low-energy, high-efficiency information processing capabilities. The formation of such structures is not artificially designed but evolved naturally. Existing human-crafted network structures may help improve performance, but difficult to escape from the inherent paradigm [14]. To enable the model to find the optimal network architecture adaptively, a field called Neural Architecture Search (NAS) is emerging [15; 16; 17; 18; 19; 20]. Most work on NAS follows the wave of deep learning to search deep network structures [21; 22; 23; 24; 25], however, to the best of our knowledge, no NAS algorithm takes into account the biologically economical small-world topology and criticality in brain.

Spiking neural network (SNN), as the third-generation neural network, not only simulates the discrete communication of biological neurons but also can combine multiple biological plasticity learning rules, which is more in line with the information processing mechanism in the biological brain [26]. In this paper, we employ a large-scale, recurrently connected SNN called Liquid State Machines (LSM) [27], a kind of reservoir network, due to its advantages of complex liquid structure, low training cost, good at processing spatiotemporal information, and more suitable for studying brain-inspired connective architecture [27; 28; 29; 30; 31; 32]. A standard LSM consists of three parts: the input information is processed by a liquid layer containing randomly fixed connections, and then abstracted by the readout neurons into the final output. The only SNN-based NAS works search for neuron operations [33; 34] and cross-layer connections [34], while the evolutionary LSM works evolve reservoir parameters, including structural parameters such as liquid density, excitatory neuron ratio, number of liquid neurons [35; 36]. [37] changes the structure of LSM by dividing a large liquid into multiple smaller liquids. These studies lack the in-depth inspiration of the unique topological characteristics of the brain, thus limiting their learning efficiency and performance.

Inspired by the small-world properties and the critical state of biological
nervous system, this paper proposes an evolutionary Liquid State Machines (ELSM) for emerging with brain-inspired small-world architecture and dynamic firing patterns. **Structurally**, the evolved liquid layer presents the network architecture characteristics of small-world, combining dense local clustering and short path length. **Dynamically**, ELSM enables the network to operate near the critical state, which is significantly more efficient and biologically plausible. The proposed multi-objective evolutionary algorithm not only brings the small-world structure and critical dynamics, but also naturally achieves higher performance and efficiency.

The main highlights of this paper can be summarized in the following three points:

* We evolve the structure of the recurrent spiking neural network to exhibit biological plausible small-world topological properties (densely local-connected hub nodes, large number of communities and long-tail degree distribution) as well as dynamic critical steady state. Brain-inspired evolutionary goals simultaneously bring about an improvement in classification accuracy.
* The proposed multi-objective evolutionary algorithm considers small-world coefficients (including shortest path length and clustering coefficients) and criticality as fitness functions, guiding the emergence of brain-inspired efficient structures.
* Our model achieves classification accuracies of 98.05%, 97.23% and 88.81% on MNIST, NMNIST and Fashion-MNIST respectively, which are comparable performance to deep SNNs. Experimental results demonstrate that adaptively evolved LSMs improve performance with biologically plausible structures and firing patterns at lower complexity. The degree distribution of the evolved network nodes exhibits the characteristics of a long-tailed distribution, similar to that found in biological brains.
## 2 Results

### The Architecture of Reservoir-based SNN

The architecture of reservoir-based SNN (LSM) is shown in Fig 1. The standard LSM model is divided into three layers: input layer, liquid layer formed of thousands of neurons sparsely connected, and readout layer. All neurons accumulate potentials according to the law shown in Eq. 2, and transmit information through spikes. The weight between the readout layer and the liquid layer is optimized by backpropagation algorithm [38], while the weights in the liquid layer are randomly fixed.

In this paper, we use the leaky integrate-and-fire (LIF) neuron as the basic unit of signal transmission, and the formula for its membrane potential update over time is:

\[\delta=\frac{I(t)-V_{\mathrm{m}}(t)}{\tau} \tag{1}\]

\[V_{\mathrm{m}}(t+1)=\left(V_{\mathrm{m}}(t)+\delta\right)\left(1-S(t)\right)+ V_{r}S(t) \tag{2}\]

Figure 1: The architecture of LSM. In the traditional definition of a reservoir, randomly connected spiking neurons receive time-varying signals from external inputs and other nodes simultaneously. The recursive connectivity enables input signals to be converted to liquid layer dynamics, which are then abstracted by the readout layer.
\[S(t)=\left\{\begin{array}{l}1,V_{\text{m}}(t)\geq V_{th}\\ 0,V_{\text{m}}(t)<V_{th}\end{array}\right. \tag{3}\]

\(V_{m}(t+1)\) and \(V_{m}(t)\) are the membrane potential at time \(t+1\) and \(t\), respectively. As shown in Eq. 1, \(\delta\) is determined by the membrane potential \(V_{m}(t)\), the magnitude of the current \(I(t)\) and the membrane potential time constant \(\tau\). When the membrane potential reaches the threshold \(V_{th}\), the membrane potential is reset to \(V_{r}\) at the same time as the spike is delivered (indicated by \(S(t)\) as Eq. 3). According to the membrane potential \(V_{m}(t)\) and \(S(t)\) at time \(t\), the update law of the membrane potential at the next time time is shown in Eq. 2. The LIF neuron model and learning rules of the proposed evolutionary LSM are based on BrainCog framework [39].

```
0: Population \(P(0)=\{C_{1},C_{2},...,C_{N_{c}}\}\);
0: Evolved individual \(C_{opt}\) ; for\(g=0\) to \(G_{th}\)do if\(g=G_{th}-1\)then \(Accuracy\) = Train (\(P(g)\), 100) \(C_{opt}\) = Max (\(P(g)\), \(Accuracy\)) return\(C_{opt}\) endif \(obj[g,0]\) = SmallWorld (\(P(g)\)) \(obj[g,1]\) = Criticality (\(P(g)\), data) \(Q(g)\)=CrossoverAndMutate (Select(\(P(g)\),\(obj\),\(N_{off}\))) \(P(g+1)\) = Merge (\(P(g)\), \(Q(g)\)) \(P(g+1)\) = Select (\(P(g+1)\),\(obj\), \(N_{c}\)) endfor
```

**Algorithm 1** The neuroevolution process of ELSM.

### Neuroevolution Algorithm

Liquid layer connectivity formed by random initialization in the reservoir will be evolved to emerge more brain-inspired structures and dynamics. The
whole neuroevolution process is presented as Algorithm 1 and Fig. 2.

#### 2.2.1 Fitness Function

**Identifying small-world topologies.**[40] describes small-world networks exhibit two properties: highly clustered and short path length. Local short paths between most nodes with hubs induce highly connected sub-networks and a few long-distance connections, enabling efficient information transmission in the brain. We refer to the quantification method of [41] to quantify the small-world characteristics, which is called small-world coefficient, and the calculation is as

Figure 2: **The neuroevolution process of ELSM.****a.** Schematic diagram of initialization. After determining the number of neurons in the liquid layer \(N\) and the initial liquid density \(\rho_{init}\), phenotypes (LSM liquid layers) are converted to genotypes (connection patterns \(C\), sizes of which are \(N*N\)). **b.** The process of crossover, mutation and selection. Initialized individuals make up \(P(g)\) in generation \(g\), and \(Q(g)\) is obtained by performing k-point crossover and bitflip mutation operations on the mating pool generated by \(P(g)\). Merge \(Q(g)\) and \(P(g)\), apply elitism method and non-dominated sorting strategy in NSGA-II to select the next generation group \(P(g+1)\), until the iteration reaches \(G_{th}\) and the evolution ends. **c.** Brain-inspired evolution objective functions. Statically, evolve the small-world coefficient of \(C\). Dynamically, a small amount of training data is input to obtain firing patterns of the liquid layer within \(T\) time period, and the branching ratio \(\mu(t)\) is calculated as the quantification of criticality.
follows:

\[\lambda=\frac{H}{L} \tag{4}\]

The clustering coefficient and short path length between nodes is measured by \(H\) and \(L\), respectively. The clustering computing of a single node is as Eq. 5, where \(o_{i}\) is the degree of neuron \(i\). If two edges that both pass through neuron \(i\) are called a pair, \(e_{i}\) is the number of all pairs of edges with \(i\) as the intermediate node. The total clustering coefficient is the average of total neurons (assuming there are \(N\) neurons in the liquid layer).

\[h_{i} =\frac{2e_{i}}{o_{i}\left(o_{i}-1\right)} \tag{5}\] \[H =\frac{\sum_{i}^{N}h_{i}}{N} \tag{6}\]

The shortest path length of the network is calculated as Eq. 7, where \(d_{st}\) represents the shortest path length between \(s\) neuron and \(t\) neuron (\(d_{st}=0\) if the path does not exist). \(V\) is the liquid neuron set.

\[L=\sum_{s,t\in V}\frac{d_{st}}{N(N-1)} \tag{7}\]

**Identifying criticality.** A commonly used concept to measure the critical state of the nervous system is called branching ratio, derived from the branching process theory [42], reflecting the spatiotemporal cascade activity of the cerebral cortex in homeostasis. The local branching ratio \(\mu_{i}(t)\) describes the tendency of a neuron \(i\) to be more or less active as spikes transmitted within the liquid layer, defined as follow [43] at time \(t\).

\[\mu_{i}(t)=\frac{\sum_{j}^{N}\sum_{l=t+\phi+1}^{t+\phi+\Delta}m_{j}(l)c_{ij}}{ \sum_{j}^{N}\sum_{l=t-\phi-\Delta}^{t-\phi-1}m_{j}(l)c_{ji}} \tag{8}\]

\[\mu(t)=\frac{\sum_{i}^{N}m_{i}(t)\mu_{i}(t)}{\sum_{i}^{N}m_{i}(t)} \tag{9}\]
\[\mu=\frac{\sum_{t}^{T}\mu(t)}{T} \tag{10}\]

Where \(t=1,2,3,...,T\). \(c_{ij}\) indicates whether there is a synapse connecting presynaptic neuron \(i\) and postsynaptic neuron \(j\). \(m_{i}(t)\) represents the firing situation of neuron \(i\) at time \(t\) (\(m_{i}(t)\) has only two values of 0 or 1, which represents not firing or firing respectively). Therefore, Eq. 10 is the ratio of the sum of the postsynaptic neuron spikes to the sum of the presynaptic neuron spikes during the simulated time \(T\). Studies have shown that the closer the value of \(\mu\) is to 1, the closer the dynamics of the network are to the critical state [43; 44; 45]. Therefore, the quantitative criticality is calculated as:

\[\mu=|\mu-1| \tag{11}\]

#### 2.2.2 Initialization

In a population of \(N_{c}\) chromosomes to be initialized, each chromosome representing the liquid connection pattern \(C\) of a \(N\) liquid neurons reservoir, as shown in Fig. 2a. We use a binary encoding method, each chromosome locus \(c_{ij}\) has two alleles of 0 or 1, indicating whether there is a synapse \(i\) connecting the presynaptic neuron and a post-synaptic neuron \(j\) (\(0<i,j<N\)). At the beginning of neuroevolution, each chromosome's connectivity is limited to be sparse, and has only \(N*N*\rho_{init}\) synaptic connections inside (the initial liquid density is is recorded as \(\rho_{init}\)).

We first generate a random matrix \(R\), the values of all elements \(r_{ij}\) in \(R\) are between 0 to 1. The Boolean matrix \(C\) is obtained by calculating the result of \(R<\rho_{init}\), as Eq. 12:

\[c_{ij}=\begin{cases}1,&\text{ if }r_{ij}\leq\rho_{init}\\ 0,&\text{ if }r_{ij}>\rho_{init}\end{cases} \tag{12}\]
#### 2.2.3 Evaluation

Since there is more than one fitness function, we will fully consider the physical topology characteristics and dynamic changes of the network when selecting excellent chromosomes as the parents of the next generation. According to the fitness function proposed in the section 2.2.1, the larger the small-world coefficient \(\lambda\), the more the model structure can reflect the characteristics of the small-world. The smaller \(\mu\) is, the closer the model dynamics are to the critical state. Therefore, this multi-objective optimization problem (MOP) according to Eq. 4 and Eq. 10 can be described as \(\mathbf{F}:\mathbf{\Omega}\rightarrow\mathbb{R}\):

\[\underset{C\in\Omega}{\mathrm{argmin}}\mathbf{F}(C)=\{f_{1}(C),f_{2}(C)\} \tag{13}\]

\[s.t.\rho_{1}\leq\rho(C)\leq\rho_{2}\]

\(\rho(C)\) is the density of the liquid layer, defined as the ratio of the number of liquid layer connections to \(N*N\). In order to keep the liquid density \(\rho(C)\) stable during the evolution process, the range of which is limited to \(\rho_{1}\) and \(\rho_{2}\). The evolution content of the two objectives is shown in Fig. 2c. The first optimization goal \(f_{1}(C)\) is to maximize the small-world coefficient:

\[f_{1}(C)=\min(-\lambda(C)) \tag{14}\]

where \(\lambda\) is calculated as Eq. 4, measuring the static LSM topology properties.

The second optimization objective \(f_{2}(C)\) is to minimize the criticality coefficient, which is formulated by:

\[f_{2}(C)=\min\mu(C) \tag{15}\]

where \(\mu\) is calculated as Eq. 11, measuring criticality of LSM dynamics.

#### 2.2.4 Selection

The elitism approach and nondominated sorting strategy of the NSGA-II algorithm [46] is used here to generate mating pools of size \(N_{offs}\) and next
generation individuals.

#### 2.2.5 Crossover and Mutation

**Crossover.** Assuming that chromosome \(C_{1}\) crosses with \(C_{2}\) to generate \(C_{3}\) and \(C_{4}\), we apply the k-point crossover operator rule [47] to select \(k\) genes as crossover points:

\[c_{a_{1}b_{1}},c_{a_{2}b_{2}},...,c_{a_{k}b_{k}} \tag{16}\]

where \(0<a_{k},b_{k}<N\). According to the selected crossover points, each gene can be divided into \(k+1\) segments:

\[\frac{(c_{0,0},c_{a_{1}b_{1}})}{S_{1}},\frac{(c_{a_{1}b_{1}+1},c_{a_{2}b_{2}}) }{S_{2}},...,\frac{(c_{a_{k}b_{k}+1},c_{NN})}{S_{k+1}} \tag{17}\]

We get the matrix \(E\) and \(D\) where \(E_{s_{1}}\), \(E_{s_{3}}\),..., \(D_{s_{2}}\), \(D_{s_{4}}\),..., are set to 1. The remaining elements in \(E\) and \(D\) are set to 0. Therefore:

\[C_{3} =C_{1}*E+C_{2}*D \tag{18}\] \[C_{4} =C_{1}*D+C_{2}*E \tag{19}\]

**Mutation.** Assuming that \(C_{1}\) is mutated into \(C_{5}\), we perform flip bit mutation [48] on \(n_{m}\) genes of \(C_{1}\). Let the mutation probability be \(m_{rate}\), generate a random number \(m_{rand}\) for each offspring after crossover, if \(m_{rand}<m_{rate}\), the mutation is accepted. Each mutation is to select a chromosome locus \(c_{ij}\) in \(C_{1}\) for inversion:

\[c_{ij}=\neg c_{ij} \tag{20}\]

Parents \(P(g)\) and population formed by crossover and mutation \(Q(g)\) and are merged, and the selection operator is applied to generate the next generation \(P(g+1)\). The whole crossover mutation and selection process is shown in Fig. 2b.

#### 2.2.6 Next Generation

The above neuroevolution process of evaluation, selection, crossover, and mutation is repeated for \(G_{th}\) generations as shown in Fig. 2b. Finally, all individuals in the last generation are trained for 100 epochs, and the individual with the highest classification accuracy \(C_{opt}\) is selected as the result of evolution.
### Datasets and Parameter Settings

Our model is validated on MNIST [49], NMNIST [50] and Fashion-MNIST [51] datasets to prove the effectiveness of the algorithm.

**MNIST.** The handwritten dataset MNIST is one of the classic machine learning datasets, which contains 70,000 grayscale images of handwritten digits 0-9, of which 60,000 examples are used for training, and the others are used for testing. The size of each image is \(28*28\) pixels.

**Fashion-MNIST.** The Fashion-MNIST dataset consists of 60,000 training samples and 10,000 testing samples, with a total of ten categories. Each sample is a \(28*28\) grayscale image.

**NMNIST.** NMNIST is a neuromorphic version of MNIST converted to MNIST images by an actuated pan-tilt camera platform. After 300 ms of signal acquisition, 60,000 training images and 10,000 test images are generated, each with a size of \(34*34*2\) (2 is the number of channels, 34 is the size of the image after offset). A preprocessing ensemble method [52] is adopted to convert the event stream into frame stream, which is then fed into the model for classification. The two channels are combined into one by summing in our experiments.

All parameter settings of the evolutionary algorithm in the experiment are shown in Table 1. All network weights, including weights from input to liquid layer, weights from liquid layer to output, and weights inside the liquid layer (the value of the weight rather than the connection mode, which is obtained by the neuroevolution) are initialized randomly. The batch size of all datasets is set to 100. The Adam optimizer [53] is adpoted to optimize the weight of the readout layer, whose learning rate's decay rate (every 50 epochs) is set to 0.0001.

In order to demonstrate the effectiveness of neuroevolution, we established four comparison models: 1) baseline LSM model with randomly generated liquid layers (marked as RLSM), 2) evolved small-world topologies LSM model after 1000 generations (marked as ESLSM), 3) evolved criticality LSM model after 1000 generations (marked as ECLSM), 4) evolved multi-object LSM model after 1000 generations (the proposed model, marked as ELSM) for the ablation studies
of fitness functions on MNIST, NMNIST and Fashion-MNIST. The evolved \(C_{opt}\) is trained as a liquid layer connection mode for 5000 epochs, and multiple random seeds are replaced to perform multiple evolutions to obtain the results in Table 3 and Fig. 3. Table 2 shows the best result among multiple evolutions.

During training, except for the liquid layer (the best individual), the other settings of the comparison model are the same, including the weight of the liquid layer, the connection between the input and the liquid layer, and the connection between the liquid layer and the readout layer to ensure the fairness of the comparison experiment.
NAS-LSM [35] on NMNIST by 4.73% as shown in Table 2. On MNIST, ELSM surpasses other evolved-architecture LSM, LSM-SHADE [36] and Multi-liquid LSM [37] by 3.55% and 2.55% respectively. Overall, ELSM exhibits better performance than all other LSM models (except NMNIST) and comparable results to deep SNN on all datasets.

\begin{table}
\begin{tabular}{c c c c c} \hline \hline Dataset & Model & Structure & Layers & Accuracy (\%) \\ \hline \multirow{8}{*}{MNIST} & Unsupervised-SNN [54] & Hierarchical SNN & 2 & 95 \\  & LIF-BA [55] & Hierarchical SNN & 3 & 97.09 \\  & Temporal SNN [56] & Hierarchical SNN & 2 & 97.2 \\  & STIDi-BP [57] & Hierarchical SNN & 2 & 97.4 \\  & SN [58] & Hierarchical SNN & 3 & 97.93 \\  & BPSNN [59] & Hierarchical SNN & 3 & 98.88 \\ \cline{2-5}  & CMA-ES-LSM [60] & LSM & 2 & 92.6 \\  & LSM-SHADE [36] & LSM & 2 & 94.5 \\  & Multi-liquid LSM [37] & LSM & 2 & 95.5 \\  & NALSM [45] & LSM & 2 & 97.61 \\  & **ELSM** & **LSM** & **2** & **98.05** \\ \hline \multirow{8}{*}{NMNIST} & DECOLLE [61] & Hierarchical SNN & 2 & 96 \\  & AER-SNN [62] & Hierarchical SNN & 2 & 96.3 \\  & BPSNN [59] & Hierarchical SNN & 3 & 98.74 \\  & SLAYER [63] & Hierarchical SNN & 3 & 98.89 \\ \cline{2-5}  & Ionic LSM [64] & LSM & 2 & 91.48 \\  & NAS-LSM [35] & LSM & 2 & 92.5 \\  & **ELSM** & **LSM** & **2** & **97.23** \\  & NALSM [45] & LSM & 2 & 97.51 \\ \hline \multirow{8}{*}{Fashion-MNIST} & SL-SNN [65] & Hierarchical SNN & 3 & 85.3 \\  & Unsupervised-SNN [54] & Hierarchical SNN & 2 & 85.31 \\ \cline{1-1}  & BS4NN [66] & Hierarchical SNN & 2 & 87.3 \\ \cline{1-1} \cline{2-5}  & NALSM [45] & LSM & 2 & 85.84 \\ \cline{1-1}  & **ELSM** & **LSM** & **2** & **88.81** \\ \hline \hline \end{tabular}
\end{table}
Table 2: Comparative performance of different LSM and SNN models on MNIST,NMNIST, and Fashion-MNIST datasets.
The multi-objective neuroevolution process we designed does not directly use the classification accuracy as a criterion for evaluating fitness, but guides the evolutionary algorithm towards an efficient direction from the perspective of physical topology and network dynamics in the brain. Surprisingly, as the evolution proceeds, individuals not only have the brain-inspired properties of small-world and critical state, but also have significantly improved classification accuracy, especially compared with other LSM models, and achieved performance comparable to other hierarchical models with lower energy consumption (minimum number of layers) as shown in Table 2.

### Ablation Study

To explore the effect of different evolution goals on individual performance, we conduct ablation experiments on each dataset using four models: RLSM, ECLSM, ESLSM and ELSM. The specific accuracy of the four models on each dataset is shown in Table 3. Randomly generated LSMs have low accuracy and large variance. ECLSM performs 97.51 \(\pm\) 0.06%, 96.7 \(\pm\) 0.24% and 88.54 \(\pm\) 0.1% on MNIST, NMNIST and Fashion-MNIST, which is significantly better than all RLSMs by 1.62%, 6.6% and 3.96%. ESLSM achieves 97.88 \(\pm\) 0.12%, 96.65 \(\pm\) 0.23% and 88.24 \(\pm\) 0.13% on MNIST, NMNIST and FMNIST, outperforming RLSMs by 1.99%, 6.55% and 3.66% respectively. Evolutionary models perform better than random LSM, among which ELSM with small-world and criticality as evolution goals performs the best on each dataset. It can be seen from Table 3 that ELSM has better accuracy and smaller tolerance than ESLSM and ECLSM, reaching 98.02 \(\pm\) 0.03%, 97 \(\pm\) 0.23% and 88.78 \(\pm\) 0.04%performance on MNIST, NMNIST and Fashion-MNIST respectively.

The comparison of ESLSM, ECLSM and ELSM on different datasets is shown in Fig. 3. In terms of criticality, different ECLSM and ELSM have been evolved for multiple times (different random seeds) for different datasets. The small-world property is independent of data, so the result of once evolution can be used for all datasets while the selected individual \(C_{opt}\) may be different for each dataset, since it is based on the training accuracy after 100 epochs.
Figure 3: Comparison of ESLSM, ECLSM and ELSM on different datasets. **a-c.** Results on evolving criticality on MNIST, NMNIST and Fashion-MNIST. The horizontal axis represents the distance between the criticality of the individual and 1, as shown in Eq. 11. **d-f.** Results on evolving small-world properties on MNIST, NMNIST and Fashion-MNIST. The horizontal axis represents the small-world coefficient of the individual as shown in Eq. 4. The green dot represents the result of 5000 epochs training of \(C_{opt}\) (the one with the highest classification accuracy after training for 100 epochs) selected in each 100 generations, and is fitted by a polynomial (blue line). The green marks with variance indicate the final single-object evolution results ECLSM and ESLSM, and the red mark indicates the multi-object evolution result ELSM.
For Fig. 3a-c, the smaller the value on x-axis, the stronger the criticality. For Fig. 3d-f, the larger the value on x-axis, the more obvious the small-world characteristics. From the results of polynomial fitting (shown in blue), it can be seen that as the evolution proceeds, the fitness of \(C_{opt}\) of each generation increases continuously, and the classification accuracy also increases. It shows a certain degree of positive correlation between the indirect time-saving evolution goal of ELSM and the classification accuracy. In Fig. 3a-f, red marks are always at the top, visualizing that ELSM exhibits better performance than single-objective evolution while ensuring maximized two evolutionary objectives.

## 3 Discussion

Current work on the performance-oriented evolution of SNN architectures is often time-consuming and fails to fully incorporate the topological properties observed in biological brains. Here, we proposes an evolutionary recurrent SNN model ELSM, which indirectly takes the more brain-inspired static small-world topological characteristics and dynamic criticality as evolution goals, replacing the time-consuming performance-oriented fitness function. ELSM achieves clas

\begin{table}
\begin{tabular}{c c c} \hline \hline Dataset & Model & Accuracy (\%) \\ \hline \multirow{4}{*}{MNIST} & RLSM & \(95.89\pm 0.97\) \\  & ESLSM & \(97.88\pm 0.12\) \\  & ECLSM & \(97.51\pm 0.06\) \\  & ELSM & \(98.02\pm 0.03\) \\ \hline \multirow{4}{*}{MNIST} & RLSM & \(90.1\pm 29.94\) \\  & ESLSM & \(96.65\pm 0.23\) \\  & ECLSM & \(96.7\pm 0.24\) \\  & ELSM & \(97\pm 0.23\) \\ \hline \multirow{4}{*}{Fashion-MNIST} & RLSM & \(84.58\pm 1.48\) \\  & ESLSM & \(88.24\pm 0.13\) \\ \cline{1-1}  & ECLSM & \(88.54\pm 0.1\) \\ \cline{1-1}  & ELSM & \(88.78\pm 0.04\) \\ \hline \hline \end{tabular}
\end{table}
Table 3: Final performance of models with different evolution goals on all datasets.
sification accuracies of 98.05%, 97.23% and 88.81% on MNIST, NMNIST and Fashion-MNIST respectively, and outperforms the best LSM models reported so far by 0.44% and 2.97% on MNIST and Fashion-MNIST, surpassing many deep SNN models and comparable to the best performance. On NMNIST, ELSM also achieved comparable results to deep SNNs and the best LSM model. The ablation experiments confirmed that the above two evolutionary goals have a certain degree of positive correlation with the classification accuracy, and the performance of the evolved model far exceeds that of the random LSM. At the same time, the multi-objective evolution model (the proposed ELSM) performs better than the single-objective model ESLSM and ECLSM.

Compared with other deep SNNs, ELSM achieves performance surpassing many deep models with the fewest layers (2 layers). We further analyze more structural brain-inspired topological features as follows:

### Hourglass Structure and Sparse Coding in Drosophila Mushroom Body

Some studies have found an hourglass-like mapping relationship in the mushroom body module of the _Drosophila_ brain: the nervous system converges from the ultra-high-dimensional signals provided by sensory cells to a small num

Figure 4: The evolvable LSM inspired by the _Drosophila_ mushroom body. Dashed box marks the hourglass structure found in the _Drosophila_ mushroom body, consistent with LSM.
ber of projection neurons (PNs), and conduct sparse encoding through a large number of Kenyon cells (KCs). Finally, lower-dimensional signals are extracted to characterize the real world. This is similar to the way that LSM processes information, as shown in Fig. 4. The dotted box marks the hourglass structure similar to that in the mushroom body of _Drosophila_, which is also the difference between LSM and hierarchical neural network.

### Emergence of Structural Properties Exist in the Brain

To study the effect of evolution on the brain-inspired topology of the model, we counted the changes in properties such as clustering coefficient, community and criticality, as shown in Table 4. Clustering coefficient is used to measure the degree of node aggregation. Communities counts the number of communities with a size of 5 in the network that can communicate through 4 common nodes.

While increasing the clustering coefficient, the shortest path length of the network is also reduced (from 0.196 to 0.195), but the total number of connections does not increase but decreases: the connection density of the random network is about 1%, while the connection density of the evolutionary network is individual 0.8%, proving that evolution does not achieve fast and efficient information transfer by increasing connections, but can reduce cost consumption (fewer connections). After evolution, due to the increase of hub nodes, the clustering coefficient of the network increases significantly, and more overlapping communities are also found, implying that the connections between sub-networks are more intricate and highly interconnected. Under limited connectivity, larger clustering coefficients, closer community connections, shorter

\begin{table}
\begin{tabular}{l c c c c} \hline \hline \multirow{2}{*}{Model} & Clustering coefficient & Communities & Criticality & Density \\  & (H) & (k=4) & (\(|\mu-1|\)) & (\(\rho\)) \\ \hline Random & 285.15 & 2 & 0.265 & 1\% \\
**Evolved** & **319.65** & **1255** & **0.096** & **0.8\%** \\ \hline \hline \end{tabular}
\end{table}
Table 4: Changes in topology properties of ELSM before and after neuroevolution.
shortest paths, and more critical connectivity patterns prove that our proposed multi-objective evolutionary algorithm can optimize LSM from both static and dynamic perspectives, the evolved network architecture is more in line with the core structural characteristics found in the human brain.

The change of the network degree distribution before and after evolution is shown in Fig. 5. Fig. 5a shows the degree distribution of random LSM, following a normal distribution, and the degrees of all nodes are concentrated between 100 and 200. After 1000 generations of evolution, the small-world properties of the network become obvious, and the degree shows an obvious long-tail distribution, as shown in Fig. 5b. A small number of hub nodes appeared: the degree of most nodes among the 8000 neurons is concentrated between 0 and 100. The greater the degree, the fewer nodes, and the maximum node degree is 1525.

Overall, evolved with more biologically plausible small-world coefficient and criticality, ELSM surpasses the best reported LSM models to date on MNIST and Fashion-MNIST, and outperforms many deep SNN models on all datasets. Through analysis, it can be found that the evolved model also has many topological structures similar to brain networks, such as hub nodes, communities, and short path lengths.

Figure 5: Degree distribution comparison. a. Degree distributions of liquid layers of random structures. b. Degree distribution of the evolved individual.
## 4 Limitations of the study

In the future, with the advancement of neuroscience, we will further in-depth explore more topological properties found in brain networks, hoping to discover more effective and energy-saving brain-inspired modular structural features on SNNs, which can be used to guide efficient evolution. In terms of application, SNN networks of various architectures (not limited to LSM or other deep SNNs) can be used to constitute multi-brain areas in the form of global self-organization and co-evolution, realize a variety of advanced cognitive functions and be applied to the exploration of transfer learning, continuous learning and other issues, not limited to classification tasks.

## 5 Resource availability

### Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Yi Zeng (yi.zeng@ia.ac.cn).

### Materials availability

This study did not generate new unique reagents.

### Data and code availability

The Python scripts of ELSM can be downloaded from the GitHub repository:

[https://github.com/BrainCog-X/Brain-Cog/tree/dev/examples/Structural_Development/](https://github.com/BrainCog-X/Brain-Cog/tree/dev/examples/Structural_Development/)

## 6 Acknowledgments

This work is supported by the National Key Research and Development Program (Grant No. 2020AAA0107800), the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDB32070100), the National Natural Science Foundation of China (Grant No. 62106261).
## 7 Author contributions

W.Pan, F.Zhao, B.Han and Y.Dong designed the study under the supervision of Y.Zeng. W.Pan, F.Zhao, B.Han Y.Zeng and Y.Dong performed the experiments and the analyses. W.Pan, F.Zhao and Y.Zeng wrote the paper.

## 8 Declaration of interests

The authors declare that they have no competing interests.
