Heterogeneous axonal delay improves the
spiking activity propagation on a toroidal
network
DOCTOR OF PHILOSOPHY
Author: Salustri Marcello(215301)
Supervisor: Professor Micheletto Ruggero

Department of Physics
College of Nanoscience
YOKOHAMA CITY UNIVERSITY
Yokohama, Kanazawa-ku, Seto, 22-2, Japan
Postal code:236-0027
Year: 2024

Contents
General Introduction

5

1

Chapter 1: Methods
1.1 Network . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1.1 Izhikevich neuron model . . . . . . . . . . . . . . . . . . . . . . . . .
1.1.2 Network structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1.3 SGE: Stochastic Grid Enhancement . . . . . . . . . . . . . . . . . . .
1.1.4 Process time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1.5 SNE: Stochastic Neuron Enhancement . . . . . . . . . . . . . . . . . .
1.2 Simulation on multidimensional networks . . . . . . . . . . . . . . . . . . . .
1.2.1 Grid dimensions and SGE effect . . . . . . . . . . . . . . . . . . . . .
1.2.2 Spiking activity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Network partitioned in equivalence classes. . . . . . . . . . . . . . . . . . . .
1.3.1 mMd induces equivalence classes on the network: Octahedrons. . . . .
1.3.2 Statistical entropy: macrostate and microstate. . . . . . . . . . . . . .
1.3.3 Macrostate and microstate of octahedrons: interspike train. . . . . . . .
1.3.4 Shannon entropy . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1
1
1
3
7
12
14
15
18
23
27
27
28
30
31

2

Chapter 2: Results
2.1 Strengthening of the SGE effect. . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Network entropy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33
33
34

General Discussion and Conclusion.

41

Compliance with Ethical Standards

44

References

45

Acknowledgement

503General Introduction
BackgroundThe study of neural network systems in the brain has gained increasing importance over the
past few decades. Human cognition and behavior are governed by the brain’s neural network
system. Analyzing this system can provide insight into how information is processed and decisions are made. Based on the brain’s neural network system, artificial neural networks (ANNs)
have been developed as a key component of artificial intelligence (AI) systems [1]. Scientists
can gain insight into information processing, learning, and adaptation by studying biologicalneural networks. A greater understanding of these insights will enable ANN design and optimization to be more precise, resulting in more sophisticated artificial intelligence systems, andpromoting advancements in technology and robotics. Neurological networks in the brain haveremarkable capabilities, including pattern recognition, sensory integration, motor control, andadaptive learning [2]. Researchers and engineers have gained a deeper understanding of how
the brain works through the study of its neural network system [3]. The emergence of thesepowerful algorithms has enabled the development of robotic systems and machine learning algorithms capable of recognizing patterns, integrating sensory information, and learning fromtheir surroundings [4]. It is in this scenario that stochastic resonance (SR) plays a critical role.
SR is a phenomenon that occurs in nonlinear systems, including neural networks, where noise
can enhance the transmission of signals [5]. In the context of neural network spiking activity,
stochastic resonance refers to the idea that noise in the system can improve the network’s information processing capabilities [6]. Neural networks are made up of interconnected neurons that
communicate via electrical signals or spikes. The generation and propagation of these spikes
are influenced by various factors, including input signals and structural properties of individual
neurons [7]. The introduction of noise into a neural network can have different effects depending on the system-specific properties. Stochastic Resonance adjusts the noise floor to optimize
network performance. When the noise is at an optimal level, it can improve the network’s response to weak input signals that would otherwise be undetectable or difficult to distinguish frombackground activity [8]. Stochastic resonance in the peak activity of neural networks has been
extensively studied to understand its impact on information processing in the brain. It has been
hypothesized that noise can enhance signal recognition, enhance the presentation and transmission of sensory information, and even affect decision-making [9]. Researchers have studied themechanisms underlying stochastic resonance in neural networks and have proposed various theoretical models to explain [10] its occurrence. These models often involve complex interactions
between the network architecture, the properties of individual neurons, and the properties of input signals and noise [11]. Experimental studies have provided evidence of stochastic resonancein neural systems [12]. These studies aim to better understand how stochastic resonance affects
neural coding and information processing. Overall, the stochastic resonance in the peak activityof neural networks illustrates the complex relationship between noise and information processing in the brain. By understanding the principles underlying stochastic resonance, researchers5hope to unravel the complex dynamics of biological systems. They may develop novel strategies
to improve information processing in artificial neural networks or to treat certain neurologicaldiseases. To study the impact of noise on neural system function, we employ a novel approach inwhich we mimic the role of noise through the incorporation of axonal and neural heterogeneity.
Axonal heterogeneity refers to the variability in the physical characteristics of neuronal axons,
such as differences in axon diameter or myelination. Neural heterogeneity, on the other hand,
encompasses the diversity in neural properties, such as neuronal excitability, firing thresholds,
and synaptic strengths, within a neural population. By introducing axonal and neural heterogeneity into our modeling framework, we effectively emulate the impact of noise on neuralnetwork behavior. This approach allows us to investigate how the inherent variability in axonal
and neural properties can influence the overall dynamics, stability, and information-processing
capabilities of neural circuits. Furthermore, our work sheds light on the potential functional advantages and robustness that might be conferred by such heterogeneity in the presence of noise,
thereby contributing to a deeper understanding of the role of variability in neural systems.

Motivation for the present research work
It is undeniable that the human brain consists of an incredibly complex biological neural network system capable of performing a wide variety of calculations. Since the brain contains
approximately 86 billion neurons and 85 billion non-neuronal cells, the heterogeneity of thesecells reveals several physiological rhythms that affect its function [13]. This cell diversity enables the brain to process and store information in unique ways, which allows a wide rangeof cognitive abilities to unfold. This complexity allows the brain to adapt to new situations,change its organization, and maintain its plasticity. By exploiting these properties, scientists
have begun to apply the principles of neural networks to design artificial intelligence systems.
Modern neuroscience is concerned with the extensive development of mathematical models
that describe the biological workings of dynamic systems that resemble the brain. Because of
the large amount of brain-related data collected over time, there has been a proportional preoccupation with mathematical computer simulations and comparing them with experimental
results [14]. From a larger perspective, the brain is a complex system whose signaling takes
place in a noisy heterogeneous environment [15]. In modern neuroscience, the considerable
development of mathematical models has facilitated the realization of dynamic systems that resemble the biological functioning of the brain [16]. Among such studies, stochastic resonance,a counterintuitive mechanism by which embedded noise [17] increases a system’s sensitivity
and enhances its performance at a finite level, is believed to play an active role in a variety ofclasses of both natural phenomena and artificially structured neural systems [18]. Several studies show how stochastic resonance produces significant improvements in signal detection [19].
In this study, we aim to contribute to this area by investigating how the combination of axonal
and neural heterogeneity enhances spike propagation in an ANN. We also want to find a simple mathematical model that explains the mechanism by which heterogeneity [20] can enhancesignal propagation within neurons. Notable examples include studies on transcranial random
noise stimulation (tRNS), in which subjects are stimulated by large electrodes with weak ran6

dom electrical stimuli that actually improve their various motor, sensory, and even cognitive
tasks [21]. Given the brain’s heterogeneity, [15], we investigated the nature and effects of two
types of heterogeneity with the realistic neural model using the pyNest simulator on a toroidal
network. The initial source of heterogeneity pertains to the diversity in axonal delay, which reflects the spatial separation between neurons. The second source of heterogeneity concerns thevariability in the dynamic characteristics of neurons. In a real brain, neurons are not arranged
in a regular grid-like in a crystalline structure. Instead, their position is affected by randomness.
Since in a model network, the distance between neurons is represented by the propagation delay,in this study we determine how the increase of axonal heterogeneity reduces the time delay ofsignaling within a simple ANN. In particular, we use the Izhikevich neuron model implemented
on a toroidal network [17, 22] to study the different time delays of a signal traveling throughthe network. We manipulate the system to achieve different levels of heterogeneity by axonal
delay manipulation. Heterogeneity was implemented using a uniformly distributed parameter.
Neurons in the brain process information in limited areas that generally have no defined ends orboundaries. In this sense, choosing a torus rather than a bounding layer or plane is natural because of its boundless nature. The use of the toroidal topology in previous studies [23] together
with the mathematical topological nature of the torus, which we might want to use for future
research, led us to design our network in a ring shape. We find that the increase in axonal heterogeneity corresponds to a decrease in the propagation time of the information. This is done froman input defined as initiator nin to an arbitrary test neuron named output no . Communication
between neurons occurs via action potentials (spikes), which travel to neighboring neurons andtrigger other spikes. Overall, this creates a wave of spikes that propagate through the network.
In addition, as numerous studies, [24–27], emphasize the utility of entropy in neuroscience, thisresearch proposes a model that links axonal and neuron heterogeneity with interspike train entropy, leveraging the algebraic concept of equivalence classes to categorize brain entropy into
distinct sets. This division into subsets simplifies both theoretical and computational investigations. The model’s objective is to streamline the analysis of information within neural networks
by grouping numerous brain cells into sets based on their equivalence class.

Neuron model
Introduced in 1952 and awarded the Nobel Prize in 1963, the Hodgkin-Huxley model uses largescale nonlinear differential equations [28] to mathematically describe all plausible action potential patterns and biological membrane properties. Over the years, models have been developed
for more complex systems. A remarkable implementation of the Hodgkin-Huxley model was
published by Eugene M. Izhikevich [29]. He combines the biological results of his model with
the integration and firing of neurons in his two-dimensional system of ordinary differential equations. Izhikevich program reproduces about 20 different responses that characterize real neurons’ spiking behavior. The model’s toroidal structure emulates a finite neural domain. Specifically, each node is connected to the nearest node, and the corresponding nodes on the opposite
edge are circularly connected. Communication between each node, in a 2D grid, can be done in
the following four directions: +x or east, -x or west, +y or north, -y or south. A torus network
7

can be defined as a graph G = (N, C). Where N (G) and C(G) are the nodes and connectionsof G [30]. The total number of nodes in a 2D torus is n2 , and a 2D torus structure naturally
leads to a von Neumann neighborhood. The simulation is started by choosing a neuron (corresponding to a node in the ring network) to act as a stimulus. This is called initiator, indicatedwith ni . It is connected to a constant external current of I = 10mA which spikes periodicallythroughout the simulation period (1000 ms for most tests). Each node has a distinct axonal delayµ = cd + δ, where cd is called central delay and is δ = nd × α. The nd parameter is called noise
delay, α = [(2x) − 1], and x is a uniformly distributed random variable in [o, 1] that implements
the axonal heterogeneity. Therefore, we analyze the signal propagation in the network to thechosen no . To collect more interesting data, we chose two input-output neurons with the largest
distance between them in the structure. Only the initiator neuron spikes at the beginning of thesimulation, after which the spike activity propagates to distant test neuron no . The simulation
calculates the time it takes for the first spike to reach the chosen output node. We ran several
simulations, increasing the noise delay using the randomly uniform distribution of α parameterx. Since the dis-homogeneity is random and uniform, we should not expect any advantage dueto the noise. Instead, we will prove that it favors faster spiking activity propagation and we will
give a simple theoretical framework to understand why.

8

Chapter 1
Chapter 1: Methods
1.1

Network

1.1.1

Izhikevich neuron model

Python and the NEST library version 2.2 [31] were used to create and connect neurons that mimic
spiking networks continuously. NEST library allows the implementation of specific functions
on large sets of nodes whose connections have a configurable delay/weight as well as parametersand state variables [32]. The neuron model, which reproduces the spiking and bursting behaviorof known types of cortical neurons, is based on the one presented by Izhikevich [29]. Among
several models for studying spiking neural networks (SNN), the type proposed by Izhikevich
offers highly plausible biological dynamics, fundamental to reproducing well the non-linear phenomena dynamics of the neural network. Without realistic neuronal dynamics, we will not beable to explain the faster spiking activity propagation caused by the increased spatial heterogeneity in the network. The choice of a spiking neural model is extremely important to efficiently
structure the network task. De facto, one of the most used prototypes to describe neural functionality in the field of neuron science is the LIF (leaky integrate-fire) model, considered simple
while implemented at low computational cost [33]. However, the Izhikevich one is widely recognized as one of the most powerful and accurate models able to simulate thousands of neurons
in real-time while offering a solid accuracy to reproduce spike patterns [34]. The model offers
highly plausible biological dynamics and reproduces a large variety of spiking models by manipulating a few parameters. We used a Runge-Kutta method to process the Izhikevich model
differential equations and a pipeline/buffer procedure to describe the delay’s behavior within
the NEST network. A weighted sum considers presynaptic input, where the weights represent
synaptic connection strength parameters. When a neuron is delayed by a specific amount oftime, its signal will continue to travel down a memory register pipe. This is so long as the delay
is expressed as several time steps. In other words, at each time step, the pre-synaptic membranepotential is inserted at the beginning of the pipe, whereas the last element is fed to the Izhikevich
differential equation. Since the pipe is rolled one step forward at each time step, and because it
is as long as the delay, this method implements the time delay between pre-synaptic and postsynaptic neurons in a seamless manner without interfering with the theoretical models of delaysin this study. Furthermore, by setting a few parameters, the Izhikevich neuron model reproduces
most spiking activities, offering solid accuracy to reproduce spike patterns [34]. The Izhikevich
neuron model describes the time evolution of the membrane potential v using a two-dimensional
system of ordinary differential equations with four parameters a, b, c, and d which characterizeboth the spiking and the bursting behavior of the neurons [29]:1v ′ = 0.04v 2 + 5v + 140 − u + Iu′ = a(bv − u)after the spike, a reset process is given by:

if v ≥ 30mV →v ← cu ← u + d(1.1.1)

where v ′ = dv
, and u′ = du.
dt
dt
The variables v and u and the parameters a, b, c and d are all dimensionless [29], t is the time.
In the mathematical model, v represents the membrane potential, and u is a membrane recovery variable, providing the negative feedback to v. After each spike, according to (1.1.1) the
membrane voltage and the recovery variable are reset down to the parameter c and to u + drespectively. I is the variable representing the external current stimulus, while the adjustment
0.04v 2 + 5v + 140 allows the membrane potential v and time t to be scaled respectively to mV
and ms. The parameters a, b, c, and d, according to their values, enable the equation to depict
various firing pattern models.
In particular: The parameters of the Izhikevich more described respectfully:• a[0.02]: recovery variable. An increment of a results in a quicker recovery for u.
• b[0.20]: underlines the strength of u to the sub-threshold fluctuations of the membranepotential v. An increment of b is translated into a stronger interrelation between u and v,
meaning possible sub-threshold oscillations and low-threshold spiking dynamics. u to thesubthreshold fluctuations of the membrane potential v.
• c[-65 mV ]: the after-spike reset value of the membrane potential v.
• d[6 mV ]: the after-spike reset of the recovery variable u.
When the membrane potential reaches a threshold value (typically around 30 mV), the neuron
produces an action potential or spike, and the membrane potential is reset to a reset value (typically -65 mV) while the recovery variable u is increased by a reset parameter (typically u += d).
The Izhikevich neuron model is known for its ability to replicate a wide range of spiking patterns
observed in real neurons, including regular spiking, bursting, and fast-spiking. By adjusting theparameters (a, b, c, d) of the model, various types of spiking behaviors can be simulated. Overall, the Izhikevich neuron model provides a computationally efficient approach to studying the
dynamics of spiking neurons and has been widely used in computational neuroscience and neuralnetwork modeling.
21.1.2

Network structure

In our network, we employ bidirectional communication between neurons, which is a communication pattern commonly observed in complex neural networks that are responsible for memoryand cognitive functions. In such networks, information flows in multiple directions, allowing
for the exchange of data and signals between neurons in both forward and backward directions.
This bidirectional flow of information enhances the network’s ability to perform sophisticatedcomputations and respond to a wide range of cognitive tasks [35]. As a result, a viable connection structure must be established for each node or neuron to ensure the proper functioning of themodel. When designing the dynamic network, one possible approach is to consider the von Neumann neighborhood type within an n-dimensional lattice, where every node is connected to its[2 × n] neighboring nodes. The configuration depicted in Figure 1.2 represents a 2-dimensionalgrid, where each neuron establishes connections with its adjacent nodes in the north, south, east,
and west directions. In the context of machine learning, the distance between nodes is quantified
using the Manhattan distance, also known as taxicab geometry. Unlike the Euclidean distance,
the Manhattan distance may yield multiple distinct paths with the same distance between two
points, as illustrated in Figure 1.1.

Figure 1.1: There are several shortest paths to go from A to B, but the Manhattan distance, which in thiscase is 12, is the same for any shortest path. Instead, in the case of Euclidean distance, there exists one
and only one shortest path from A to B

We generate, therefore, a matrix grid of neurons connected according to the von Neumann neighborhood model. In the two-dimensional case, the first row is up-connected with the last row (andvice-versa), and the first column is left-arrow associated with the last column and vice-versa, as
3

in figure 1.2. The edges of the two-dimensional matrix grid are theoretically bent into a threedimensional torus where a network representing a neural aggregate is arranged. This leads tosimilarities in toroidal neural networks (TNN) [36]. The toroidal structure opens the study ofseveral appealing geometric interpretations.

Figure 1.2: Toroidal network of 16 neurons plotted in 4 × 4 square matrices. To create a biologicallyplausible continuous system, that is to avoid the presence of boundaries or borders in the structure, thefirst row is connected to the last row and the first column to the last one.

For example, in Figure 1.3, the evolution of spike rates on a toroidal surface tends to approacha geodesic line, which is a local length-minimizing curve. From a topological perspective, thetorus exhibits significant homeomorphic properties, which pave the way for further explorationin computational topology to enhance our comprehension of brain function (Coli, Tozzi).
After framing the network within a toroidal model and defining neurons according to the Izhikevich archetype, we conducted simulations. In the initial phase of the program, we constructed amatrix where nodes represent neurons interconnected in a pattern to form a toroidal grid consisting of 400 neurons (20 × 20 arrays). To maintain simplicity, all neurons are considered identical,
devoid of inhibitory synapses and thalamic currents. Our primary focus lies in examining theimpact of non-uniform axonal propagation delay (Madadi), so we investigate the model underpristine conditions. Specifically, we focus on regular spiking (RS) neurons, as they typify themost common neuron class in the cortex. To characterize the delay and connection strengthfor each neuron, we employ the following terminology. Each neuron is defined by its constantsynaptic strength connection denoted as w, which remains consistent across the network andis not included in our calculations. The intrinsic axonal delay for each neuron is generally described by the variables µ = cd + δ, where cd represents the central delay, signifying the time
taken for a spike to propagate along the axon, and δ represents random variation around thisvalue (noise). Two main loops are involved in the process. The external loop iterates throughthe integer j, which controls the neuron’s central delay (cd), while the internal loop iterates
through i, which modulates the stochastic noise delay component (δ). We initially assign values
to the variables ns (representing noise steps) and (cd)max , which respectively denote the number
4

Figure 1.3: In A1 and B1 we run three simulations on the toroidal network of figure 1.2. Neuron 3 isconsidered the input of the network, being the only one connected to the external input and to initiallyspike. The three plots in A1 show the spike rate calculated in three different simulations of 1000 ms. In B1other three simulations are shown, but this time the input neuron is in position 1. On B1 the path 1-5-9-1314-15 shows higher spiking frequency on the nodes 1-5-9-13, compared to paths with the same Manhattan
distance as, for instance, the sequences 1-2-3-7-11-15 or 1-2-6-10-11-15. In A1, if we don’t consider the
direct connection 3 to 15, the shortest path from input(3) to output(15) is 3 steps. Nevertheless, a highspike frequency tends to lay on the path 3-7-11-15. Alternative paths with the same Manhattan distanceare 3-2-14-15 and 3-4-16-15. In the three-dimensional torus, if we consider the Euclidean distance, thesequences in red in A2 and B2 (1-5-9-13) tend to lay on the longitudinal geodesic.

5

of iterations in the internal i-loop and the external j-loop as follows:
1 ≤ i ≤ ns , and 1 ≤ j ≤ (cd)max , i, j ∈ N
When i reaches ns, the central delay (cd) increases by 1. For simplicity, we consider j = (cd)j toi,jassume the integers values (1, 2, 3, 4... (cd)max ). We indicate with µn,n
, the axonal transmissionkdelay for a spike to travel between two neurons n, nk . It is calculated at each simulation by afixed delay equivalent for each neuron, referred to as central delay (cdj ), and its stochastici,jcontribution δf,f
as follows:
p

z

i,j
δf,f

p

j

}|

{

(cd)
i,j
i,j
µf,f
= (cd)j + i
× αf,fpp
ns(1.1.2)

where
i,j
i,j
∗ αf,f= [2xf,f− 1]
pp∗ (cd)j = central delay; 1 ≤ (cd)j ≤ (cd)max , j ∈ N
j

i,j
i,j
i,j
≤ (cd)j= i (cd)such that: − (cd)j ≤ δf,f
∗ δf,f× αf,fpppnsi,jwhere 0 ≤ xf,f
≤ 1 is an aleatory variable with uniform distribution in [0, 1] assigned topi,j
the connection of neuron nf to nfp for each i-loop and j-loop, such that −1 ≤ αf,f
≤ 1.
pFigure 1.4: Biological representation of neuron

In Figure 1.4, on the right side, the biological representation of a neuron and its axonal delaytransmissions are denoted as µ. The central delay, labeled as cd, simulates the neuron’s axon,
while δ simulates the axon terminal. The total distance comprising these two segments varies
for each connection between nf and nfp , as shown on the left side of the figure. The simulation
starts with j=1 (cd=1). At each i-loop, each couple of neurons is assigned with an intrinsic
i,jdelay µf,fwhere, respectively, j represents the cd (same for each couple of neurons), and ipthe increment of axonal heterogeneity. This framework realizes a network of neurons having ai,j. In the whole network, only
different (heterogeneous) intrinsic axonal transmission delay µn,nk
nin is stimulated by an external current I = 10mA, as suggested by Izhikevich [29]. Once thisspikes regularly, its activity propagates to its neighbors with a certain delay, and so on. Hence,
we measure the propagation of spike activity between two defined neurons. The initiator nin andrandom neuron no . The subscripts in and o stand loosely for initiator and output respectively.
6

Consequently, we define:i,j• foutas the time of the first spike of no .

• fini,j as the time of the first spike of nin .
And, accordingly the spike propagation delay from nin to no as:i,j∆fgi,j = fout
− fini,j(1.1.3)In the context where g represents the minimum Manhattan distance (mMd) from nin to no , among
the two neurons [37] [36], a specific ∆fgi,j is given for each i. For every complete iteration of
the i-loop, a collection of ns values is assembled. All ∆fgi,j within one i-loop share the identicalcentral delay (cd)j . As i progresses to ns, the central delay (cd)j is increased by 1, initiatinga new cycle of i-loops. Upon j reaching (cd)max , the simulation concludes. The primary aim
is to assess whether the increment of axonal heterogeneity impacts spike signal transmission.
Specifically, the focus is on determining if as i increases it corresponds to a reduction in ∆fgi,j .
To achieve this, we calculate the regression line of the ∆fgi,j for each i-loop.

1.1.3

SGE: Stochastic Grid Enhancement

To elucidate the impact of axonal heterogeneity on spike propagation, our initial analysis involves tracking the temporal evolution of spike activity propagation within the network while
maintaining a constant delay value. In Figure 1.5, we set the delay parameter to (cd) = 50 milliseconds and calculate δ = (cd)/4 = 12 milliseconds. We monitor the ’time of the first spike’
for all neurons in the network, with only neuron nin exhibiting initial spiking activity. The propagation of the spiking signal commences from neuron nin and propagates through the network.
As anticipated, we observe a linear relationship between the Manhattan distance from the test
neuron and the initiator neuron. In our specific scenario, the network comprises a grid of 20 ×
20 neurons. We identify each neuron with an index based on its position within the grid. By
placing the initiator in the center of the second row, it is assigned an index of 30 (20 neurons in
the first row and an additional 10 neurons for the central position in the second row). Following
this indexing convention, we investigate the time of arrival of the spike-wave for each neuron
along the geodesic path toward the output neuron 230, which is also positioned on the geodesic
line in the central row. As evident from the simulation results depicted in Figure 1.5, the propagation delay in milliseconds exhibits a linear increase with distance, aligning with our intuitive
expectations.
We now proceed to assess whether the introduction of axonal heterogeneity has a positive impact
on the propagation of spiking activity throughout the neural network. In Figure 1.6, we present
two maps depicting the toroidal network under varying noise levels. Each pixel within thesemaps corresponds to an individual neuron, with its color indicating the time of the first spike inmilliseconds. Notably, neuron nin is situated at coordinates x=10, y=2 on the grid, indexed as
30. The remaining neurons exhibit spiking activity that follows the axonal propagation delay
7

Figure 1.5: The spike activity propagation delay on the torus. Here the first spike happens at time
t = 0 at the node located at index 30, the delay is calculated when the neuron at the node in the positionindicated on the horizontal axis has its first spike. The tested neurons are those located at the geodesic of
the toroidal grid. In these tests, the intrinsic delay and delay heterogeneity are kept constant at cd = 50
ms and δ = 12.5 ms respectively.

across the network. When we reduce the level of axonal heterogeneity, as illustrated in Figure
1.6 (A), we observe a corresponding decrease in the number of neurons reached by the spiking
activity, represented by deep purple regions on the map. Conversely, an increase in axonal
heterogeneity level, as depicted in Figure 1.6 (B), results in an expansion of the area covered byspiking activity. Our findings indicate that elevating the degree of axonal heterogeneity leads toa reduction in propagation delay. Furthermore, in the context of Equation (1.1.3), as the axonal
heterogeneity parameter (i → ns) is heightened, the value of ∆fgi,j tends to decrease.
To describe the effect of axonal heterogeneity from an analytical point of view, we reckon thegeneral equation for the regression line [38]:cx̄ + αbȳ = m

where

Pn
c=
m

σx,y
i=1 (xi − x̄)(yi − ȳ)
= 2
Pn
2
σx
i=1 (xi − x̄)(1.1.4)

Here we introducex̄ and ȳ as the average of the xi and yi .
σx and σy as the standard deviations ofx and y.
σx2 and σx,y as the variance and covariance.
j

We pose i (cd)
= (nd)i,j ; such that in our model (xi,j ; yi,j ) = ((nd)i,j ; ∆f i,j ).
ns
We now try to give a mathematical formulation for m and prove that since from our data m < 0,then we must have, in general, ∆f (i,j) < ∆f (e;j) when e ≤ i. In other words, when the noise
δi,j increases, the delay ∆fgi,j decreases. From (1.1.4) we calculate the slope for the linear
8

Figure 1.6: In this analysis, we present the temporal delay of spikes for each neuron within a 20x20
toroidal network. Specifically, we provide the timestamp at which the first spike occurs for every neuron
in the network, signifying the moment when each neuron first receives spiking information. The neuron
nin located at index 30 with coordinates x=10 and y=2 on the map, commences spiking at t=0 in response
to a constant input current of I0=10 mA. The color bar is utilized to represent the timing of the first spikefor each neuron in milliseconds. The areas depicted in dark blue on the map signify regions where neurons
fail to spike within the given time frame of our test simulation, which spans 500 ms. This occurs when
the spiking activity initiated by the ’initiator’ neuron nin does not propagate to these regions promptly.
Notably, in the first panel, where the noise level is minimal, a significant portion of the network remains
untouched by the spiking activity. In contrast, in the second panel, where the noise level is higher, thecirculation of spikes is accelerated, resulting in a greater number of neurons being influenced by the nin ’s
activity.

regression 1 , where:

Pns

i=1 (xi − x̄)(xi − x̄)

σx2 =

ns

Pns

i=1 (xi − x̄)(yi − ȳ)

σx,y =

ns(1.1.5)(1.1.6)our respective values are:(ns+1)
j; x̄ = (cd)j2ns• xji = (nd)i,j = i (cd)ns•yij = ∆f i,j

Pns

; ȳ =∆f i,jns

i=1

To simplify the symbolism from now on we reckon cdj = cd, and ∆f i, j = ∆f i . We omit theindex j since cdj is constant during the calculation of a single m while only i ranges from 1 tons. By substituting xi , x̄, yi , ȳ respectfully in (1.1.5) and (1.1.6), we get:
Pns

σx2 =

and
σx,y =12
i=1 (xi − x̄)ns=(cd)2 (ns2 − 1)
ns2
12

ns
ns
(cd) X(ns + 1) Xi[
∆f i ](i∆f)−
ns2 i=12
i=1

It is just reported a part of the demonstration.

9

which we can write as:
ns−1
2

2

σx,y = (cd)ns

X
i=1

|(i −

ns + 1)|(∆f (ns+1−i) − ∆f i )2

From (1.1.4) we can formulate:
P ns−12m=i=1|(i − (ns+1))|(∆f (ns+1−i) − ∆f i )2(cd) · (ns12−1)
2

We have m < 0 only if the numerator is negative, therefore, for most of the ∆f i , it must follow:∆f (ns+1−i) < ∆f i

(1.1.7).
The (1.1.7) demonstrated how, by increasing the index i (i.e., noise), the delay of a signal isreduced compared to the preceding one. As a result, we can conclude that increasing noise
reduces the delay of a signal from input to output for a fixed central delay (cd)j .
We varied the intrinsic axonal delay parameter (cd)j across a range of values, spanning from
1 to 100 milliseconds. For each of these values, we conducted experiments involving forty
different noise levels, ranging from (cd)j/40 to the maximum delay value (cd)max = 100ms.
The top panels of figures 1.7 and 1.8 illustrate the distribution of spike arrival times in the output
neuron no . In figure 1.7, with a relatively large central intrinsic delay of 93 ms, only 20 arrival
times are observed. Our simulation has a duration of 1000 ms, and depending on the stochastic
δi, j factor’s variability, in certain cases, neuron no does not spike, as the information fromneuron nin doesn’t arrive in time. Conversely, in figure 1.8, when the central delay is around
31 ms, the histogram is complete, displaying 40 arrival times concentrated at a shorter delay of
approximately 350 ms. The central panels of these figures depict the first spike time in no forincreasing values of (nd)i,j . The linear regression’s dark dotted line illustrates the significant
trend of time decreasing as noise levels increase. Notably, when examining the bottom panelsof both figures, the regression slope m consistently remains negative, and the minimum value
occurs around cd = 35 ms (with the y-axis reversed). The bottom panel of figure 1.7 presents
the values of m for each cd, spanning from the first to the last simulation (i.e., 2 ≤ cd ≤ 93.07).
We were able to replicate the minimum m value around cd = 35 ms in repeated simulations
with different random seeds. The maximum value of m consistently occurred within a global
minimum range of 30 ≤ cd ≤ 35.
The consistent negativity of the slope m, which signifies how axonal heterogeneity diminishes
the propagation speed within the network, persists across all conducted simulations. The replicable nature of this behavior underscores the robustness of this feature in our model neuralnetwork
10

Figure 1.7: In the top panel, we plot the histogram of the time difference between the first
spike of the output neuron no and the first spike of the initiator nin . These times represent the times of arrival of spiking activity from an initiator neuron placed at position 30 toa test or output neuron placed far away in the toroidal network (in position 230). We performed 40 simulations with increased noise levels, from cd/40 ms up to cd = 100 milliseconds. The resulting histogram is centered around a delay difference ∆fgi,j of about 800 milliseconds.
In the central panel, we plot the same arrival time differences in the function of the noise (nd)i,j .
We have 40 increasing values of (nd)i,j from 1 to (cd)j (cd is fixed in these plots to (cd)j =93.07). We notice that for a small value of heterogeneity ((nd)i,j between 0 and 23), the delay difference doesn’t change much, settling on values greater than 800. Eventually, by increasing the noise, ∆fgi,j decreases accordingly, reaching values below 800 ms. The dotted line represents linear interpolation with an inclination of m = −2.15, as indicated in the figure’s title.
In the bottom panel, we plot all the values of m obtained by previous simulations, where each has adifferent neuronal intrinsic delay (cd)j (ranging from 5 to 100 ms in 1 ms steps). While the top andcentral panels refer only to the value of (cd)j = 93.07, the bottom panel reports all the previous valuesof m up to the current (cd)j . A negative m indicates that the delay value is decreasing.

11

Figure 1.8: Same as figure 1.7. Precisely, we present the value of m corresponding to the global minimum for m=-4.17 (in the graph, minimum appears as a maximum as the ordinate are inverted for clarity).
The biggest effect of delay difference was observed in all of our testings for intrinsic delay values in therange 31 ≤ cd ≤ 35. When the intrinsic delay is centered at 35 ms, the signal from nin propagates faster.

1.1.4

Process time

The following section offers a theoretical interpretation of the negative aspects associated withthe variable denoted as m. Our investigation delves into the dynamic characteristics of neurons,
necessitating a comprehensive comprehension of the phase portrait of the membrane potential
and recovery system. The phase space features various conditions, including equilibrium, limit
cycles, attractors (as described by Izhikevich [29]), and other properties that aid in elucidating
the spiking behavior from a dynamic systems perspective. In this analysis, we examine two
distinct phases, commencing with the scenario in which a neuron resides in an equilibrium state,characterized by a resting membrane potential [39]. In Figure 1.9.1 of the phase space plotshown in Figure 1.9.2, input currents induce a minimal presynaptic potential (PSP) [40]. In thisphase diagram, the membrane potential exhibits slight deviations from the equilibrium point(represented as a black spot or an attractor corresponding to the resting potential). Subsequently,
after a few milliseconds, it returns to its resting state.
Instead, in figure 1.9 panel 3, we show the effect of two pre-synaptic stimuli. A little signal
PSP (A) causes a slight change in the equilibrium, whereas a bigger signal PSP (B) causes theneuron’s intrinsic dynamics to spike after a short time. The period between the arrival of thestimulus and the realization of the spike is referred to as process time. Now we’ll look at what
happens when a PSP signal A reaches a post-synaptic neuron while that same neuron is stillfiring as a result of a prior PSP signal B. Consider the case where a spike leaves a neuron nf
and travels in the direction of neuron nf 1 .
Firstly we define:12

Figure 1.9: Panels 1 and 2: a stable equilibrium. A small pre-synaptic signal (PSP) goes from neuron
nf to nf1 . The membrane potential will fluctuate around the equilibrium point before returning to its
resting state (black dot, the equilibrium point). Panels 3 and 4 depict the responses of two different presynaptic stimuli. The small PSP (A) doesn’t induce a spike, whereas a bigger PSP (B) does. The time ittakes for the neuron to achieve a spike following the PSP is referred to as process time, and it is measured
in milliseconds as in the Izhikevich regular spiking neurons employed in our model.

• Process time τfi,j: that represents the time for neuron nf to process a pre-synaptic spikep
and deliver its post-synaptic spike, the two indexes i and j represent the current central
delay (cd)j and its stochastic variation δ i,j .
• Axonal transmission delay µi,jf,fp : the time for a spike to go from neuron nf to one of its
four adjacent neuron nf p,(p=1,2,3,4) , as in eq.(1.1.2).
i,j:• Total transmission Delay df,f
p

i,j
i,j
df,f
= µi,jf,fp + τfp
p(1.1.8)i,j, characterizes the delay transmission from nf to nf p , and a process
The random variable xf,f
p
i,j
i,j
. Once everything is in place, the network
time τfp . Therefore, each i-simulation defines df,fp
operates as expected: a spike reaches a specific input neuron nf . The neuron nf will process the signal with a processing time τfp . Subsequently, the spike leaves nf and reaches four
adjacent neurons nf p,(p=1,2,3,4) with four different axonal transmission delays µi,f, jfp defined inequation(1.1.2). We now want to focus on the case when :i,jµi,jf,fp < τfp(1.1.9)Following equation (1.1.9), the spike generated by neuron nf is anticipated to reach neuron nf1
promptly, even as the post-synaptic neuron is still in the process of handling the signal from
B ′ . The amalgamation of these two membrane potentials, as illustrated in the phase diagram
depicted in Figure 1.9.4, serves to expedite the generation of a spike, thereby reducing the processing duration denoted as τfi,jp . Consequently, the overall time interval spanning from the13

Figure 1.10: A sketch representing how the variability of the axonal delay influences the overall spike
activity propagation along with the network. Two cases are shown: in panel 1. the variability of theaxonal transmission µi,jf,fp is wide (from zero to 50ms in this example where the central value of the delay
is 25ms). Since this value can become smaller than the processing time τfi,j, this will positively influencep
the phase space of the neuron’s dynamic, resulting in a faster post-synaptic response. In panel 2 it isshown the opposite situation in which the variability of the axonal delay is smaller and never realizes
this condition. Since a longer than average, µi,jf,fp doesn’t increase post-synaptic process times, whereasa shorter than average does, this indicates how a degree of randomness helps the overall circulation of
spike activity on the toroidal network.

initial pre-synaptic stimulation initiated by A′ to the eventual post-synaptic spike, labeled asd;i,jf,f1 , undergoes a reduction. In simpler terms, owing to the stochastic variability inherent inaxonal transmission delays denoted by µi,j
f,f p , the temporal delay may fluctuate randomly, either exceeding or falling short of the mean value. In instances of shorter delays, the arrival ofa second spike at the post-synaptic neuron leads to decreased processing times. Conversely,i,jwhen µi,j
f,f p is higher, the processing times τfp remain unaffected. This inherent asymmetry inthe system accounts for the observed phenomenon in our simulations. As the degree of axonal
heterogeneity in the network increases, the likelihood of this condition occurring across multiple
neurons rises, resulting in accelerated propagation of spiking activity and enhancing the overall
flow of information throughout the network. For a visual representation of this phenomenon,
please refer to Figure 1.10.
In Figure 1.11, we provide a schematic illustration of the operational principles of the SGE
effect concerning the temporal aspects of processing time τ depicted in red, and transmission
time µ represented by the blue line. Our theoretical framework postulates that the presenceof processing time τ is a prerequisite for the occurrence of the SGE effect. In this context,transmission time µ, given a fixed central delay, signifies the duration required for an impulse
to traverse the spatial gap between two neurons, whereas τ denotes the interval during whicha neuron emits a post-synaptic spike after the arrival of pre-synaptic stimulation. The very
existence of this mechanistic framework leads us to propose that during the τ interval, while aprior signal is undergoing processing, as denoted by (PSPs)B in the nf region of the plot, there
exists the potential for another spike to propagate from the same source, nf, thereby expediting
the emission process and ultimately facilitating the manifestation of the ”SGE” effect.

1.1.5

SNE: Stochastic Neuron Enhancement

The initial type of heterogeneity manifests as non-uniform spatial gaps between neurons, achieved
through the introduction of random axonal delays. This results in accelerated information trans14

Figure 1.11: Representing how the variability of the axonal delay µ, triggers the SGE effect influencing
the overall network’s spiking activity. Two cases are shown: Process time > Transmission time (panel
1.,3.) and Transmission time > Process time (panel 2.,4.). In panel 3 we see how increasing the randomness results in a wider range of variability of µi,j
f,f p from 0ms to 50ms, with a central value of 25ms. In
this case, µ can become smaller than the processing time τfi,jp , and positively influence the phase spaceof the neuron’s dynamic, provoking a faster postsynaptic response. In panel 4, for a small randomnesslevel, the variability of µ never gets below the value of τ and the SGE effect won’t take action.

mission, denoted as SGE effect. The second type of heterogeneity revolves around a neuralmodel that lacks uniformity. In this scenario, we manipulate neuron cell types by adjusting theparameters of the Izhikevich model, thereby altering various differential characteristics of theneurons and influencing their spike patterns and dynamics (see Figure 1.12). The simulations
demonstrate that heterogeneity amplifies the overall number of spikes throughout the network,suggesting an augmentation in information robustness. Our findings align with recent notablestudies, indicating that heterogeneity confers various advantages to the brain, including improved coding efficiency, enhanced reliability, superior working memory, and other functional
attributes [20]. Henceforth, we shall refer to this phenomenon, stemming from the stochastic
nature of the neuronal model, as ”Stochastic Neuron Enhancement” (SNE effect).

1.2

Simulation on multidimensional networks

We will now explore the impact of the combination of SGE-SNE effect on various lattice grids,
spanning dimensions two, three, and four, along with their associated toroidal frameworks. For
instance, when considering a 3D lattice grid, we observe the emergence of a 4D toroidal network,where each neuron is surrounded by six adjacent neighbors, as opposed to four. To achieve
the 4D torus, we achieve this by connecting opposite surfaces of the cube, creating a toroidal
network that exists in a four-dimensional Clifford torus [41].
By employing a comparable idea, we can expand this to encompass a four-dimensional grid,15

Figure 1.12: To different values of the parameters, a, b, c, d, correspond different types of Izhikevich
neuron models as described by equations (1.1.1). In our simulation, the constant values of a, b, c, d=(0.02,
0.2, -65, 8) are congruent to a regular spiking (RS) model, while random changes of c, d correspond toa continuous variability of a large range of neuron models [29] such as IB (intrinsically bursting), CH(chattering), FS (fast-spiking), LTS (low-threshold spiking), TC (thalamocortical), RZ (resonator). The
variability of c, d according to eq.(1.2.1), guarantees the neural heterogeneity of the model.

resulting in the formation of a torus situated within a five-dimensional space. In the following
illustration, we present an overview of the connection between the dimensional grid and thecorresponding toroidal network.
• 2D Lattice grid =⇒ Torus in 3D space• 3D Lattice grid =⇒ Torus in 4D space• 4D Lattice grid =⇒ Torus in 5D space
In the context of the three-dimensional topological arrangement [42], each node establishes
connections with its six adjacent neighbors, enabling communication in six distinct directions:+x (east), -x (west), +y (forward), -y (back), +z (north), and -z (south). The four-dimensional
torus structure naturally results in a von Neumann-shaped neighborhood [43], with cells being
surrounded by six input connections [22]. In both scenarios, we introduce an external current of
I=10 mA to an arbitrary neuron denoted as nin , following the methodology presented in section(1.1). This current stimulates the neurons within the network, initiating the SGE effect. We
subsequently examine the occurrence of this phenomenon on 2D, 3D, and 4D grids, investigate
the impact of neural heterogeneity, and explore the influence of SNE. Finally, we consider the
combined effects of both types of heterogeneities on the network.
In our model, we will keep a, b at constant values of 0.02 and 0.2, while the parameters c, d areheterogeneous, as in figure 1.12, will be subject to the form:c = −65 + (15 × x1 × x1 )d = 8 − (6 × x2 × x2 )0 ≤ x1 , x2 ≤ 1: uniformly distributed random variables [29].
16(1.2.1)

Figure 1.13: We present four distinct depictions of a 3D Torus derived from a 2D grid model. In thisi,jcontext, the network’s neurons serve as nodes, and their edges are the axonal transmission delay µf,fpdefined in eq.(1.1.2). The x-axis signifies the extent of neural diversity denoted as H, which varies from 0
to 1. When H is at 0, all neurons are identical to each other, exhibiting neural homogeneity. Meanwhile,the y-axis illustrates the range of axonal heterogeneity, specifically the axonal delay δ i , which spansi,j.
from 0 to the central delay value cd. As δ i varies, it results in differing values for the connections µf,fp
A: In the specific scenario where H equals 0 and all δ i are 0, the network grid exhibits regularity. As a result, all neurons are identical, denoted by the presence of red nodes, specifically of the R.S.
Izhikevich type, and the axonal delay between them is uniform at the value of cd within a regular 3D
layout.
B:Hypothesis H states that δ i is equal to 0. Our approach focuses on altering neural heterogeneity
exclusively by adjusting the parameters c and d in the Izhikevich neuron model. As depicted in Figure
1.12, this manipulation leads to the creation of a stochastic typology of neuron models. Importantly, theinter-neuron spacing remains fixed at the value of cd, while the varying colors of each neuron represent
distinct Izhikevich neuron model types.
C: Under the conditions where H equals 0 and δ i equals 1, all the neurons in the system are modeled
as R.S. Izhikevich neurons. The variability in axonal delay between neurons is induced by axonal
heterogeneity, as described by equation (1.1.2).
D: H: 1, and δi : 1. The illustration represents the utmost variability in structure, showcasing two formsof heterogeneity. It provides a visual representation of the most likely distribution of neuron populations
in a network of diverse cells, where the inter-neuronal connection distances vary.

17

Heterogeneity of Izhikevich neuron modelParameters
a:
b:
c:
d:

RS
0.02
0.2
-65
8

IB
0.02
0.2
-55
4

CH
0.02
0.2
-50
2

FS
0.1
0.2
-65
2

LTS
0.02
0.25
-65
2

TC
0.02
0.25
-65
0.05

RZ
0.1
0.26
-65
2

Table 1.1: Values of the parameters a, b, c, d used for each model of figure 1.12.

In Figure 1.13 we show the effects of the combined heterogeneity on a 3D toroidal model network. Where the node’s color variation represents the neural heterogeneity, the different distances of the connections illustrate the axonal delay heterogeneity.
We established three distinct networks in two, three, and four dimensions, each corresponding toa three-dimensional, four-dimensional, and five-dimensional torus, respectively. Subsequently,
we conducted comparative simulations to scrutinize the outcomes. Initially, we constructed athree-dimensional matrix composed of 343 neurons, arranged in a (7 × 7 × 7) grid, interconnected using the von Neumann neighborhood model, creating a toroidal structure by connecting
opposing surfaces [36]. Once the network framework was established, we configured the neurons according to Izhikevich’s model to initiate the simulations. Our primary focus in this study
is to investigate the impact of non-uniform axonal propagation delay [44]. To begin, we examined the model in its simplest form, assuming that all regular spiking (RS) neurons wereidentical, without considering inhibitory synapses or thalamic currents. Subsequently, we willcompare these results with those obtained from a model where neurons of the same size exhibit
heterogeneity, and we will analyze the variations in the effects of heterogeneity.

1.2.1

Grid dimensions and SGE effectWe will now extend our analysis of the SGE effect phenomenon to various dimensions. We
illustrate how this phenomenon becomes more pronounced when we increase the network’s
dimensionality. To begin, we explore the number of minimum paths with the same mMd between
two random neurons in a lattice grid configured in a toroidal shape. The mMd from neuron nin
to neuron no represents the minimum number of steps required to reach no from ni n, as depicted
in Figure 1.14.
In a previous study involving a 2D lattice network (as described in [45]), each neuron had four
adjacent neighbors. For any pair of neurons [nin ; no ], there were a specific number of paths
with the same mMd that connected nin to no . Given the network’s structure, the most efficient
way for a signal to reach the no was in mMd steps. In two dimensions, there were only two
possible pre-synaptic connections to no , as shown in Figure 1.15(a). Initiating a path from oneof the other two connections would require more steps than mMd. The strength of the SGE
effect depends on the number of paths with the same mMd, which increases by manipulating thenumber of dimensions of the network and the nodes’ coordinates as shown in the following.
18

Figure 1.14: In this case, the mMd from nin to no is 4. The red dotted lines represent all the possible
paths with mMd, where all the others (as the blue dotted line) have a Manhattan distance > mMd. In a
3D lattice grid, each node/neuron has six connections. The nodes on the border surface are connectedwith their respective counterparts on the opposite face. For clarity, the brown dotted line represents oneof those cyclic connections. This structure is equivalent to a torus in 4D space.

Mathematical proof
Let’s consider the following 3D cubic lattice (4D torus):• side: ln = 2 × sl
• coordinates of [(nin , no )]: [(x0 , y0 , z0 ); (x1 , y1 , z1 )]We can formulate :
X=

Y =

Z=



|x0 − x1 |; if : |x0 − x1 |≤ slln − |x0 − x1 |; if : |x0 − x1 |> sl
|y0 − y1 |; if : |y0 − y1 |≤ slln − |y0 − y1 |; if : |y0 − y1 |> sl
|z0 − z1 |; if : |z0 − z1 |≤ slln − |z0 − z1 |; if : |z0 − z1 |> sl(1.2.2)(1.2.3)(1.2.4)In this context, the variables X, Y, and Z represent the minimum number of steps required alongeach dimensional axis to go from nin to no . This definition aligns with the concept of the minimum Manhattan distance, which can be expressed as mMd = X + Y + Z. For instance, suppose
our calculations yield the values (X, Y, Z) = (3, 4, 5). Consequently, the minimum number of19

Figure 1.15: Illustration of a network in both 2D and 3D is presented. When we contrast these two network representations, it becomes evident that as we elevate the dimension, the quantity of pathways withidentical mMd also increases proportionally. Moreover, the presynaptic activity experiences an uptick in
the higher-dimensional network, consequently amplifying the likelihood of observing the SGE effect.
In panel (a), within the two-dimensional grid, there exist three distinct routes from nin to no , each with a
mMd value of 3. Under these circumstances, there is a possibility of having up to two presynaptic signals,
denoted as s1 and s2 , arriving at no almost simultaneously, thus instigating the SGE effect.
In panel (b), the three-dimensional grid creates six distinct routes, all with an identical mMd value of 3,
connecting nin to no . In this scenario, the presence of a greater number of routes allows for the simultaneous convergence of up to three presynaptic signals, denoted as s1 , ; s2 , ; s3 , at no , thereby enhancing
the likelihood of the SGE effect becoming significant.

20

steps to traverse from nin to no is mMd = 3 + 4 + 5 = 12. To elaborate, one could start at nin ,take 3 steps along the x-axis, followed by 4 steps on they-axis, 5 steps on the z-axis, and thenconclude with 2 more steps on the y-axis. This path from nin to no illustrates an mMd of 12.
Our objective is to determine the total number of paths with an mMd of 12. Intuitively, this canbe achieved by permuting the partitions of X, Y, and Z, and the precise solution can be obtainedusing the multinomial coefficient [46]:P =(X + Y + Z)!
X! Y ! Z!

(1.2.5)

P represents here the number of paths with the same mMd=X+Y+Z from nin to no in the 3D
lattice grid. By keeping nin in the same position, if we increase the number of dimensions of
the grid, we may expect the number of paths with the same mMd to increase accordingly (as infigure 1.15). This is not always true; instead, it depends on the choice of no . Let’s simplify theconcept by analyzing the number of paths from nin to no with the same mMd from a 2D to a 3D
lattice grid. In the 2D case we pose mMd = X ′ + Y ′ , in the 3D case mMd = X + Y + Z. Since
we want mMd to be the same, it must be X ′ + Y ′ = X + Y + Z. In a 2D grid (10 × 10) let be(nin , no ) of coordinates [(1,1) ; (3,4)], it follows X = 3 - 1 = 2, Y = 4 - 1 = 3, and mMd= 3 + 2 =5. By applying the eq.(1.2.5) in 2D, we can calculate the number of path with mMd=5 as:P′ =5!
(X ′ + Y ′ )!
=
= 10
′
′
X !Y !
12

In the 3D grid (10 × 10 × 10) we assign to nin the coordinates (1,1,0). Depending on the choiceof no , the number of paths with the same mMd=5 will change. For no of coordinates (2,2,3), in
fact, we have X,Y,Z=(1,1,3), mMd=1+1+3=5, and the number of paths is:P =(X + Y + Z)!
5!
=
= 20
X! Y ! Z!
1! 1! 3!

In this case, the number of paths in 3D increased consistently, but if we had chosen no of coordinates (0,1,4), we would have:P =5!
(X + Y + Z)!
==5
X! Y ! Z!
0! 1! 4!

which shows fewer paths compared to the 2D case (P ′ =10). Therefore, the choice of no would
affect the number of paths which influences the SGE effect. We propose, hence, a way to determine the coordinates of no such as to benefit the SGE. Let’s define the ratio α=P/P ′ , andconsider X ′ + Y ′ = X + Y + Z = mM d. We want to find the coordinates of no such thatα > 1, which means P > P ′ , i.e., the number of paths in 3D are greater than those in 2D withthe same mMd. Since X + Y + Z = X ′ + Y ′ = mM d we can write:α=(X ′ ! Y ′ ! )(X ′ ! Y ′ ! )P(X + Y + Z)!
×==P′
X! Y ! Z!
(X ′ + Y ′ )!
X! Y ! Z!
21

For α > 1:

(X ′ ! Y ′ ! )(X ′ ! Y ′ ! )> 1 → Z! <X! Y ! Z!
X! Y !

It follows:
X + Y + Z = mM d


(1.2.6)
Z! < (X ′ !Y ′ !)X!Y !

Let’s clarify everything with an example by using the same coordinates in 2D for (nin , no ) =[(1,1) ; (3,4)]. We know mMd=5, P ′ =10, X ′ =2, and Y ′ =3. We now choose some values for
X, Y and calculate Z since we know X+Y+Z=mMd. Consequently, we just need to verify if Z
satisfies the eq.(1.2.6). For instance, if we chose the values for (X,Y )=(1,2), we have Z=2, which
satisfies eq.(1.2.6):

X + Y + Z = mM d




→





Z! < (X ′ !Y ′ !)
X!Y !




1+2+2=5







2! < (2′ !3′ !)1!2!

Therefore (X,Y,Z)=(1,2,2) solves for α > 1, and by substituting (1,2,2) in eq.(1.2.5), we expect
to have P > P ′ :5!
(X + Y + Z)!
== 15 > P ′ = 10
P =X! Y ! Z!
1! 2! 2!
From the eq.(1.2.2, 1.2.3, 1.2.4), by knowing [nin :(1,1,0)] and [(X,Y,Z)=(1,2,2)] we can calculate 2 points for no =[no,1 : (0,3,2); no,2 :(2,3,2)]. The points have both mMd=5 to nin :(1,1,0)
whose number of paths is P = 15, greater than the number of paths P ′ =10 for 2 points having
mMd=5 in the 2D network. Hence, by carefully choosing the coordinates of the no (no,1 or no,2 ),we can strengthen the SGE in the 3D network.
Let’s evaluate two cases.
• A: All distances are equals X=Y=Z, it follows:
P =

(3X)!
>1
(X! )3Since the fraction is bigger than 1, for n → ∞, considering the distances all equals, wehave:(n × X)!
P = limn→∞→∞(X! )n.
• B: In the 3D case where the distances are not equal, we can assume there is a minimum
between them. For instance, let’s be X < Y < Z (the order is not relevant). We candevelop eq.(1.2.5) by expressing Y, Z in the function of X. Since X, Y, Z > 0 we canwrite Y = X + A and Z = X + B, with 0 < A < B. We can write eq.(1.2.5) as:
22

P =

[(X + B) + (X + A) + X]!
[3X + B + A]!
=(X + B)! (X + A)! X!
(X + B)! (X + A)! X!

The numerator can be developed in 3 terms : N = S × L × R :S = (3X + B + A) · ·(3X + A + 1) ; L = (3X + A) · ·(3X + 1) ; R = (3X · ·X! )
|

{z

}

|

{z

Bterms}|Aterms{z}

3Xterms

The denominator can be developed in 3 terms as: D = S’ × L’ × R’
S ′ = (X + B) · ·(X + 1) ; L′ = (X + A) · ·(X + 1) ; R′ = (X! X! X! )
{z

|

}

|

Bterms{z}

Aterms|{z}

3Xterms

Let’s consider the fraction N/D by comparing the fractions of the B, A, 3Xterms .
Bterms =S
LR; Aterms = ′ ; 3Xterms = ′′SLRAll the 3 fractions are greater than 1. If we increase the dimension n, P can be written as
sum of n fractions all greater than 1, and for n → ∞:P = limn→∞

n
X
Ti
′
i=1 Ti→∞

In any case, by increasing the dimension n, P increases.

1.2.2

Spiking activity

Single spike propagation
To gain a deeper understanding of the process time τ , we conducted a simulation in a 3D latticegrid network consisting of (7 × 7 × 7) identical RS neurons. The graphical representation in
Figure 1.16 illustrates a simulation spanning 1000 milliseconds during which a solitary spiketraverses the quickest route from the input neuron nin :12 to the output neuron no :155 (witheach neuron’s number corresponding to its position in the grid). In this simulation, an initiatingneuron emits a single spike at time t = t0 , which in turn stimulates its neighboring neurons.
This action generates post-synaptic spikes that transform into pre-synaptic signals transmitted
to adjacent neurons, and the process continues to propagate across the grid.
By categorizing all impulses with positive membrane potential in Figure 1.16 as ”spikes,” weobserve that the second spike (occurring at the first post-synaptic neuron 61) occurs approximately at 100 milliseconds. Subsequently, the adjacent neuron 154 spikes around 110 milliseconds. The initial spike does not propagate entirely to the six adjacent neurons; instead, its
potential fragments into smaller post-synaptic potentials (PSPs) that spread along the interconnected neurons in the system. A single neuron can be influenced by up to six PSPs, each arriving
at approximately the same time, and eventually trigger a spike (as indicated by the green dot).
23

Figure 1.16: The membrane potential (vertical axis) of 6 adjacent neurons representing the propagationof a single spike in a 2D network of regular spiking (RS) Izhikevich neurons. The central delay is 10 witha noise of 6.43 ms, meaning the transition time µ varies in the range (13.57;26.43) ms. The plot shows
the behavior and development of the single-impulse spike (green dot) from the initiator neuron nin =12
across the path [12,61,62,111,105,154,155], with 155 being the output neuron no . We can observe how
the initial spiking activity at neuron 12 induces connected neurons to spike after a brief delay of durationnear the average axonal delay value (cd=20 milliseconds in the simulation shown).

24Remarkably, a single impulse spike has the capacity to propagate and generate significant action
potentials throughout the entire network within a 1000-millisecond timeframe.

RS vs Heterogeneous model
In this study, we investigate and compare the correlation between the randomness of neuron
types and the average number of spikes that occur in the entire network over a specific duration.
Initially, we examine homogeneous and heterogeneous models using a 3D model (4D torus).
We conduct multiple simulations while varying the level of heterogeneity from 0 to 1. A heterogeneity level of 0 signifies that all neurons are of the same type (RS Izhikevich), while a levelof 1 represents a completely heterogeneous network of Izhikevich neurons, as defined in reference [29]. In all these simulations, a constant current of 10 mA triggers the input (initiator) withnin set to 12. It’s important to note that increasing axonal delay and heterogeneity accelerates
the propagation of spikes, leading to a logical expectation that the number of spikes will riseas heterogeneity increases. In contrast, an examination of Figure 1.17 uncovers an unexpectedphenomenon, where the occurrence of spikes decreases as the level of variability increases. The
plot depicts the results of 15 different simulations, each lasting 1000 ms, conducted at six different levels of heterogeneity. In each simulation, the central delay is held constant at 22 ms, whileδ i ranges from 0 to 20 ms. The first simulation provides the average number of spikes in the
network with a central delay of 22 ms and a maximum variability of 0 ms. In this scenario, allneurons exhibit the same transition time, with µ equal to 22 ms. Furthermore, since all neurons
are of the RS type, their processing time for action potential signals, denoted as τ , is approximately uniform. In Figure 1.16, assuming constant values for µ and τ , we observe a smoothsynchronization of signals in the network. This synchronization implies that signals, such ass1 , s2 , and s3 in Figure 1.15(b), are likely to reach no simultaneously, given the constancy of µand τ . This synchronization results in the PSPs clustering together at approximately the sametime, leading to the generation of spikes.
On the contrary, when the level of stochastic variation in axonal delay amplifies, the transmissiontime, denoted as µ, exhibits fluctuations as described in equation (1.1.2). As a result, the synchrony of signals diminishes, causing them to reach their destination with varying time intervals,
thereby diminishing the probability of postsynaptic potentials (PSPs) converging to generate aneuronal spike. As the degree of noise escalates, the level of synchronization decreases, whereas
the incidence of spike occurrences rises.
Surprisingly, it is evident that as axonal heterogeneity value changes, the spiking activity of the
network maintains a consistent shape (see Figure 1.17) for all values of neural heterogeneity.
To elaborate, while the heterogeneous axonal delay primarily influences the spatial separation
between neurons, neural heterogeneity is fundamentally associated with the neuron model, affecting the variability in τ specifically, how distinct neurons process pre-synaptic signals withvarying process time.
25

Figure 1.17: The study investigates the relationship between axonal delay randomness, which ranges
from 0 to 20 milliseconds, and the number of spikes. This analysis focuses on a central delay of 22
milliseconds and encompasses six different levels of neuronal heterogeneity across three distinct network
dimensions. The count of spikes and the associated standard deviation are computed across hundreds ofneurons, and the outcomes are influenced by the network’s size. The graph illustrates the results, with ablue dot (1) representing complete heterogeneity and a yellow dot (0.0) symbolizing zero heterogeneity(all RS neurons). Across all three networks, the number of spikes for H=1.0 is approximately twice
as high as that for H=0.0. It’s worth noting that as spatial variability increases, the number of spikes
decreases, and this effect is more pronounced in models with higher levels of heterogeneity. The righthand panels display the linear regression coefficients (top) for the six heterogeneity levels, ranging from
0.0 to 1.0. The central values show the mean number of spikes for each level of heterogeneity (center)and the corresponding standard deviation (bottom).

26

Nt.size
5×5×5×5

S.D.
Mean
Slope

H.:0.0
55.3
3091
-201

H.:0.2
55.3
3101
-190

H.:0.4
55.7
3340
-228

H.:0.6
56.6
3636
-231

H.:0.8
57.7
4543
-345

H.: 1.0
47.7
6166
-392

7×7×7

S.D.
Mean
Slope

55.7
1104
-81.2

51.1
1069
-72.7

53.5
1133
-80.0

55.6
1228
-91.3

57.3
1445
-113

52
1978
-144

11 × 11

S.D.
Mean
Slope

28.4
223
-8.02

28.2
210
-7.25

36.4
226
-10.6

34.5
243
-10.9

38.6
265
-12.9

42.7
311
-18.3

Table 1.2: The table summarizes the data of figure 1.17. For each of the three networks are specified
the value of standard deviation (S.D.), the average number of spikes of the network (Mean), and theinterpolation line (Slope) for all the heterogeneity values (from H:0.0 to H:1.0)

1.3

Network partitioned in equivalence classes.

1.3.1

mMd induces equivalence classes on the network: Octahedrons.

We commence by partitioning the network based on the equivalence relationship denoted as
mMd. This process results in the formation of neuron groups that share similarities in their connections and characteristics, facilitating the subsequent examination of the network by scrutinizing these groups. To begin, we revisit the definitions of equivalence classes and the equivalence
relationship, as presented by ϕ [47] within the context of a given set F. A relation ϕ defined ona set F is classified as an equivalence relation if and only if it satisfies three essential criteria:
1. It is reflexive: for any a ∈ F , a must be equivalent to itself: a ϕ a
2. It is symmetric: For any couple (a,b) ∈ F , if a ϕ b → b ϕ a
3. It is transitive: Let a, b, and c ∈ F .
Then, if a ϕ b and b ϕ c → a ϕ c
If ϕ is an equivalence relationship on F, we can define the equivalence classes on F as:[a] = {x ∈ F | x ϕ a}
It can be shown that any two equivalence classes are either equal or disjoint, hence the collectionof equivalence classes forms a partition of F, called quotient set of F, indicated as F , divided
into non-empty subsets, in such a way that every element of F is included in exactly one subset
of F [48]. The quotient set is defined as:
F/ϕ = F = {[a] | a ∈ F }(1.3.1)Hence, we build the same mathematical structure in the network by proving mMd to be an equiv27

alent relationship in Nt. This induces a quotient set of equivalence classes on it. We randomly
choose a neuron nin as the initiator. We consider the neurons ni , nj to be in relationship ϕ,
ni ϕ nj , if they have the same mMd to nin . We refer to mMd for nin , as the minimum number of steps to reach nin from a node nj , as explained in figure (1.18). In order for ϕ to be an
equivalence relationship on Nt we need to prove the following three points:• reflexive: ∀ ni ∈ N t → ni ϕ ni .
The reflexivity here is trivial since ϕ is an equality relationship. If ni has mMd to nin it isobviously in a relationship with itself.
• symmetric: ∀ ni , nj ∈ N t if ni ϕ nj → ni ϕ nj
For any couple of neurons ni , nj ∈ N t, if ni ϕ nj , it means ni , nj have the same mMd tonin for the commutative property of equality relationship it follows → nj ϕ ni .
• transitive: ∀ ni , nj , nk ∈ N t :if ni ϕ nj and nj ϕ nk → ni ϕ nk .
If ni ϕ nj → ni and nj have the same mMd to nin . If and nj ϕ nk → nj and nk have thesame mMd to nin . From the transitivity properties of equality, it simply follows ni and nk
have the same mMd to nin , that is to say, ni ϕ nk .
Therefore, given a random nin , ϕ is an equivalence relationship on Nt inducing equivalenceclasses of the form:[nmin ] = {ni , nj ∈ N t | ni ϕ nj }, with: 1 ≤ m ≤ g(1.3.2)where g is the maximum mMd from nin . There would be, hence, g equivalence classes describedin eq.(1.3.2).
In Figure 1.18, the neurons belonging to a class denoted as [nmin ] are positioned along a geometricstructure that can be described as an ”octahedron.” These classes take their names from theshape itself. The class [nmin ] encompasses all neurons nj that share the same mMd value as
the initiating neuron nin . The quotient set N t/ϕ serves as a valid partition of the set N t andrepresents the collection of all octahedrons. These octahedrons possess the properties of having
non-overlapping intersections (meaning that a neuron n cannot be part of more than one distinct
octahedron) and collectively forming the entirety of the set N t.

1.3.2

Statistical entropy: macrostate and microstate.

Statistical entropy holds significant relevance within the domain of statistical mechanics. As
expounded by Boltzmann’s research [49], entropy can be construed as a quantification of themultitude of conceivable microscopic states (microstates) from which a given macrostate of asystem may emerge [50]. If we denote Ω as the count of microstates contributing to the realization of a specific macrostate denoted as A, then, in accordance with Boltzmann’s equation, thelikelihood of this particular macrostate manifesting is directly proportional to the value of Ω:
28

Figure 1.18: An illustration featuring the mMd and the corresponding equivalence classes [n1in ] and[n2in ] is presented. The central green node represents the initiator neuron, denoted as nin . In Figure A,
we depict the class [n1in ], comprising a collection of six neurons, all with a mMd value of 1 with respect
to nin . The mMd is represented by a red dotted line. In Figure B, we turn our attention to the class[n2in ], which consists of eighteen neurons, each having an identical mMd of 2, measured from the samereference neuron, nin . The red dotted lines in this context illustrate the paths required to reach nin inexactly two steps. These sets are comprised of 18 elements, all maintaining equivalent relationships witheach other. Hypothetically connecting these nodes would result in the neurons lying on the surface of anoctahedron.

S(A) = ln2 (Ω)(1.3.3)We have chosen to omit the inclusion of Boltzmann’s constant in the current context, as it proves
unnecessary for the purposes of our study; nonetheless, it has been computed in the course of oursimulation. The function under consideration quantifies the total number of feasible arrangements for a given macrostate. The quantity denoted as S(A), corresponding to the entropy ofmacrostate A, serves as a measure of the information required for describing the system. The
probability of encountering the system in a particular state, or macrostate, hinges upon the microstates that compose that specific macrostate; in other words, it is directly proportional to thecount of possible microstate configurations. To elucidate this relationship among microstates,
macrostates, and entropy, we offer an illustrative example related to equation (1.3.3). Whenrolling two dice, the microstate associated with macrostate A=2 is singular, while the number ofmicrostates for macrostate A=7 amounts to six, reflecting the six distinct combinations of dicethat result in a total of seven pips. As depicted in Figure (1.19), the macrostate with the highest
likelihood of occurrence, or, equivalently, the highest entropy, is macrostate A=7, characterized
by an entropy of S(7) = ln(6). This observation aligns with the fundamental principle that themacrostate possessing the greatest number of microstates corresponds to the macrostate with thehighest entropy
29

Figure 1.19: In panel A, the system can assume one of the eleven possible macrostates from 1 to 12 on
the x-axes. The macrostate 7 is the one with the highest probability to happen since it has the greatest
number of microstates Ω = 6, y-axes. Therefore, it is the macrostate with the highest entropy and themost likely state the network will tend to.
In panel B, the picture represents the number of microstates for the macrostate of a network of (11×
11 × 11) neurons. For the initiator neuron n12 , mMd assumes values from 1 to 15, formulating fifteen
octahedrons and, therefore, fifteen macrostates. According to eq. (1.3.3) the macrostate with the highest
8 with 180 microstates whose entropy is S(I 8 ) = log(180).
entropy is the macrostate I12
12

1.3.3

Macrostate and microstate of octahedrons: interspike train.

In the previous sections, denoted as 1.3.1 and 1.3.2, we have provided definitions for octahedrons as equivalence classes within the network and statistical entropy in terms of microstateand macrostate. Our objective now is to incorporate these concepts into our network framework.
Firstly, we should revisit the concept of interspike training, which pertains to the time intervals
between successive action potentials (also known as spikes) of a neuron within a period denoted
as T. An interspike train is represented by nin and encompasses the activities of other neurons
in the network after the time T has transpired. Given our division of Nt into g-octahedronsm[nmin ]; (1 ≤ m ≤ g), we designate In as the interspike train of neuron n ∈ [nin ]. Additionally,
we define the interspike train of the octahedron [nmin ] as the union of the interspike trains In ,Smmm
denoted as Inin = In | n ∈ [nin ]. Inin signifies the collection of interspike frequencies, or
spiking activity, within the octahedron [nmin ], representing a possible state of network activity
after the time T. Inmin constitutes a macrostate, representing a potential configuration of spiking
activities within Nt. Consequently, octahedrons with a larger number of neurons are expected to
exhibit more substantial spiking activity compared to those with fewer neurons. Thus, they are
more likely to approximate network spiking activity over a simulation period T. Once we havemlinked each octahedron [nmin ] to its respective macrostate [Inin ], we need to establish the conceptof microstates Ω. To simplify the concept and facilitate calculations, we can proceed as follows.
Since a microstate within [Inmin ] corresponds to the interspike train In generated by a single neuron n, the macrostate with the highest number of microstates is associated with the octahedron[nmin ] that has the most neurons. If an octahedron possesses a greater number of neurons, its
associated macrostate Inmin will have a larger number of interspike trains (microstates) In . Consequently, as the number of neurons in an octahedron increases, its entropy also increases. The
cardinality of an octahedron can be employed to represent the number of microstates. Finally,30

in accordance with equation (1.3.3), we can express this as follows:

S([Inm]) = ln(card([nmin ]))
in(1.3.4)

Where:• S([Inm]) = the entropy of Inmin (octahedron [nmin ]).
in• Ω = card([nmin ]) .
As in figure 1.19 panel B, for Nt = (11 × 11 × 11) with nin = 12, the macrostate In812 has thehighest entropy and we expect it to best approximate the interspike distribution (spiking activity)of the network Nt. Given the toroidal shape of the models, we can estimate the octahedron withthe highest number of neurons as that one formed by a mMd=g/2, where g is the maximum
mMd.
We refer to such distance as ζ: where:ζ = (g/2) + 1; if g is an even numberζ = (g + 1)/2; if g is an odd number(1.3.5)

We refer to the octahedron nζin as the class with the macrostate Inζ in with the highest entropy.

1.3.4

Shannon entropy

The mathematical representations of entropy as formulated by Boltzmann and Gibbs exhibit
similarities with Shannon’s concept of information entropy, which was developed during the1940s [51]. Gibbs entropy, which pertains to the probability distribution of microstates, serves
as a generalization of Boltzmann’s entropy [52]. It serves to quantify the extent of disorder
or randomness within a physical system. In a parallel vein, Shannon entropy, also referred
to as information entropy, provides a metric for the average information content or uncertainty
associated with a random variable. The computation of Shannon entropy is expressed as follows:
H(X) = −nX[pi (x) × log2 (pi (x))](1.3.6)i

The logarithm base 2 results in entropy measurements in bits [53]. It is maximum when all
possible outcomes are equally likely (a uniform probability distribution). It is zero when there
is only one possible outcome (a deterministic distribution). In the context of neuron activity,
Shannon entropy can be used to quantify the level of information or variability in neurons’ firing
patterns within a population. To calculate the Shannon entropy of neuron activity, we would
first discretize the ISI distribution of the neurons in each octahedron. Then, we can computethe probability distribution of spike occurrences during the simulation time T by counting thenumber of spikes in that interval and normalizing it by the total number of spikes observed. Once31

we obtain the probability distribution, we can calculate the Shannon entropy using the formula:H=−X(p([niin ]) × log2 (p([niin ])))

(1.3.7)

The pi represents the probability of spike occurrence in the octahedron [niin ]. A higher Shannon
entropy value indicates a more diverse or unpredictable firing pattern. This suggests a larger
amount of information contained in neuronal population activity. Conversely, a lower entropy
value implies a more regular or predictable firing pattern [54].

32Chapter 2
Chapter 2: Results
2.1

Strengthening of the SGE effect.

We now compare the evidence of the SGE effect on three different networks and confirm themathematical formulation expressed in eq.(1.2.5). To verify the SGE effect from one dimension
to another, the choice of nin and no must be done accordingly. As pictured in figure 1.15, in
panel (b) we chose no to be in a different plane than its counterpart in panel (a), where, even if
the mMd is 3 in both cases, the number of paths duplicates in the 3D grid model. We consider
three network grids (2D, 3D, and 4D) of regular firing (RS) neurons. First, we examine thesimulation’s spike timing data. Depending on the dimension of the model, the heterogeneous
spatial structure of the grid reduces the time of propagation of spiking information.

Figure 2.1: Comparison of regression slopes for three-dimensional toroidal networks in 3D (blue), 4D(red), and 5D (green). The simulations show the variation of the central delay (cd)j from 1 to 71 msec
in 2 msec steps. For each delay value, we executed 40 simulations. By changing the grid’s spatial
randomness, we measure the first spike’s arrival time. These 40 values are linearly interpolated, and thenegative value of the slope reveals the SGE effect; in other words, it shows that the interval between thefirst spike of nin and the first output of no decreases as the spatial heterogeneity of the grid increases. In
all three models the nin is subjected to the same current I=10 mA. All the neurons in the network havesynaptic strength w=18. We used the Izhikevich model implemented in pyNest. For each network grid
(2D, 3D and 4D), we chose (nin ; no ) to have the same Manhattan distance mMd. By maintaining thesame mMd and increasing the network dimension, the number of paths with the same mMd will increase
and strengthen the SGE effect.

We changed the intrinsic axonal delay (cd)j from 1 up to 71 milliseconds and tested twenty
different levels of noise delay (nd)i for each (cd)j . In figure 2.1 we compare the regression
33

T.D.
3D
4D
5D

I
10
10
10

Nt. size
11×11
7×7×7
5×5×5×5

w
18
18
18

nin
12
12
12

no
116
155
296

mMd
6
6
6

Paths
6
20
180

Slope
-1.09
-2.05
-3.42

Table 2.1: The characteristics of the multidimensional neural networks used in the simulations reported
in figure 2.1

OCTAHEDRONS
11
12
13
14
12
Octahedron [n112 ] [n212 ] [n312 ] [n412 ] [n512 ] [n612 ] [n712 ] [n812 ] [n912 ] [n1012 ] [n12 ] [n12 ] [n12 ] [n12 ] [n15 ]Cardinality 6
Spike
136

18
364

38
752

66
102 140 168 180 176 156 120801257 1811 2544 2988 3162 2994 2633 1964 1309

48
788

Table 2.2: Display of the data of the histograms in figure 2.2. For each octahedron, ’Cardinality’ represents the number of neurons that belong to the octahedron, while ’Spike’ represents the number of spikes
within the corresponding octahedron. The class [n112 ] has 6 elements that perform 136 spikes during simulation time T.

coefficient for the delay difference time ∆f in the case of three different networks: 3D, 4D, and
5D. Every single point of a line represents the slope resulting from 20 simulations with a fixed
central delay. The (nd)i varies from 0 to cd.
The negative value indicates a faster propagation of the signal. The three networks have thesame initiator neuron nin in position 12, driven by a constant 10 mA current. Across the wholenetwork, all neurons are connected with constant synaptic strength w = 18. However, the output
no has been chosen to have the same Manhattan distance mMd from nin in all three different
models. We notice that while the dimension increases, the number of paths increases according
to (1.2.5). Respectively for the networks of 2D (Torus 3D), 3D (Torus 4D), and 4D (Torus 5D),
no is set up to be at positions 116, 155, and 296. This results in 6, 20, and 180 paths with thesame mMd=6, for each model. Noticeably, the average slope value decreases as the dimension
rises, from -1.09 (3D Torus) to -2.05 (4D Torus) and -3.42 (5D Torus). Based on our discussionof equation 1.2.5, information propagates faster in higher-dimension networks. In other words,
increasing the dimension increases the number of paths from nin to no , and increases the numberof spikes arriving at no , resulting in a higher probability for the SGE effect to occur.

2.2

Network entropy

We consider the equivalence relationship mMd on Nt which subdivides the network into octahedrons, equivalent classes. The partition of Nt in its quotient set N t = N t/mM d, the set of
all octahedrons, permits us to compare the system entropy with the entropy of every octahedron
as the axonal heterogeneity (noise) increases. To do this we recall the functions of Boltzmann
and Shannon, as well as the equations eq.(1.3.3) and eq.(1.3.6), to calculate the entropy for eachoctahedron and of the network. We compare the entropies of their respective interspike trains
during their simulations and plot the results.
34

24
403

8
155

Interspikes time intervals

Figure 2.2: Illustration of the fifteen octahedrons of a network of size 11 × 11 × 11 (1331 neurons)of full neural heterogeneity. The initiator neuron nin =12 is supplied with a current of I= 10.0 mA andsynaptic strength w=18. The network is subdivided into octahedrons (equivalence classes), according to
the equivalence relationship mMd which ranges from 1 to 15. The number in the legend of each histogram
represents the corresponding octahedron class. For instance, in the top left panel the number 1 stands forthe octahedron [n112 ], and so on. The macrostate (octahedron in this case) with the highest entropy is
the class [n812 ] (panel number 8), with a number of microstates =180 and entropy S([n812 ]) = ln(180),
which we expect to best describe the network’s entropy. The class with the highest entropy represents
the macrostate towards which the network is most likely to move. As shown in the picture 1.19 and thedata in table 2.2, the behavior of the macrostates’ interspike frequency is proportional to their cardinality.
Class [n812 ] has the greatest cardinality (180 neurons) and the strongest spiking activity.

Boltzmann entropy.
The equivalent relationship mMd induced on Nt, a quotient set whose classes are named octahedron. According to the analytical expectation of Boltzmann’s equation for statistical entropy,
the macrostate with the highest entropy is the one with the highest value of microstates. The
macrostate of the system is defined as the set Inmin of the interspike trains of all of the neurons of
the octahedron [nmin ] whose cardinality is the number of neurons of the octahedron itself. Consequently, by extrapolating the cardinality of each octahedron we can predict the macrostate withthe highest entropy. This is, equivalently, the interspike-train distribution the system is most
likely to assume. Given the toroidal shape of the model, we can estimate the octahedron withthe maximum number of neurons as shaped by a mMd=ζ as in eq.(1.3.5).
Accordingly, we refer, respectfully, to [ninζ ] and Inζin as the octahedron with the highest numberof neurons and the macrostate with the highest entropy.
After running a number of simulations, we compare the results with analytical expectations. We
have constructed a neural heterogeneous network of size (11 × 11 × 11) with Izhikevich neurons
in figure 2.2. The initiator neuron nin =12 was connected to a constant I= 10.0 mA and synaptic
strength w=18. Specifically, figure 2.2 depicts the study of the interspike train, spiking activity,
for a central delay (cd) of 20 ms and noise delay (nd) of 16 ms. The maximum mMd (g) is35

Figure 2.3: The plot refers to a heterogeneous neural network. Izhikevich’s neurons are not identical.
From the network described in figure 2.2, the top panel (blue) depicts the interspike frequency of the8 ], the macrostate
network, and the bottom panel (green) depicts the interspike frequency for the class [n12with the highest entropy. The histograms refer to four different simulations of time T: 1000 ms, where
the heterogeneous axonal delay is calculated over a cd= 20 ms and four different values of the noise delay
δ, respectively 0, 4.7, 10.0, and 17.0 ms. As each legend reports the number of total spikes, we can seethat by changing δ, the interspike frequency distribution of the octahedron [n812 ] tracks that one of the
networks. As δ increases, the interspike frequency distribution follows that of the network. This indicatesthat noise delay has a large impact on the overall activity of the octahedron.

Figure 2.4: Similarly, in figure 2.3, we compare the frequency distribution of the class [n812 ], bottom
panel (green), and of the network, top panel (blue). Neural homogeneity case. All neurons spike regularly(RS). Also in this situation, we observe that by changing δ, the interspike frequency distribution of theoctahedron [n812 ] tracks the distribution of that network.

36

Figure 2.5: Heterogeneous neural network. The plot shows the distribution of the average entropy ofi ] over 20 simulations for a cd=20 ms as the noise ranges from 0 to 20 ms.
the equivalent classes [n12We determine the probabilities associated with each class, pi = p([ni12 ]), and calculate their entropy,
according to H([ni12 ]) = −pi × log2(pi ). We take the average entropy for each class across the 20
simulations centered on a central delay of cd=20 ms. Each value is a term of the average entropy of theP
network of equation H(X) = − ni [pi × log2 (pi ] whose values are plotted in figure 2.6. We can observe
how the highest value belongs to the class [n812 ], which has the highest number of neurons. Having thehighest entropy means that the class exhibits the highest level of unpredictability/information content. In
the legend, the value 3.51 is the average entropy of the network for a cd=20 ms noise range. It corresponds
to the sum of the average entropies of the classes (plotted values).

Entropy value of heterogeneous Nt:11 x 11 x 11 for cd: 20 ms.
3.54

3.53

Average entropy

3.52

3.51

3.50

3.49

3.48

3.47

0.0

2.5

5.0

7.5

10.0

12.5

15.0

17.5

Central delay: 20 ms

Figure 2.6: Heterogeneous neural network. As in figure 2.5, we consider the average entropy over 20
simulations for a cd=20 ms as the noise ranges from 0 to 20 ms. In this case, each value is the averagevalue of the network‘s entropy. The legend reports a final average of 3.51, as we expect from figure 2.5.

37

Figure 2.7: Heterogeneous neural network. The relationship between the network’s average number of
spikes and its average Shannon entropy is reported as the central delay cd ranges from 10 ms to 90 ms.
For each cd, as in figure 2.6, we calculated the average value, both for the number of spikes and entropy,
over 20 simulations as the noise ranges accordingly.

calculated as [(ln/2) − 0.5] × 3. Since the side ln=11 and the network dimension is three, itfollows g=15.

The quotient set N t, in this case, is subdivided into 15 octahedrons and ζ=8. Figure 2.2 depicts
the interspike-train distribution during a simulation time T, for each of the fifteen octahedrons.
Table 2.2 reports the number of neurons and the total number of spikes in every octahedron
after a time T. We can observe how the interspike-train distribution varies by class and shows
the maximum number of spikes in the octahedron [n812 ]. The macrostate I12ζ of the octahedron[n812 ], is the one with the highest entropy. It shows a number of microstates =180 and entropy
S(I12ζ ) = ln(180).
We aim to conduct a comparative analysis of the interspike-train distribution of I12ζ with that of
the network Nt. Figure 2.3 illustrates the interspike-train distributions for Nt (depicted in blue atthe top) and the macrostate I12ζ (depicted in green at the bottom). This analysis is conducted under conditions of neural heterogeneity and involves four distinct values of axonal delay (noise)with a fixed value of cd equal to 20 ms. When altering the noise levels, both the interspike-train
distributions for Nt and I12ζ exhibit proportional changes over time T. In Figure 2.4, we present
the interspike-train distributions for the same axonal delay parameters, but within a homogeneous neural network, where all neurons follow the Izhikevich regular spiking (RS) pattern. In
this scenario, the proportion of interspike-train distributions for both Nt and I12ζ remains relatively constant, irrespective of whether the number of spikes decreases significantly. Thus,
even in a homogeneous system, the proportion of interspike-train distributions remains consistent, regardless of the level of noise. This observation suggests that while noise can enhance the38spiking activity of the system, it does not significantly impact the relative proportion ρ of theinterspike-train distributions. This underscores the robustness of the subdivision of the system
into octahedrons, maintaining stability even in the presence of noise. In both cases denoted as(H) and (RS), as anticipated based on the equation (1.3.3), the octahedron structure closely approximates the interspike frequency distribution of Nt. An analysis of the interspike frequency
distribution reveals that in both cases, the octahedron structure provides the closest match to
the Boltzmann distribution, which is the expected distribution for a system in equilibrium. Thisfinding confirms that the octahedron structure represents the most energetically favorable configuration under conditions of high entropy.
Shannon entropy
Shannon entropy serves as a metric for gauging the level of uncertainty or informational content contained within a dataset. Information theory enables the quantification of the average
information necessary to describe or forecast an event. In the context of neural activity and theconcept of equivalence classes, Shannon entropy offers valuable insights into the variability anddiversity of neural responses. The resultant Shannon entropy value symbolizes the average levelof information or uncertainty linked to the equivalence classes of neural activity. Higher entropy
values signify an elevated degree of diversity and variation among these classes, whereas lower
entropy values imply more predictable or uniform patterns.
In our model, we employ the computation of Shannon entropy to evaluate the interspike train
activity of equivalence classes and compare it with the entropy associated with the network. As
for Boltzmann entropy, we consider the equivalence classes determined by the ”mMd” criterion. The interspike train of each neuron within an equivalence class is analyzed to estimate its
interspike activity, reflecting a temporal simulation of the class. Subsequently, the probability
distribution of the equivalence classes is calculated by normalizing the data points attributed
to each class with the interspike activity of the entire network denoted as Nt. Using the probabilities derived in the preceding step, we calculate the Shannon entropy of the network, as
expressed in equation (1.3.7). Here, p([niin ]) represents the probability of each equivalenceclass, and the summation encompasses all equivalence classes. Each component of this sum,−(p([niin ]) × log2 (p([niin ]))), is regarded as the entropy of the corresponding equivalence class.
Figures 2.5, 2.6, and 2.7 illustrate the Shannon entropy outcomes associated with equivalenceclasses and the activity of the neural network. Figure 2.5 illustrates the entropy distributionof the octahedrons and aligns with the Boltzmann entropy of equivalence classes. Notably, theequivalence class [n812 ], featuring the highest Boltzmann entropy value, also exhibits the highest
Shannon entropy value, indicating a high degree of diversity and variability among the classes.
Conversely, lower entropy values suggest more predictable or uniform patterns. Figure 2.6
presents the average network entropy value for a central delay of 20 ms across 20 simulations,with noise varying from 0 to 20 ms. In Figure 2.7, we examine the relationship between the
average Shannon entropy of the network and the average interspike network activity. In thissimulation, the central delay varies from 10 ms to 90 ms, and for each ”cd,” the noise varies from39

Figure 2.8: Brain neural heterogeneity, H. Keeping the parameters of the network in figure 2.2, the plot
depicts the values of ρ for each octahedron [ni12 ] (1 ≤ i ≤ 15) by varying the noise delay δ in 20 steps
from 0 to 20 ms for a cd= 20 ms.
The value of ρ represents the ratio between the total number of spikes of the network and the total numberof spikes of the class [ni12 ]. We can notice how the ratio remains quite constant for those octahedrons
close to the class with the highest entropy [n812 ].

0 to 20 ms. When dealing with systems characterized by spikes, each spike represents a specificoutcome or state. A reduction in the number of spikes indicates a decrease in system diversity
and the likelihood of certain outcomes. Consequently, the system becomes more predictable,
and entropy decreases proportionally. High entropy values signify increased variability, while
low entropy values indicate more predictable or uniform patterns. Consequently, the selectionof equivalence classes and grouping criteria must be made judiciously to ensure the capture ofdesired data patterns.
Conversion factor ρ.
Furthermore, we report the ratio ρ between the number of spikes of Nt and its macrostate, as in
2.2. In figure 2.8 (neural heterogeneity), we plot the value of ρ for each octahedron of figure
2.2 as the value of noise increases from 0 to 20 ms, with a cd= 20 ms. We can notice how
the closer we get to the octahedron at the highest entropy, [n812 ], the value of ρ tends to beconstant while fluctuating for those far from it. The macrostates at higher entropy, therefore,
offer a reasonable approximation of network spiking activity when axonal heterogeneity (noise)and neural heterogeneity are taken into account. This is because macrostates at high entropy
are more stable, meaning fluctuations in neuronal activity due to heterogeneity and noise areminimized. The macrostates at higher entropy more accurately reflect the network’s spiking
activity.

40

General Discussion and Conclusion
We constructed three different model networks framed on a toroidal grid in multiple dimensions. This realization helps to mimic realistic brain patches where formal limits and edges are
not defined, constraining data traffic in borderless areas. In the first part of the manuscript, we
generalized what was found in previous studies related to the influence of heterogeneous axonal
delay on various dimensions. This was to facilitate information transmission within a modelneural network. We connect the initiator neuron nin to an external current I, forcing it to regularly spike during the simulation time T and, consequently, to deliver signals to the adjacentnodes on the grid. The signal diffusion varies according to the grid structure, based on a von
Neumann neighborhood model in 2D, 3D, and 4D dimensions. In all the simulations undertakenin this study, for all the three different dimensions of the network, the negative value of the regression line’s slope m, figure 2.1, indicates how heterogeneous axonal delay (noise) induces
faster propagation speed within the network. The robustness by which we find negative m as
the studied models vary regarding dimensions and size, confirms the solidity of what we call thestochastic grid enhancement SGE effect. We extended the properties of this counter-intuitive result, generalizing the concept for any network dimension, showing that information contained in
spikes travels faster under higher heterogeneous axonal delay randomness. Through an analytical demonstration and several simulations, we have shown that at any dimension the increaseof heterogeneity, represented as a noise level in the axonal propagation delay, helps a signal(information) to spread within the network. Moreover, we investigated the role of neuron type
heterogeneity in the same set of spiking neural networks. We did this by manipulating the parameters c, d of the Izhikevich model through a uniformly distributed random variable. The
level of neural heterogeneity ranges from 0 to 1, with 0 representing a perfectly homogeneous
regular spike (RS) neuron model network and 1 a heterogeneous one as defined in equations
1.2.1. In all cases, we set synaptic weights at w = 18, and the initiator neuron nin input current
at 10 mA.
Recognizably as shown in figures 2.9, by keeping constant the initial value of current, synaptic
strength, and the number of neurons of the model, introducing cell type heterogeneity results ina large increase in firing activity a phenomenon that we called here stochastic neuron enhancement, SNE effect. The information created from the spiking activity was abundant in heterogeneous networks since the networks could increase spiking activity in the same period and withequal beginning values for the number of neurons, current, and synaptic strength. As a result, weconclude that the combination of axonal and neuronal heterogeneity is a metabolically efficient
brain strategy. Heterogeneous networks have no additional cost in terms of neurons or synapses
and outperform homogeneous networks with the same order of initial parameters. This neuralcell heterogeneity is therefore a biologically and computationally effective strategy for the brain
to improve signal robustness.
In addition, we introduced the concept of equivalent classes in a neural network, proving theminimum Manhattan distance (mMd) to be an equivalence relation on the network. Neurons are41

considered equivalent if they have the same mMd from nin neuron. These equivalent classes arereferred to as octahedrons, and interspike trains are defined for each octahedron. The Boltzmann
and Shannon entropy of the network are analyzed based on these equivalence classes, allowing
for a higher-level understanding of the network’s behavior and the identification of common
optimization techniques. The theoretical model partitions the network into classes associatedwith their entropy, providing an alternative approach to evaluating network entropy and spike
information. This partitioning helps calculate the entropy of each state and identifies the mostprobable state, optimizing computational performance by avoiding unnecessary calculations.
These findings contribute to understanding the fundamental characteristics of information transfer in the brain, providing a conceptual model for an explanation of other paradoxical phenomena
in which neural heterogeneity and heterogeneous axonal delay appear to favor computational
processes enhancing perception and other performance in human subjects. While in this study
we kept parameter manipulation to a minimum to better clarify the effect of heterogeneity on
information transmission. We also suggested further statistical analyses such as synchronization
factors, discharge frequency, and other possible statistical analyses to deepen the understandingof spike information transmission.

42

Figure 2.9:
A: A comparison of the SGE effect over two models of neural heterogeneity (H=1, blue) and a regularspiking model (H=0, R.S., red). According to panel B, the H=1 line shows the greatest amount of spiking
activity, while the R.S. is the minimum. Nevertheless, the SGE effect acts similarly on both networks.
All the models have the same constant current input of 10 mA, a synaptic strength of 18, and the same
number of neurons. Noticeably, the heterogeneity offers the most suitable model regarding robustness(highest spiking activity) while keeping the benefit of improved spike propagation times (SGE effect).
B: The panel assesses the variations in spike counts at various central delays from 1 to 71 ms and acrosssix different levels of heterogeneity. At each central delay value cd, the degree of randomness is adjusted
in 15 incremental steps, as depicted in Figure 1.17. The total spike count is computed across all neurons
in the network during the simulation period. Notably, when the neuronal heterogeneity is at its maximum,
the model exhibits the most resilient spiking activity, which we refer to as the SNE effect43Compliance with Ethical Standards
Conflict of Interest: Marcello Salustri declares that he has no conflict of interest. Ruggero
Micheletto declares that he has no conflict of interest.
Ethical approval: This article does not contain any studies with human participants performedby any of the authors.

44

List of Publications
Referred journals:
1. Salustri, Marcello and Micheletto, Ruggero, “Heterogeneous Axonal Delay Improves the
Spiking Activity Propagation on a Toroidal Network”, Cognitive Computation (2022),
Springer. https://doi.org/10.1007/s12559-022-10034-2.
2. Salustri Marcello and Yoshida Shunra and Micheletto Ruggero, “Neural and axonal heterogeneity improves information transmission”, Physica A: Statistical Mechanics and its
Applications (2023), Springer. https://doi.org/10.1016/j.physa.2023.
128627.

Under submission:
1. Salustri, Marcello and Micheletto, Ruggero, ‘Partitioning of entropy in equivalence classes.’,Cognitive Computation (2023), Springer. .
Submitted on July 31, 2023.

Papers communicated in referred journals:
1. Sun, Zhe and Cutsuridis, Vassilis and Caiafa, Cesar F. and Solé-Casals, Jordi, “Brain Simulation and Spiking Neural Networks”, Cognitive Computation (2023), Springer. https://doi.org/10.1007/s12559-023-10156-1.

45

Bibliography[1] J.F. Pagel and P. Kirshtein. Chapter six - neural networks: The hard and software logic.
pages 83–92, 2017.
[2] A. Tavanaei, M. Ghodrati, S.R. Kheradpisheh, T. Masquelier, and M. Maida. Deep learningin spiking neural networks. Neural Networks, 111:47–63, 2019.
[3] S. Feng and S. Zheng. How can artificial neural networks approximate the brain? Frontiersin Psychology, 13, 2023.
[4] T. Hongwei, Z. Yifan, Quanzheng, R. Johanna, and S. T. van Dijken. Bioinspired multisensory neural network with crossmodal integration and recognition. Nature Communications,12, 2021.
[5] E. Manuylovich, D. Argüello Ron, M. Kamalian-Kopae, and S. Turitsyn. Stochastic resonance neurons in artificial neural networks. 2022.
[6] M. Kawaguchi, H. Mino, and D. Durand. Stochastic resonance can enhance informationtransmission in neural networks. IEEE transactions on bio-medical engineering, 58:1950–
8, 03 2011.
[7] S. Radhakrishnan, S.and Sebastian, A. Oberoi, D. Sarbashis, and D. Saptarshi. A
biomimetic neural encoder for spiking neural network. Nature Communication, 12, 4 2021.
[8] M. D. McDonnell and D. Abbott. What is stochastic resonance? definitions, misconceptions, debates, and its relevance to biology. PLoS computational biology, 5, 4 2009.
[9] O. van der Groen, M. F. Tang, N. Wenderoth, and J. B. Mattingley. Stochastic resonanceenhances the rate of evidence accumulation during combined brain stimulation and perceptual decision-making. PLOS Computational Biology, 14(7):1–17, 07 2018.
[10] Harikrishnan N.B. and Nithin Nagaraj. When noise meets chaos: Stochastic resonance inneurochaos learning. Neural Networks, 143:425–435, 2021.
[11] Huixia Liu, Lulu Lu, Yuan Zhu, Zhouchao Wei, and Ming Yi. Stochastic resonance:
The response to envelope modulation signal for neural networks with different topologies.
Physica A: Statistical Mechanics and its Applications, 607:128177, 2022.
[12] W.C. Stacey and D.M. Durand. Stochastic resonance in simulated and in vitro hippocampal
ca1 cells. 1:364 vol.1–, 1999.
[13] F. A. Azevedo, L.R. Carvalho, L. T. Grinberg, J. M. Farfel, R. E. Ferretti, R. E. Leite,
W. Jacob Filho, R. Lent, and S. Herculano-Houzel. Equal numbers of neuronal and nonneuronal cells make the human brain an isometrically scaled-up primate brain. J. Comp.
Neurol., 10(513):532–541, April 2009.
46

[14] AT. Winfree. The geometry of biological time. New York: Springer, 2001.
[15] C. Koch and G. Laurent. Complexity and the nervous system. Science, 284:96–98, 1999.
[16] S. G. Tewari, M.K. Gottipati, and V. Parpura. Mathematical modeling in neuroscience:
Neuronal activity and its modulation by astrocytes. Frontiers in Integrative Neuroscience,
10, 2016.
[17] A. Aldo Faisal, Luc P. J. Selen, and Daniel M. Wolpert. Noise in the nervous system. Nat
Rev Neurosci, 9, APR 2008.
[18] Z. Budrikis. Forty years of stochastic resonance. Nature Reviews Physics., 3(771), 2021.
[19] WC. Stacey and DM. Durand. Stochastic resonance improves signal detection in hippocampal ca1 neurons. J Neurophysiol, 83(3), MAR 2000.
[20] N. Perez-Nieves, V.C.H. Leung, and P.L. Dragotti. Neural heterogeneity promotes robust
learning. Nat Commun, 12(5791), October 2021.
[21] B. Moret and M. et al. Donato, R.and Nucci. Transcranial random noise stimulation (trns):
a wide range of frequencies is needed for increasing cortical excitability. Scientific Report,
9, 2019.
[22] J. von Neumann. Probabilistic logics and the synthesis of reliable organisms from unreliable components. Automata Studies., 34:43–99., 1956.
[23] Richard J. Gardner, E. Hermansen, M. Pachitariu, Y. Burak, N. Baas, Benjamin A. Dunn,
M. Moser, and Edvard I. Moser. Toroidal topology of population activity in grid cells.
Nature, 602, 2022.
[24] G. Tononi, A.R. McIntosh, D.P. Russell, and G.M. Edelman. Functional clustering: Identifying strongly interactive brain regions in neuroimaging data. Neuroimage., 7:133–149,
1998.
[25] T.O. Sharpee, A.J. Calhoun, and S.H. Chalasani. Information theory of adaptation in neurons, behavior, and mood. Curr. Opin. Neurobiol., 25(25):47–53, 2014.
[26] Z. Huang, J. Zhang, J. Wu, G.A. Mashour, and A.G. Hudetz. Temporal circuit ofmacroscale dynamic brain activity supports human consciousness. Sci. Adv., 6(eaat0087),
2020.
[27] A. Demertzi, E. Tagliazucchi, S. Dehaene, G. Deco, P. Barttfeld, F. Raimondo, C. Martial,
D. Fernández-Espejo, B. Rohaut, H.U. Voss, and et al. Human consciousness is supported
by dynamic complex patterns of brain signal coordination. Sci. Adv., 5(eaat7603), 2019.
[28] CJ. Schwiening. A brief historical perspective: Hodgkin and huxley. J Physiol., 590(11),2012.
[29] E.M. Izhikevich. Simple model of spiking neurons. IEEE Transactions on neural networks,
14(6), 2003.
47

[30] Chang Wu, Yubai Li, and Song Chai. Design and simulation of a torus topology for network
on chip. Journal of Systems Engineering and Electronics, 19(4):694–701, NOV 2008.
[31] J. Eppler, M. Helias, E. Muller, M. Diesmann, and M.-O. Gewaltig. Pynest: a convenientinterface to the nest simulator. Frontiers in Neuroinformatics, 2(3), 2009.
[32] MO. Gewaltig, A. Morrison, and Plesser H.E. Nest by example: An introduction to theneural simulation tool nest. Computational Systems Neurobiology, pages 533–558, 2012.
[33] Soni Chaturvedi. Comparison of lif and izhikevich spiking neural models for recognitionof uppercase and lowercase english characters. CiiT International Journal of Digital Image
Processing, 6, 07 2014.
[34] S. Nobukawa, H. Nishimura, and T. Yamanishi. The importance of neighborhood scheme
selection in agent-based tumor growth modeling. Scientific Reports, 8, 2018.
[35] L. Fang, N. Zhao, L. Caudal, H. Chang, R. Zhao, C. Lin, N. Hainz, C. Meier, B. Bettler,
W. Huang, A. Scheller, F. Kirchhoff, and X. Bai. Impaired bidirectional communication
between interneurons and oligodendrocyte precursor cells affects social cognitive behavior. Nature Communications, 13, 2022.
[36] M. Coli, P. Palazzari, and R. Rughi. The toroidal neural networks. 2000 IEEE International
Symposium on Circuits and Systems (ISCAS), 4:137–140 vol.4, 2000.
[37] R. Suwanda, Z. Syahputra, and E. M. Zamzami. Analysis of euclidean distance and manhattan distance in the k-means algorithm for variations number of centroid k. Journal of
Physics, 1566:696, NOV 2019.
[38] A. Schneider and M. Hommel, G.and Blettner. Linear regression analysis: part 14
of a series on evaluation of scientific publications. Deutsches Arzteblatt international,
107(44):776–782., 2010.
[39] SM. Chrysafides, S. Bordes, and S. Sharma. Physiology, resting potential. NCBI National
Center for Biotechnology, 2021.
[40] Stephanie B. Aldrich. The use of multiple neurotransmitters at synapses. Synaptic Transmission, pages 449–480, 2019.
[41] A. Tozzi and F Peters. James. Towards a fourth spatial dimension of brain activity. Cognitive neurodynamics, 10(3), 2016.
[42] B. Sonal and G. Mahendra. Design of mesh and torus topologies for network-on-chip
application. International Journal of Reconfigurable and Embedded Systems (IJRES), 2,
07 2013.
[43] D. A. Zaitsev. A generalized neighborhood for cellular automata. Theoretical Computer
Science, 666:21–35, 2017.
48

[44] M. Madadi Asl, A. Valizadeh, and P. Tass. Dendritic and axonal propagation delays determine emergent structures of neuronal networks with plastic synapses. Scientific Report, 7,
2017.
[45] M. Salustri and R. Micheletto. Heterogeneous axonal delay improves the spiking activity
propagation on a toroidal network. Cogn Comput, June 2022.
[46] M. Sato and T. Sado. Lattice paths restricted by two parallel hyperplanes. Bulletin of
informatics and cybernetics, 21(3/4):97–105, 1985.
[47] T. Britz, M. Mainetti, and L. Pezzoli. Some operations on the family of equivalence relations. Springer Milano, 1 2001.
[48] O. Ore. Theory of equivalence relations. Duke Mathematical Journal, 9(3):573 – 627,
1942.
[49] C. Cercignani. The boltzmann equation. Springer New York, 40(103):279–300, 1988.
[50] K. Uchiyama. A fluctuation problem associated with the boltzmann equation for a gas ofmolecules with a cutoff potential. Japanese journal of mathematics. New series, 9(1):27–
53, 1983.
[51] Yeliz Karaca and Majaz Moonis. Chapter 14 - shannon entropy-based complexity quantification of nonlinear stochastic process: diagnostic and predictive spatiotemporal uncertainty of multiple sclerosis subgroups. pages 231–245, 2022.
[52] P. Županović and D. Kuić. Relation between boltzmann and gibbs entropy and example with multinomial distribution. Journal of Physics Communications, 2(4):045002, apr
2018.
[53] Alejandro A. Torres-García, Omar Mendoza-Montoya, Marta Molinas, Javier M. Antelis,
Luis A. Moctezuma, and Tonatiuh Hernández-Del-Toro. Chapter 4 - pre-processing andfeature extraction. pages 59–91, 2022.
[54] H. Nakahama, M. Yamamoto, K. Aya, K.vand Shima, and H. Fujii. Markov dependency
based on shannon’s entropy and its application to neural spike trains. IEEE Transactions
on Systems, Man, and Cybernetics, SMC-13(5):692–701, 1983.

49

ACKNOWLEDGMENTS
It is a great pleasure for me to express my respect and a deep sense of gratitude to my supervisor
Professor Micheletto Ruggero, Department of Physics Nanoscience College, Kanazawa-Hakkei,
Yokohama City University, for his wisdom, vision, expertise, guidance, enthusiastic involvement and persistent encouragement during the planning and development of this research work.
I also gratefully acknowledge his painstaking efforts in thoroughly going through and improving
the manuscripts without which this work could not have been completed.
I am highly obliged to all members of the Faculty of Science of Yokohama City University, forproviding all the facilities, help, and encouragement for carrying out the research work.
I am obliged to my parents Salustri Mario and Fondi Ecle for their moral support, love, encouragement, and blessings to complete this task. I am especially thankful to my wife Sato Marina,
my daughters Arianna and Sofia, and my son Leo for her patience, love, and encouragement
during this journey.
I also would like to express my deep and sincere thanks to my friends and all other persons
whose names do not appear here, for helping me either directly or indirectly in all even and odd
times.
I am also thankful to the anonymous reviewers of my research publications. Their comments
and suggestions were very helpful in shaping my research work.

Salustri Marcello50