Blended Glial Cell’s Spiking Neural Network

Spiking Neural Networks (SNNs), the third generation of artificial neural networks, have been widely employed. However, the realization of advanced artificial intelligence is challenging due to the dearth of efficient spatiotemporal information integration models. Inspired by brain neuroscientists, this paper proposes a novel spiking neural network - Blended Glial Cell’s Spiking Neural Network (BGSNN). BGSNN introduces glial cells as spatiotemporal information processing units based on neurons and synapses, and also provides four new network dynamics connection models which extend the information processing dimension, enhance the network global information integration in the spatiotemporal domain, as well as the plasticity of neurons and synapses. In this paper, a BGSNN application - Sudoku solver is designed and implemented on the “WenTian” neuromorphic prototype. On the Easybrain dataset, the BGSNN solver achieves 100% accuracy, outperforming the same structure SNN solver by 97% at the Evil difficulty level, and has faster converges speed compared with the SOTA Sudoku solver LSGA. On the kaggle dataset, the BGSNN solver achieves over 99.99% accuracy, outperforming the publicly available optimal DNN solver under this dataset by 3.82%. In addition, BGSNN exhibits good parallelism and sparsity, decreasing computation by at least 92.9% compared to serial solvers and reducing sparsity by 88% compared to the equal fully dense DNN. BGSNN improves the expression, feedback, and regulation capabilities of neural networks while maintaining the advantages of SNN parallel sparsity, making it simpler to implement advanced artificial intelligence.


I. INTRODUCTION
The biological brain is a highly intelligent system existing in nature. SNNs adopt the biological neural computing model of the biological brain [1], which is more biointerpretable than the Deep Neural Network (DNN) [2]. The information transmission method of SNNs take the spike form, and enables it to obtain the unique temporal information processing capability, while the sparse firing also gives it low power consumption [3]. Therefore SNNs are considered as network models to achieve more advanced artificial intelligence. However, SNNs that rely solely on neuron dynamics models and synaptic connections have a limited information capacity, and The associate editor coordinating the review of this manuscript and approving it for publication was Qichun Zhang . may have spike errors or even spike disappearance in complex networks due to the decay of spikes during the transmission process, resulting in disconnected neurons not being able to effectively transmit distal information and thus making it challenging to achieve advanced brain-like intelligence.
As brain neuroscience research progresses, researchers have discovered and demonstrated that the biological brain relies on more than just neuronal and synaptic networks to process and integrate information as well as achieve advanced intelligence functions. Glial cells, which were once thought to only provide support for the biological activities of the brain, are likely to be one of the key players in the realization of higher brain intelligence functions [4]. The feasibility of such approaches has been corroborated by researchers who have drawn on the regulatory mechanism of brain astrocytes • Proposed a novel spiking neural network, BGSNN, with a global information interaction mechanism and a diverse network structure. It supports global spatiotemporal domain information integration and can plasticize neurons and synapses. This enables feed-forward and feedback optimization of the network, enhancing its expressiveness. To verify the performance of BGSNN, we designed and implemented a Sudoku solver application for BGSNN based on the ''WenTian'' neuromorphic prototype, evaluated on the Easybrain and the kaggle dataset. The experiments show that the BGSNN has excellent parallelism and sparsity performance with the Sudoku solving accuracy exceeding 99% over all datasets.
This paper is organized as follows: Section II introduces the related work; Section III introduces the method in detail; Section IV reports the experimental results and analysis and the conclusion is given in Section V.

II. RELATED WORKS A. GLIAL CELL BIOLOGY RESEARCH
Glial cells are the collective name for a class of neuron cells in the brain. Existing research suggests that the glial cells are effectual for the brain far more in the aspect of offering structural support of metabolic substances for the brain activities.
Many studies have found that there is information interaction and regulation between glial cells and neurons through messenger substances such as molecules, ions, or proteins. For example, the receptors of glial cells will affect the synaptic transmission activity of neurons in activated state; thereby affecting the long-term memory of neurons in the hippocampus [15]. Neuronal activity can significantly affect the activity of glial cells, and the cholesterol complexed to apolipoprotein E-containing lipoproteins, and adenosine provided by glial cells can in turn regulate the growth and development of brain synapses [16], [17]. It shows that the whole process of the growth and development of the brain neural network is the result of the interaction between the activity of neurons and glial cells.
In addition, researchers have found that information interaction between glial cells also exists through messenger substances such as molecules and ions. Glial cells enhance intercellular communication through the diffusion of Ca 2+ in the glial space [18]. The communication between glial cells is the basis for the brain to distinguish different cognitive properties and is a key element for the brain to achieve cognitive functions [19], [20].
The researchers found that brain glial cells and neuronal synapses are physically connected and affect the function and efficiency of synapses. Neuronal synapses have three-dimensional connection structures of glial protrusions [21], [22], which adhere to the synapses through Ca 2+ permeable glutamate receptors to maintain structural stability [23]. Glial protrusions are capable of transmitting a variety of molecules and ions, such as K + , amino acid transmitters, tumor necrosis factor α (TNFα), ATP, and adenosine [24], [25], [26], which can affect the function and efficiency of synapses, in order to regulate neuronal activity and integrate information in a spatiotemporal domain complementary to neurons [27], [28], [29], [30].

B. SNN PLASTICITY
Most of the existing SNN learning algorithms are only related to synaptic parameters, such as synaptic weights, whereas the neuron related parameters in the network are often defined as hyperparameters, which to some extent limits the expressiveness of SNNs. It has been pointed out that there are differences in membrane time constants of neurons in different brain regions [31], [32], [33], and these differences play a crucial role in learning and memory work [34], [35]. Fang et al. proposed to learn membrane time constants along with SNN synaptic weights to achieve enhanced SNN learning [36]. This shows that SNN plasticity can be reflected not only in synapses but also in neurons.

C. GLIAL CELLS REGULATE NEURAL NETWORKS
The application of glial cells in neural networks has gradually been noticed by researchers, Tang et al. implemented a triple synapse model based on glial protrusions on Intel's neuromorphic chip Loihi, and some special ways of biological brain information processing are realized, like synchronous excitation and local plasticity [5]. Ivanov and Michmizosproposed a neuron-astrocyte liquid state machine (NALSM) [6], which addresses the low performance of LSM through self-organized near-critical dynamics and verifies that the brain-inspired machine learning methods have the potential to achieve comparable performance to deep learning while supporting robust and energy-efficient edge computing. At present, there is no precedent for using glial cells and their network structure in SNN for parameter learning and global spatiotemporal information integration in SNN. Therefore, the purpose of this paper is to apply the glial cells network structure in SNN and elaborate the subjected research.

D. NEUROMORPHIC COMPUTING PLATFORM
Neuromorphic computing has become a new generation of research boom, and several excellent neuromorphic computing platforms such as SpiNNaker [7], Loihi [8], TrueNorth [9] and Tianjic [10] have emerged in the market. SpiNNaker, which designed by the University of Manchester, is a massively parallel, many-core supercomputer architecture based on a spiking neural network, able to simulate the human brain, providing new possibilities for neuroscience, robotics, and computer science. ''Tianjic'', developed by the Tsinghua University, is the world's first heterogeneous fusion brain-like computing chip, adopts a many-core architecture, reconfigurable components, and a streamlined dataflow with hybrid coding schemes, can support both machine learning algorithms and existing brain-like computing algorithms. This chip, integrating neuroscience and computer science, is expected to enhance the capability of each system, promote the research and development of artificial general intelligence (AGI), and provide a hybrid collaborative development platform for AGI technology.
In this paper, the ''WenTian'' neuromorphic prototype developed by Nanjing Institute of Intelligent Technology is selected for the simulation design and implementation of BGNN. The prototype consists of thirty FPGA boards installed on the cabinet in three layers of 10 boards each, forming a 3 × 10 torus (circular surface) topology network. Each individual node in the network, i.e., each board, is a many-core system with the ability to work independently. The ''WenTian'' neuromorphic prototype can support 480 Cortex M4 processors running in parallel, enabling the simulation of millions of neurons and billions of synapses. It also has a supporting programmable simulation softwae system, which only needs to add new network nodes and connection types based on SNN, and adapt the original neuron model to complete the BGSNN simulation environment. Therefore, it has the simulation ability of BGSNN.

III. BLENDED GLIAL CELL'S SPIKING NEURAL NETWORK
Existing artificial neural networks are composed of two basic units: neurons (network nodes) and synaptic connections (network connections). Inspired by the advanced characteristics of glial cells, a Blended Glial Cell's Spiking Neural Network (BGSNN) is proposed in this paper. In addition to neurons and synapses, BGSNN introduces glial cells and four corresponding network dynamics connection models, which expand the original single information processing dimension of SNN. Glial cells can not only act as both spatiotemporal information processing units to process and transmit information but also as modification units to modify neuronal and synaptic states in the network based on supervised signals and global spatiotemporal information.

A. NETWORK STRUCTURE OF BGSNN
Organization of all nodes and connections in BGSNN is diagrammed in Fig. 1.
Inspired by biological brain research, BGSNN adds a special information processing unit called the glial cell, as well as four types of connections between neurons or synapses: neuronal ion connection, glial ion connection, glial gap connection, and glial protrusion connection.
Neuronal ion connections are the channels working between neurons and glial cells, where neurons can release messenger factors to glial cells, and neuronal ion receptors can receive messenger factors transmitted by neuronal ion channels [15], [17]. Glial ion connections are the channels also working between neurons and glial cells, that glial cells can release regulatory messenger factors via channels to neurons, and glial ion receptors can receive the messenger factors transmitted by glial ion channels [4], [27], [30]. Glial gap connections are the channels working between glial cells that there has the mutual transmission of messenger factors between glial cells, and glial gap receptors can receive the messenger factors transmitted by the glial gap channels [18], [19], [20]. Glial protrusion connections are the channels working between glial cells and neuronal synapses that glial cells can transmit regulatory messenger factors to neuronal synapses [21], [22], [23], [24]. Fig. 2 illustrates the general structure of the BGSNN, including the BGSNN neuronal cell populations and connections relationships. in BGSNN, each glial cell population is both a recognizer and a sharer of local area information.

B. DYNAMICS MODELS OF BGSNN
As scientists make further study on glial cells, corresponding dynamics models of glial cell were proposed which mainly started from molecular and ionic changes [37]. These models, however, oversimplified the mechanism of glial cell activity and hardly reflected the advanced functions of glial cells in information processing. Therefore, the current glial cell biodynamic model is not suitable for directly applied to SNNs, and necessary to build dynamics models by combining the biological properties of glial cells with the engineering properties of neuromorphic computing.
As shown in Fig. 3, this paper gives a general glial cell dynamics model framework, which contains three parts: the neuronal spike event processing module, the glial gap event processing module, and the spatiotemporal information   integration processing module. The main structure of this framework is consistent with the biological structure of glial cells, and the designer can customize the dynamic algorithm of each module according to different engineering requirements. Considering that BGSNN adds the structure of glial cells, the framework of neuron dynamics model will change accordingly, as shown in Fig. 4. Based on the original neuronal spike event processing module, the neuron dynamics model framework adds three parts: a glial ion event processing module, a glial protrusion event processing module, and a spatiotemporal information integration and processing module. The neuron dynamics model and the glial cell dynam-ics model constitute the computing model of BGSNN. The details of two models will described in the next subsection.
The glial cell dynamics model is a time-driven model, different from the event-driven model of neurons. Glial cells perform corresponding operations, such as making state judgments and spike firing, within a specific time cycle. Its information processing and responding time cycle is shown in Fig. 5. A complete cycle T g is divided into four subperiods: T ng (information processing period from neural network to glial cells), T gg (information processing and response period from glial network to glial cells), T gn (information response period from VOLUME 11, 2023  glial cells to neural network), T re (glial cell non-response period).
During the T ng the period, the glial cells only process the input neuronal spiking information, i.e., obtain the activity information of local neuron populations; during the T gg period, glial cells diffuse valuable activity information of local neuron populations through glial gaps, and receive activity information of neuron populations diffused from glial cells in other regions; during the T gn period, glial cells integrate the activity information of the local neuronal population with those of other regions, and regulate the local neuron population activity accordingly; during the T re period, glial cells are resting.
Taking a single-layer feedforward neural network as an example, as in Fig. 6, the BGSNN network computation model can be expressed mathematically as follows according to dynamic models.

a. Neuronal Spike Event Processing Module
This module operates during the T ng period and is responsible for processing the neuronal spike events received by glial cells through the neuronal ion connection channel. And the glial cell obtain the input current in T ng period through this connection channel, as shown in (1): where I ng (t) denotes the input current value of neuronal ion channels in the glial cell; w ng denotes the weight of the neuronal ion connection; and S i denotes the output spike of neuron.

b. Glial Gap Event Processing Module
This module operates during the T gg period and is responsible for processing glial gap events received by the glial gap connection channel, and the glial cell obtain the input current in time period T gg through this connection channel, as shown in (2): where I gg (t) denotes the input current value of the glial gap channel of the glial cell; w gg denotes the weight of the glial gap connection; and O gg denotes the glial gap output event of glial cell, which is the output by the spatiotemporal information integration processing module of the glial cell. c. Spatiotemporal Information Integration and Processing Module The module can output three different events in the form of spike, including glial ion events, glial gap events, and glial protrusion events. The output results can be expressed in Fig. 7 as (3): where O gg denotes the output of glial protrusion channels, i.e., glial protrusion events; O gn denotes the output of glial ion channels, i.e., glial ion events; O gp denotes the output of glial gap channels, i.e., glial gap events; I ng denotes the input current value of the neuronal ion channel in glial cells; and I gg denotes the input current value of the glial gap channel in glial cells. where: • Glial ion events: generated during the T gn period, glial cells transmit glial ion events to neurons through glial ion connection channels, which can be used to modulate neuronal plasticity.
• Glial gap events: generated during the T gg period, glial cells diffuse glial gap events between glial cells through glial gap connection channels, which can be used for global information interaction.
• Glial protrusion events: generated during the T gn period, glial cells transmit glial protrusion events to neurons via glial protrusion connection channels, which can be used to modulate the plasticity of neuronal synapses.

2) NEURON DYNAMICS MODEL
a. Neuronal Spike Event Processing Module As the spatiotemporal information computing unit, after the neuron receives the spike event. The input of the neuronal spatiotemporal information integration and processing module is obtained through calculation,   which can be expressed as (4): where w i represents the weight of the ith fan-in synapse, S i representing the input spike event.

b. Glial Ion Event Processing Module
After the neuron receives the glial ion event through the glial ion connection channel, it will correspondingly change its state parameters to realize neuron plasticity, as shown in (5): where P n (t) denotes dynamic parameters of the neuron adjusted by glial ions at time t; w gn i denotes the weight of the glial ion connection; O gn i denotes the glial ion output event of the glial cell, given by the glial ion event in the spatiotemporal information integration processing module of the glial cell model. c. Neuron Spatiotemporal Information Integration and Processing Module This module receives the processed neuron spike input and neuron dynamic parameters adjusted by glial ion events, then integrates and processes their information. As shown in Fig. 8, outputs neuronal spike events can be expressed as (6):  where O n (t) represents the neuronal spike event output by the neuron; P n (t) represents dynamic parameters of the neuron adjusted by glial ions; N n (t) represents the neuron spike input; F represents the neuron model.

d. Glial Protrusion Event Processing Module
After the neuron receives the glial protrusion event through the glial protrusion connection channel, it will adjust the weight of the neuron synapse, i.e., the plasticity of the neuronal synapse, as shown in Fig. 9 can be illustrated by (7): where W (t) represents the synaptic weight of the neuron at time t; W (t −1) represents the synaptic weight of the neuron at time t-1; O gp denotes the glial protrusion event of the glial cell, given by the spatiotemporal information integration processing module of the glial cell model. w gp denotes the weight of glial protrusion connection.

C. SPATIOTEMPORAL INFORMATION PROCESSING INTEGRATION AND PLASTICITY PRINCIPLE OF BGSNN
This section will describe the mechanism of BGSNN for information processing and integration in the spatiotemporal domain, and explain how BGSNN implements the neural node and connection modification in the network. The spike transmission of BGSNN is shown in Fig. 10. Without glial cells, spikes can only be transmitted from presynaptic neurons to postsynaptic neurons via axons, synapses, and dendrites, and feedback transmission is impossible. With the introduction of glial cells, spikes can be secondary processed by glial cells and transmitted by feedforward or feedback networks, giving the networks the ability to integrate spatiotemporal information and plasticity.
In feedforward networks, spikes can be transmitted to fan-out synapses or postsynaptic neurons by feedforward networks1 or 2 via glial gap connections and glial cells, giving the BGSNN the ability for secondary deep processing of spatiotemporal information and feedforward optimization or plasticizing of the network.
In feedback networks, spikes fired by the postsynaptic neuron can be transmitted to the fan-in synapse via feedback network 1, which is composed of glial cells and glial protrusion connections, to adjust the parameters of the fan-in synapse; or to the presynaptic neuron via feedback network 2, which is composed of glial cells and glial gap connections, to adjust the parameters of the presynaptic neuron.
Therefore, BGSNN has the capabilities of feedforward and feedback optimization or plasticizing of the network. The processing flowchart of the BGSNN for spike events is shown in Fig. 11. For a complete glial cell, the neuronal ion channel is the direct sensing channel for short distance coupling of spatiotemporal information; the glial gap channel is the indirect sensing channel for distant coupling of spatiotemporal information; the glial cell is the smallest unit for integration and processing global spatiotemporal information of the network; the glial ion channel or glial protrusion channel is the direct optimization channel of glial cells to the neuronal network.
To summarize the BGSNN, it adds the glial cell as a special information processing unit and four new types of connections between neurons or synapses. These changes enable BGSNN to have the capabilities of global spatiotemporal information integration, feedforward and feedback network optimization, enhancing the expression ability and the coupling of networks. Because of the addition of the glial cell structure, BGSNN supports not only the general SNN training method but also the direct network adjustment through network feedback, so that inference can be performed directly without training.

IV. EXPERIMENTS AND RESULTS
In this section, an application example is designed and implemented to demonstrate the spatiotemporal information integration and plasticity of BGSNN. According to the application requirements, the glial cell dynamics model and neuron dynamics model are instantiated, and the details can be found in Appendix A.

A. BGSNN APPLICATIONS-SUDOKU PUZZLE SOLVER
Sudoku is an NP-complete problem [38], as shown in Fig. 12. The Sudoku board is composed of 9 blocks, each of which in turn consists of 9 cells. The rule is to infer the numbers of all remaining spaces based on the known numbers on the 9 × 9 board and to satisfy that the numbers in each row, column, and block contain 1-9 without duplication. Every qualified Sudoku puzzle has one and only one answer, and the inference is based on this; any unsolved or multiple-solution puzzle is ineligible.
When solving Sudoku problems, people try to fix a number in a cell choosing from 1 to 9 according to the rules. This process can be thought of as accumulating each number's possibility in the cell, and the cell will be filled with the number with the highest probability. This process of accumulation fits well with the computation mechanism of spiking neural networks. At the same time, Sudoku problems are a class of problems that require local and global synergy, which can highlight the performance improvement of glial cell structure brought by global spatiotemporal information integration and plasticity. Therefore, we designed a Sudoku puzzle solver based on BGSNN as a case study to demonstrate the implementation of BGSNN for intelligent decision-making and its effectiveness. The software and hardware of the neuromorphic platform are designed for spiking neural networks with good parallelism, so to ensure the optimal performance of the BGSNN solver, we deployed the BGSNN solver on the ''WenTian'' neuromorphic prototype for experiments.
The population structure of the BGSNN Sudoku solver network is shown in Fig. 13A, which consists of a total of three populations: 1) the Sudoku initial input neuron population; 2) the Sudoku board neuron population; 3) the Sudoku analysis and decision-making glial cell population.
The Sudoku initial condition input neuron population establishes a one-to-one excitatory connection (one_to_one_excite(input)) to the Sudoku board neuron population, allowing the initial condition to be represented on the Sudoku board. For the Sudoku board neuron population, describes rules of the Sudoku game by establishing inhibitory connections (sudoku_rule_inhibit) within the populations. For the Sudoku analysis and decision-making glial cell population, enables glial cells to make a rule-based analysis of the current state of the Sudoku board by establish excitatory connections (sudoku_rule_excite) internally that satisfy the rules of Sudoku. The one-to-one excitatory connection (one_to_one_excite) established from the FIGURE 13. BGSNN Sudoku solver. A) Sudoku solver network population structure. B) Sudoku solver network structure (with the number ''1'' as an example): the Sudoku initial input neuron population will input the initial condition ''1'' into the known cell, then the ''1'' in the same block, line and column will be inhibited, and the remaining numbers in the cell except ''1'' will also be inhibited. The ''1'' in this cell transmits excitatory signals to the corresponding glial cell neurons, causing the glial cells to adjust accordingly, i.e., the ''1'' in the same block, line, and column is excited, and the rest of the numbers in the cell except for the ''1'' are also excited. And it also receives excitatory signals from the corresponding glial cells to make adjustments to the Sudoku board.
Sudoku board neuron population to the Sudoku analysis and decision-making glial cell population can transmit the state of the Sudoku board to glial cells for the perception of the local state of the Sudoku board by glial cells. The oneto-one excitatory connection (one_to_one_excite (adjust)) established from the Sudoku analysis and decision-making glial cell population to the Sudoku board neuron population enables to regulate the state of the Sudoku board neuron population. The Sudoku solver network structure (with the number ''1'' in a cell as an example) is shown in Fig. 13B. The BGSNN Sudoku solver can directly infer to obtain solution results without network model training. Since there is no published research data related to SNN solving Sudoku, this paper selects Sudoku solvers of traditional algorithms, deep learning algorithms and DNN methods as experimental control groups for comparison experiments.

1) EASYBRAIN DATASET
The EasyBrain website provides five difficulty levels of Sudoku puzzles: easy, medium, hard, expert, and evil.   Sudoku puzzle solving process: A) Sudoku puzzle initial conditions. B) spiking activity of Sudoku puzzle solving. Until a solution to the Sudoku puzzle is solved, glial cells will fire spikes with Tg as the interval. The neurons that correlate to the board will stay active throughout the solving process, while glial cells will help the equilibrium state arrive more quickly. When the equilibrium state is attained, only the neurons matching the solution's corresponding number will fire, leaving the other neurons inactive. C) solution of Sudoku puzzle. solver; and to verify the superiority of parallel computing of BGSNN.

a. Glial Cells Necessity Evaluation
A SNN Sudoku solver, as a control group, was designed based on the network structure of the BGSNN Sudoku solver. This SNN network only lacks the glial cell module compared to the BGSNN network, and the schematic diagram of the network structure is shown in Fig. 14.
We took a Sudoku puzzle of hard difficulty level as an example as shown in Fig. 15. The Sudoku anal-ysis and decision-making glial cell population were adjusted for 15 iterations to finally bring the Sudoku board neuron population to the expected state, i.e., the correct solution.
The experiment uses the Easybrain dataset to run the BGSNN solver and the SNN solver respectively, then recorded the average computational accuracy of Sudoku solvers at different difficulty levels, as shown in Table 1.
It can be seen from the Table 1 that the SNN Sudoku solver is difficult to solve correctly at the Hard level and above, while the BGSNN Sudoku solver can solve correctly at all difficulty levels under the same conditions with a stable accuracy of over 99%. Evil had an average of 1.71 fewer clue values than Expert difficulty level, but the SNN Sudoku solver accuracy was reduced by 20.41%. It is proved that the accuracy of SNN solver without global regulation depends greatly on the number of initial clues. The difficulty in transmitting information throughout the entire board to achieve overall balance without global regulation may be the cause of this situation since there are too many cells to fill in and a great distance to cover.
The only difference between the two Sudoku solvers was the glial cell structure. Therefore, the experiments showed that the global information processing and the plasticity of neurons brought by glial cells can effectively help the network converge to the correct answer, which verified the necessity of the glial cell structure. b. Evaluation of accuracy over runtime The BGSNN Sudoku solver is a time-driven network solver, and the network completes the solution task through periodic iterations. To show how the accuracy of the BGSNN Sudoku solver varies with runtime, we set the runtime to ten intervals, each interval is 30 in length, and a runtime has 1ms step size. Using the Easybrain dataset, for each difficulty level, the accuracy with which the BGSNN solver can fully solve for the correct answer in each interval is counted. The BGSNN Sudoku solver was tested on the dataset with 30 puzzles per difficulty level, for a total of 750 trials (each puzzle was repeated 5 times). Fig. 16 shows the accuracy of the BGSNN solver in each interval that the answer can be completely and correctly solved. The accuracy of the BGSNN solver increases with the increase of runtime, and the most difficult level, Evil, requires at most 300 runtimes to complete the solution. Since the probabilistic stochas-tic down-sampling function in the glial cell dynamics model used in the experiment (see Appendix A for details), there is stochastic in each initialization, and poor initialization will lead to an increase in solution time. Therefore, the accuracy curve is not completely smooth, but maintains an upward trend. As with other iterative algorithms, the accuracy of BGSNN increases with the number of iterations until the task is completed. c. the BGSNN Sudoku Solver Performance Evaluation The BGSNN Sudoku solver is compared with SOTA's Sudoku algorithm for performance evaluation. LSGA [41] is an evolutionary algorithm to solve Sudoku and performs well in terms of accuracy and convergence speed. The dataset used by LSGA in the original paper is WebSudoku [42], which has four difficulty levels: Easy, Medium, Hard, and Evil, and the mean values of clues for each difficulty level are 36.25, 30.75, 27.63, and 25.62, respectively. According to the difficulty levels, this dataset can correspond to the first four levels of Easybrain.
In the original paper, the number of generations needed by LSGA to obtain the optimal solution is used as an evaluation criterion, and this data is independent of the hardware. BGSNN is a time-driven network solver that iterates according to the period, and each iteration updates the network state and parameters until the correct solution is found, while the number of iteration periods is hardware-independent. In the comparison experiment, we consider one complete cycle of BGSNN as a generation and evaluate the performance with 100% accuracy guaranteed.
In this experiment, the BGSNN solver was tested on the dataset for a total of 750 times (5 repetitions per puzzle). The test results were compared with the published test results of the LSGA solver, as shown in Table 2. From the comparison results, it can be seen that the LSGA solver requires a higher number of generations to find the correct solution than the BGSNN Sudoku solver at all difficulty levels, and the BGSNN can solve puzzles at higher difficulty levels. This proves that the BGSNN Sudoku solver can converge faster with guaranteed accuracy.

d. Parallelism Evaluation
The BGSNN Sudoku solver experiments run on the ''Wentian'' neuromorphic prototype, which supports parallel computing. Parallel computing has significant advantages over serial computing in multi-loop complex problems, saving significant runtime and reducing algorithmic time complexity. The backtracking algorithm uses depth-first search to solve the Sudoku puzzle, traversing all the cells until the correct solution is obtained, which is a typical serial logic algorithm; BGSNN supports parallel computing and can calculate multiple neuron spikes at the same time. In this paper, a set of control experiments were designed to compare the number of loops required by the solver when solving the same Sudoku puzzle to verify the advantages of parallel computing in time complexity. There was no random function in the backtracking algorithm solver, and there was only one loop count when a puzzle was completely solved. Since there was a probabilistic random downsampling function in the spatiotemporal information integration and processing module of the glial cell dynamics model used by the BGSNN solver, it may have a different number of loops for the same puzzle. To make the results comparable, for one puzzle, the experiment was repeated 10 times using the BGSNN solver, and the average number of loops was compared with the results of the backtracking algorithm solver. In particular, the glial cell dynamics model in BGSNN was time-driven, and one timestep could correspond to one loop in the backtracking algorithm during hardware execution. Experimental results are shown in Table 3.
In this experiment, to ensure the accuracy and stability of the solution results, we assumed that the BGSNN Sudoku solver correctly solved this puzzle when the same result could be output stably for three consecutive glial cell response cycles (T g ). Taking this experiment as an example, T g was set to 900ms, so the number of BGSNN in the table was equal to the number of loops required for the BGSNN solver to obtain the correct solution, plus the length of the two additional response cycles, that is, 1800. As can be seen from Table 3, BGSNN can reduce the computation by at least 92% through parallel computing, which helps improve the computing speed and reduce the computing cost.

2) KAGGLE DATASET
This dataset contains 1 million Sudoku puzzles with corresponding solutions. After statistics, the maximum clue value for this dataset was 37, the minimum was 29, and the average was 33.81279. In this paper, three sets of controlled experiments are conducted to verify the accuracy of the BGSNN Sudoku solver on a large scale dataset; to verify the performance of BGSNN compared with DNN on the same scale dataset; to verify the superiority of BGSNN sparse computing.
a. Accuracy Evaluation of the BGSNN Sudoku Solver The experiment ran the BGSNN Sudoku solver on this 1 million scale dataset, recorded the result of each solution and the solution accuracy. The correct solution rate is shown in Fig. 17.
According to the experimental findings, the solution accuracy of the BGSNN Sudoku solver can reach 99.9999% with this dataset. In the 1 million-scale dataset, the solution accuracy of two puzzles was less than 100%, which were 90.1235% and 85.1852% respectively. In this experiment, the upper limit of the running time was artificially set to 20000ms. If the upper limit was reached, the solver would automatically terminate. Meanwhile, there is a probabilistic stochastic down-sampling function in the glial cell dynamics model used in the experiment (see Appendix A for details). When the initialization is poor, the solution time will increase accordingly.
In response to the above, three additional control experiments were conducted for the two puzzles, and the statistical solution accuracy was 100%. Therefore, it was analyzed that in the first large-scale validation experiment, the case where the accuracy does not reach 100% may be due to the poor random initialization, resulting in the runtime exceeding the upper limit. b. the different Approach Sudoku Solver Performance Evaluation Since the BGSNN Sudoku solver is a network approach solver, several representative Sudoku solvers are chosen from the various Sudoku solving algorithms released in the kaggle dataset to evaluate the performance of the BGSNN solver.  The accuracy of the backtracking solver is up to 100% [43] with O(n 7 ) time complexity, but its computation is substantially greater than that of the BGSNN solver. This analysis has been described in subsection c of the Easybrain dataset.
''Solve Sudoku with CNN acc 97%'' [44] is currently published the DNN solver with the highest accuracy except for the backtracking algorithm. As the control group, it used an 11-layer network structure, including 8 convolutional layers, 2 fully connected layers, and 1 softmax classification layer. The solver used 12 epochs for training, each epoch size was 12500, and the final training accuracy is 0.9617. The compared results are shown in Table 4.
The table demonstrates that BGSNN outperforms DNN in terms of accuracy, the number of layers, and parameter quantity. The final solution accuracy of the DNN model depends on the training data set and training algorithm. If the data set or training algorithm is inappropriate, the final accuracy will decrease. BGSNN is different from DNN, the model supports directly infer, so the final solution accuracy does not depend on the training process. When the data samples are insufficient, BGSNN can effectively resolve the problem that the model cannot be trained, and provides a new solution for small sample scenarios.
In addition, we chose ''Sudoku Solver -97.45% Accuracy on 1M Games'' [45], which employs a systematic approach and achieves 97.45% accuracy without using neural networks, and its time complexity is O(n 6 ). However, this is still less than the 99.99% accuracy of the BGSNN Sudoku solver. c. Sparsity Evaluation BGSNN transmits information in the same form of spikes as SNN, and the sparse spike firing makes its low computation. The BGSNN Sudoku solver runs on the ''Wentian'' neuromorphic prototype, which supports sparse matrices instead of full matrices in computing. Ran the BGSNN Sudoku solver on this 1 million size dataset, and recorded all its sparse matrices involved and the corresponding full matrices. The difference between the average computation of the two was counted as the reduction of computation by sparse computing of BGSNN compared with the equal fully dense DNN. The experimental results are shown in Table 5.
According to the experimental findings, BGSNN can reduce roughly 88% of the computation through sparse computing, which contributes to lowering computing costs and increasing computing speed.

V. CONCLUSION
In Summary, this paper proposes a novel Spiking Neural Network, called Blended Glial Cell's Spiking Neural Network (BGSNN), which is inspired by the glial cells. BGSNN introduces glial cells and four corresponding network dynamics connection models, thus realizing the modification of neural node and synaptic connection, improving the capability of the network in deep local information processing and global information integration, and enhancing the expression abilities of the network.
In experiments, based on the ''WenTian'' neuromorphic prototype, this paper instantiates dynamics models of glial cell and neuron based on the application scenario, and design the BGSNN Sudoku solver and series of experiments to evaluate on the Easybrain and kaggle datasets. The solution accuracy of the BGSNN solver exceeds 99% on all datasets. With 100% accuracy, the BGSNN solver improves accuracy by over 97% over the same structure SNN solver at the evil difficulty level in the Easybrain dataset, demonstrating the need for glial cells. In the performance evaluation, compared with the SOTA Sudoku solver LSGA, the BGSNN Sudoku solver has a faster converge speed. It also achieves over 99% accuracy on the million-level dataset, which is 3.82% better than the publicly available optimal DNN solver on the same dataset. BGSNN performs well in sparsity and parallelism experiments. BGSNNN improves sparsity by about 88% over the equal fully dense DNN and reduces computation by at least 92.9% compared to the serial logic algorithm.
During the experiments, it has been demonstrated that in the absence of multi-dimensional information processing mechanisms and the global information interaction mechanism, neurons can only communicate with the neurons connected to them. When the network size is large, this kind of transmission process limits disconnected neurons from efficiently transmitting the distal information since there is attenuation in the transmission of spikes. BGSNN has the ability of global regulation and secondary processing of information while solving the bottleneck of the single-dimensional information processing mechanism, which can significantly improve the quantity and quality of information. And it has a more diversified network structure while maintaining the advantage of SNN parallel sparsity, which can effectively improve the expression ability of the network and realize more advanced artificial intelligence.

APPENDIX A INSTANTIATION OF BGSNN
The BGSNN variant GLIF model of the LIF (Leaky integrateand-fire model) model [46] was chosen as the neuronal dynamics model applied in the experiments of this paper, and a mean disbursement strength encoding bandpass conduction error feedback glial dynamics model (mEIC_bPT_eFB) was proposed. The parameter setting of this model is shown in Table 6.

A. GLIF MODEL 1) NEURONAL SPIKE EVENT PROCESSING MODULE
As the spatiotemporal information computing unit, after the neuron receives the spike event, the spike input of the neuronal spatiotemporal information integration and processing module is obtained through calculation, which can be expressed as (8): where w i represents the weight of the ith fan-in synapse, S i representing the input spike event.

2) GLIAL ION EVENT PROCESSING MODULE
After the neuron receives the glial ion event through the glial ion connection channel, it will correspondingly change its state parameters to realize neuron plasticity, as shown in (9): where P n (t) denotes the dynamic parameters of neurons adjusted by glial ions at time t; P n (t − 1) denotes the dynamic parameters of neurons adjusted by glial ions at time t-1; w gn i denotes the weight of the glial ion connection; O gn i denotes the glial ion output event of the brain glial cell, given by the glial ion event in the spatiotemporal information integration processing module of the glial cell model.

3) NEURON SPATIOTEMPORAL INFORMATION INTEGRATION AND PROCESSING MODULE
The mathematical expression of the module is shown in (6), where the neuron model is chosen as the LIF model, whose first order differential equation is described in (10): where τ m = R m C m is called the membrane time constant; I is the sum of the synaptic currents generated by the firing behavior of the individual presynaptic neurons. When the membrane potential V is greater than or equal to the threshold potential, the neuron immediately generates excitation and fire a spike that accompanies the conduction of the action potential, while resetting the membrane potential to V reset , and holding it during the absolute refractory period (T r ).

4) GLIAL PROTRUSION EVENT PROCESSING MODULE
After a neuron receives a glial protrusion event through the glial protrusion connection channel, the weight of the neuron synapse changes as shown in (11): where W (t) represents the synaptic weight of neurons at time t; W (t − 1) represents the synaptic weight of neurons at time t-1; O gp denotes the glial protrusion events of glial cells, given by the glial protrusion events in the spatiotemporal information integration processing module of the glial cell model. w gp denotes the glial protrusion connection weights.

B. mEIC_bPT_eFB MODEL
Based on the framework of brain glial cell model and combined with practical engineering experience, this paper proposes a mean disbursement strength encoding bandpass conduction error feedback glial dynamics model (mEIC_bPT_eFB). The driving method of this model adopts a time-driven approach, i.e., the state judgment and event firing are made only when a specific time period arrives.

1) GLIAL PROTRUSION EVENT PROCESSING MODULE
I ng denotes the input current value of neuronal ion channels in glial cells, as shown in (12): where t denotes the current time, T ng denotes the information processing period from the neural network to the glial cells; T gg denotes the information processing period from the glial network to the glial cells; T gn denotes the information processing period from the glial cells to the neural network.

2) GLIAL GAP EVENT PROCESSING MODULE
I gg denotes the input current value of glial gap channels in glial cells, as shown in (13): I gg = 0, t%T g < T ng ∪ t%T g > (T ng + T gg ) II (t, T gg ), t%T g ≥ T ng ∩ t%T g ≤ (T ng + T gg ). (13) where t denotes the current time, T ng denotes the information processing period from the neural network to the glial cells; T gg denotes the information processing period from the glial network to the glial cells; T gn denotes the information processing period from the glial cells to the neural network.
The channel input current values II (t, T ) in (12) and (13) are used to extract information about the input spike events, which are encoded by the average release intensity as shown in (14): where ⃗ s(t) is a 1 × k-dimensional vector representing the k fan-in connections of the glial cell at the current moment of binary spike input, where the convention is ''0'' for no spike and ''1'' for the opposite. ⃗ W is a k × 1-dimensional vector representing the weight size of the k fan-in connections of the glial cell, and T indicate the time window size.

3) SPATIOTEMPORAL INFORMATION INTEGRATION AND PROCESSING MODULE
The mathematical expression of the mEIC_bPT_eFB model is (15), as shown at the bottom of the page, where ⃗ O glial is the output vector of glial cells, O gg for the output of glial gap channels, O gn for the output of glial ion channels, and O gp for the output of glial protrusion channels.
V m (I ng , I gg ) in (15) is the information integration function of glial cells, which is used to represent the membrane potential of glial cells, as shown in (16), where R ng is the neuronal ion channel input resistance and R gg is the glial gap channel input resistance. T %T g > (T ng + T gg ) mSIC(x) in (15) represents the average pulse release intensity encoding function, as shown in (17): where x denotes the value to be encoded at the current moment of the type of channels, and pS(x) is a function that determines with probability whether to fire a spike or not, and is used to encode the information to be delivered in the form of a spike event, as shown in (18). Where V sa denotes the saturation threshold and RAND(0, 1) is a function that generates a random number in the range (0,1).
pS(x) = 0, x/V sa < RAND(0, 1) 1, x/V sa ≥ RAND(0, 1). (18) bP(x) in (15) is a bandpass filter function for processing local spatiotemporal domain information i.e., to obtain valid local neuronal population activity information and diffuse it to other regions of the glial cells. The trapezoidal filter function is used here, as shown in (19), by setting the a, b, c, d parameter, the glial cells can select or filter the specific neuron activity information.
In (15), eFB(i, v) is the error feedback calculation function, which is used to process the global spatiotemporal information, and in turn adjust the local neuron population activity, as shown in (20); where S(i) is the switching function, as shown in (21), I th is the effective current threshold of the glial gap input; where D(v) is the error correction function, as shown in (22), P is the feedback polarity coefficient, and V e is the expected average spike firing intensity value, V threshold is the minimum effective information threshold.