Physica A 618 (2023) 128627

Contents lists available at ScienceDirect

Physica A
journal homepage: www.elsevier.com/locate/physa

Neural and axonal heterogeneity improves informationtransmission✩
Salustri Marcello, Yoshida Shunra, Micheletto Ruggero∗

Graduate School of Nanobioscience, Cognitive Informatics Laboratory, Yokohama City University, 22-2 Seto,
Kanazawa-ku, 236-0027, Japan

articleinfo

Article history:
Received 2 November 2022
Received in revised form 2 February 2023
Available online 2 March 2023
Keywords:
Heterogeneity
Neural
Noise
Delay
Spike

a b s t r a c t
The complexity of the brain lies in an intricate, not homogeneous structure; however,
it is not clear how important is the heterogeneity at neural and structural levels andhow it may play an important constituent in the brain’s functionality. Here we show
the role of cellular and axonal delay heterogeneity in the brain’s performance. We
simulated, at different scales, the spiking activity on a toroidal network realized in
multiple dimensions with varying degrees of heterogeneity on a network of Izhikevich
neuron models. We found that increasing the heterogeneity and network dimension
improved the robustness and propagation speed of the spiking activity. Our results
demonstrate that the behavior of the spiking activity depends on both the cellularneural and axonal delay heterogeneity. We developed a simple theoretical frameworkcompatible with the results of the simulations, putting forward a novel method tostrategically analyze any similar networks.
© 2023 Elsevier B.V. All rights reserved.

1. Introduction
It is still unclear whether or not the heterogeneity of the brain can aid the process of the brain‘s performance [1], asmost previous studies often employed models based on regular structures. Given the heterogeneity of the brain, [2], we
investigated the nature and impact of two types of heterogeneity with the realistic neural model using pyNest simulator ona toroidal network; the first heterogeneity refers to the axonal delay variability, a proxy for inter-neuron distance, and thesecond one to the variability of Izhikevich neuron type. The results of this study suggest that the observed heterogeneity
in the brain may be an important factor in the capacity to optimize its overall performance. Actually, we proved thatheterogeneity at all levels enhances performance in terms of signal propagation speed and increases the spike activity inthe net, in other terms the robustness of information circulating in the brain.
In the last decades, mathematical models to study the propositions that govern cognitive abilities of the nervous
system, helped significantly to unveil unknown mechanisms of critical behaviors in brain activities and related neuralnetwork models [3]. Additionally, the support of computational simulations describing biologically plausible neurons
and their dynamics, validated mathematical models of artificial neural networks that aim to capture neurons’ biologicalarchitecture and behavior in neuroscience and in many relevant fields [4]. In modern neuroscience, considerable
development of mathematical models facilitated the realization of dynamic systems resembling the biological functioningof the brain [5]. Among such studies, stochastic resonance, a counter-intuitive mechanism by which an embedded noise [6]✩ Supported by the Yokohama City University.

∗ Corresponding author.

E-mail address: ruggero@yokohama-cu.ac.jp (M. Ruggero).
https://doi.org/10.1016/j.physa.2023.128627
0378-4371/© 2023 Elsevier B.V. All rights reserved.

S. Marcello, Y. Shunra and M. Ruggero

Physica A 618 (2023) 128627

at a finite level enhances the sensibility of a system empowering its performance, is considered to play an active rolein a variety of classes of both natural and artificially structured neural systems [7]. Furthermore, given the potential
significance of this phenomenon, and in light of the previous study of how random noise stimulation (heterogeneous
axonal delay) improves cognitive performance [8], this article will provide a theoretical formulation of two distinct
heterogeneity effects. The first type of heterogeneity introduced is variable distances between neurons obtained byrandomizing axonal delay. This is leading to faster information propagation speed, a phenomenon that we call stochastic
grid enhancement (SGE). The second heterogeneity tested is concerning a non-homogeneous neural model. In this case, we
variate the neuron cell type, manipulating the Izhikevich model parameters that change the neuron’s various differential
properties resulting in spike patterns and dynamics.
We do not have direct evidence on the real effect of heterogeneity on a biological brain, however, our simulation
using the biological realistic spiking Izhikevich model has shown that heterogeneity enhances the number of spikes along
the whole network, and this suggests that information robustness is enhanced. Other recent prominent studies are also
compatible with our finding, as it seems that heterogeneity gives the brain various advantages, including efficient coding,reliability, better working memory, and other functional specialization [1]. Hereafter we will call this effect caused by the
randomness of the neuronal model stochastic neuron enhancement (SNE).
In this study firstly, we consolidate the SGE effect, already discussed extensively in a previous work [8], testing it
on multidimensional toroidal grids, secondly, we provide novel SNE tests and results and we discuss the mathematicalconnection between the functional role of neural type heterogeneity and the robustness of spiking activity within a
network [9]. To implement the study of heterogeneity we choose the Izhikevich neuron model. Despite many studies
advocating the use of the LIF (leaky integrate-fire) model, the Izhikevich model is considered one of the most accurate,
but at the same time is not computationally heavy and allows to simulate of a large number of neurons [10]. The model,
throughout two dimensionless variables u and v, and four dimensionless parameters a, b, c, and d, naturally fits therequirements to pursue the study of neural heterogeneity. In this article, therefore, we want to bring light to the roleof these two combined heterogeneities in the improvement of spiking activity and robustness of information propagatingin the neural network [11].
2. Methods
We examine the toroidal grid in two, three, and four dimensions. The construction of a 2D torus grid is realized byconnecting each neuron to its four von Neumann neighbors and connecting each edge of the grid to its opposite edge
in a cyclic manner realizing a torus that exists in a three dimensional space. This is a structure that does not alter the
topological properties of the original grid.
The 4D torus comes from the 3D grid where nodes (neurons) are connected in a similar von Neumann threedimensional structure, with six adjacent neighbors instead of four. The 4D torus is realized by merging the opposite
surfaces of the cube, obtaining again a toroidal network (this ‘‘torus’’ cannot be visualized in three dimensions; it isactually a Clifford structure in four dimensions [12]. With a similar scheme, we will also realize a four-dimensional grid
resulting in a torus placed in five-dimensional space. In the following scheme we summarize the relationship between
the dimensional grid and the corresponding toroidal network.

• 2D Lattice grid H⇒ Torus in 3D space• 3D Lattice grid H⇒ Torus in 4D space• 4D Lattice grid H⇒ Torus in 5D space
Biologically speaking, neurons in the brain live in confined areas without specified borders [13]. Accordingly, we choose
a borderless structure as a torus to confine the network in a model which may portray the most plausible representationof the brain. Particularly, we use the neurons model presented by Eugene M. Izhikevich [10], which accurately resembles
the spiking behavior of cortical neurons, reproducing a lot of variety of neuron spiking behaviors.
Our three-dimensional topological structure [14], is defined as a graph G = (N , C ), where N(G) and C (G) are respectively
the nodes and connections of G [15], each node is connected to its six adjacent neighbors and the communication of eachnode can take place in six directions: +x (east), -x (west), +y (forward), -y (back), +z (north), and -z (south).
The four-dimensional torus structure leads naturally to a von Neumann-shaped neighborhood [16] where cells are
surrounded by six inputs [17].
We chose arbitrarily an initiator neuron, which is regularly spiking due to its connection to a constant external currentof I = 10 mA that induces the network to spike for the duration of the simulation. All the other neurons in the grid are
not connected to external currents and may spike only because stimulated by spikes coming from neighbor neurons. Eachcouple of nodes represents a different distance between neurons, so it characterizes a delay or the transmission time fora spike signal to reach its neighbor. We refer to it as µ = cd + δ (cd: ‘‘central delay’’) and δ = nd × α (nd:noise delay),
where α = [(2x) − 1] ∼ U(−1, 1), and x∼ U(0, 1) being uniformly distributed random variables used to implement
the axonal propagation delay heterogeneity, the trigger of the SGE effect. Once proved this phenomenon on 2D, 3D, and4D multidimensional grids, it will be introduced and analyzed the effect of neural heterogeneity on the overall network
spiking activity testing the SNE effect on the grid. Finally, we will consider the bonding effect of the two heterogeneities
on the network and elucidate the results and conclusions of the study.
2

S. Marcello, Y. Shunra and M. Ruggero

Physica A 618 (2023) 128627

Table 1
The table reports the actual values of the parameters a, b, c, d used for each model of Fig. 2.
Heterogeneity of Izhikevich neuron modelParameters

RS

IB

CH

FS

LTS

TC

RZ

a :
b :
c :
d :

0.02
0.2
−65
8

0.02
0.2
−55
4

0.02
0.2
−50
2

0.1
0.2
−65
2

0.02
0.25
−65
2

0.02
0.25
−65
0.05

0.1
0.26
−65
2

Table 2The table summarizes the data of Fig. 7. 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 the interpolation line (Slope) for all the heterogeneity
values (from H:0.0 to H:1.0).
Nt.size

H.:0.0H.:0.2

H.:0.4

H.:0.6

H.:0.8

H.: 1.0

5×5×5×5

S.D.
Mean
Slope

55.3
3091
−201

55.3
3101
−190

55.7
3340
−228

56.6
3636
−231

57.7
4543
−345

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

2.1. Simulation
2.1.1. Neuron model
We used Python and the pyNEST library version 2.2 [18], to simulate spiking neurons network in continuous time
whose connections are adjustable parameters for weight and delay. 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 neuralfunctionality in the field of neuron science is the LIF (leaky integrate-fire) model, considered simple while implemented at
low computational cost [19]. 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 [20].
The model offers highly plausible biological dynamics, which realistically reproduce the non-linear phenomena dynamics
caused by the increased noise level in the network, and it has an outstanding capacity to reproduce a large variety ofspiking models by manipulating a few parameters. Therefore, since it naturally fits the requirements to pursue the studyof neural heterogeneity, it becomes our first choice. Izhikevich’s model, uses two dimensionless variables u and v, andfour parameters a, b, c, and d, is a two-dimensional system of ordinary differential equations which describes the time
evolution of the membrane potential v and membrane recovery u [10], described in the following equation:{ ′v = 0.04v 2 + 5v + 140 − u + I
u′ = a(bv − u)(1)with after the spike, a reset process is given by:
if v ≥ 30 mV →

{
v←c
u←u+d

(2)

where v ′ = ddtv , and u′ = du
. The variables v and u and the parameters a, b, c and d are all dimensionless, t is the time indtmilliseconds [10].
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 takes the presynaptic neuron’s input into
consideration, where the weights represent the synaptic connection strength parameters. When a neuron is delayed bya specific amount of time, its signal will continue to travel down a memory pipe register long as the delay is expressed
as a number of time steps.
The parameters of the Izhikevich model describe respectfully (see Tables 1–3):• a: time scale of the recovery variable u.
• b: the sensitivity of the recovery variable u to the subthreshold fluctuations of the membrane potential v.
• c: the after-spike reset value of the membrane potential v.
• d: The after-spike reset of the recovery variable u.
3

S. Marcello, Y. Shunra and M. Ruggero

Physica A 618 (2023) 128627

i,jFig. 1. Four different representations of 3D Torus from a 2D grid model. The nodes are the neurons of the network; the connections are the µn,nk
of Eq. (4). The x-axis scales the variation of neural heterogeneity H, from 0 to 1. H:0 every neuron is identical to each other (neural homogeneity).
On the y-axis ranges the axonal heterogeneity, the axonal delay δi , ranging from 0 to the value of the central delay (cd). To the variation of δii,jcorresponds a different value for the connection µn,nk .
A: Particular situation for H: 0, and δi : 0. The network grid is regular (axonal delay is homogeneous and neurons are of the same type). It resultsin having all neurons of the same type (in this case all red nodes indicate all the neurons are identical, specifically of R.S. Izhikevich type) and all
the axonal delay between them are homogeneous to the value of cd (regular 3D tours).
B: H: 1 and δi : 0. We change only the neural heterogeneity by manipulating the parameters c, d of the Izhikevich neuron model, resulting, according
to Fig. 2, on a stochastic typology of neurons model. The distance between neurons, hence, remains constant to the value of the cd, while thedifferent colors of each neuron correspond to a different type of Izhikevich neuron model.
C: H: 0, and δi : 1. All neurons are R.S. Izhikevich neurons, while the axonal heterogeneity provokes the axonal delay to vary between neurons, as
formalized by the Eq. (4).
D: H: 1, and δi : 1. The sketch depicts the most variable structure possible: the two heterogeneity, portraying the most plausible distribution ofneuron populations where the distance of the connection between neurons vary within a network of heterogeneous cells.

In our model we will keep a, b at constant values of 0.02 and 0.2, while the parameters c, d are heterogeneous, as inFig. 2, will be subject to the form:{c = −65 + (15 × x1 × x1 )(3)

d = 8 − (6 × x2 × x2 )0 ≤ x1 , x2 ≤ 1: uniformly distributed random variables [10].
In Fig. 1 we show the resulting effects of the combined heterogeneity on a 3D toroidal model network. Where thenode’s color variation represents the neural heterogeneity, the different distances of the connections illustrate the axonal
delay heterogeneity.
2.1.2. Network structure
A natural choice of the network framework, to secure the SGE effect on a 3D lattice grid, falls in the three-dimensionallattice von Neumann neighborhood (6-neighborhood model as in Fig. 3), where each cell is connected with its six adjacentcells located on its north, south, east, west, in forward and back position, where the distance between nodes is calculatedwith the Manhattan distance [16].
4S. Marcello, Y. Shunra and M. Ruggero

Physica A 618 (2023) 128627

Table 3The characteristics of the multidimensional neural networks used in the simulations reported in Fig. 9.
T.D.

Nt. size

I

w

nin

no

mMd

Paths

Slope

3D
4D
5D

11 × 11
7 × 7 × 7
5 × 5 × 5 × 5

10
10
10

18
18
18

12
12
12

116
155
296

6
6
6

6
20
180

−1.09
−2.05
−3.42

Fig. 2. To different values of the parameters, a, b, c, d, correspond different types of Izhikevich neuron models as described by Eqs. (1) and (2). 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, dcorrespond to a continuous variability of a large range of neuron models [10] 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. (3), guarantees the neural heterogeneity ofthe model.

Fig. 3. 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 areconnected with their respective counterparts on the opposite face. For clarity the brown dotted line represents one of those cyclic connections. This
structure is equivalent to a torus in 4D space.

5

S. Marcello, Y. Shunra and M. Ruggero

Physica A 618 (2023) 128627

We created 3 different networks of two, three, and four dimensions associated, respectfully, with a three, four, and fivedimensional torus. Hence, we ran comparative simulations to analyze the results. Initially, we consider a three-dimension
matrix of 343 neurons (7 × 7 × 7) connected according to the von Neumann neighborhood model such that the surfaces
are connected with their opposite to form a toroidal structure [21]. Once we have framed the network, we configured
the neurons according to Izhikevich’s archetype to start the simulations. As we aim to focus on the effect of non-uniform
axonal propagation delay [22], we first study the model in its simplest form considering all RS (regular spiking) neurons
to be identical, where no inhibitory synapses nor thalamic currents were examined. Hereafter we will compare the resultswith the same-size heterogeneous model and analyze the differences in the effect of heterogeneity.
As established in the previous study [8], the heterogeneous axonal delay (representing the variable time it takes thei ,jspike to travel the length of the axon) µn,nk between two neurons n, nk is calculated at each simulation by its averagevalue, that we call here central delay (cdj ) and its stochastic contribution noise delay (nd)i as follows:δi  
µni,,nj k = (cd)j + (nd)i αni,,nj k(4)where i and j are the index of the loops over (cd) and (nd) respectively. In other words, j represents a particular central
delay (that is the average delay between nodes) and i represents the cd randomness level added to each of the neuronsof the grid for that simulation.

∗ αni,,nj k = [(2xi,j ) − 1]
∗ (cd)j = central delay; 1 ≤ (cd)j ≤ (cd)max , j ∈ N∗ δi = (nd)i · αni,,nj k ; (δi ≤ (cd)j ) is the noise delay(cd)∗ for each j: (nd)i = i ns j , 1 ≤ i ≤ ns, i ∈ N where 0 ≤ xi,j ≤ 1 is an aleatory variable with uniform distribution
i ,jassigned to the connection of neuron n to nk at the i-loop and j-loop. In this way αn,nk is always within the −1, +1interval.
i ,jFirstly, we proceed to study the 3D lattice grid representing the network’s different (heterogeneous) axonal delay µn,nk .
Only one neuron called nin is stimulated by an external current I = 10 mA, as suggested by Izhikevich [10], whose spiking
activity stimulates its neighbors inducing secondary spiking to propagate along the model. We compute the propagation
speed between nin and a neuron called no representing an arbitrary output:i,j• foutas the time of the first spike of no .
i,j• fin as the time of the first spike of nin .

Hence, the spike propagation delay from nin to n0 is:i,j∆fg i,j = fo i,j − fin(5)with g being the minimum Manhattan distance (from now on referred to as mMd) from nin to no , between the twoi,jneurons, [23]. All the ∆fg of one i-loop, have the same central delay (cd)j = j which is increased at the end of one
i ,ji-loop. For the same (cd)j = j, by increasing the noise (heterogeneity of axonal delay), we previously proved that ∆fg
decreases. We refer to this phenomenon as the stochastic grid enhancement effect (SGE) and we expect its presence inmultidimensional models.
2.2. Grid dimensions and SGE effect
Here we want to generalize the stochastic grid enhancement on a larger scale, showing the phenomenon to appear
more prominent by increasing the number of dimensions of the network. Firstly, we discuss the number of minimum
paths with the same mMd between two random neurons in a lattice grid framed in a toroidal shape, where the mMdfrom nin to no , is the number of minimum steps to reach no from nin as pictured in Fig. 3.
In a previous 2D lattice network study [8], for each neuron there were four adjacent neighbors, and for any couple ofneurons [nin ; no ] (initiator; output), there were a certain number of paths with the same mMd connecting the nin to no .
Given the structure of the model, the fastest way for the signal to reach the output no is in mMd steps. In two dimensions
there are only two possible pre-synaptic connections to no (as in Fig. 4(a)), since starting a path from one of the othertwo connections, will take a number steps greater than mMd.
The strength of the SGE effect depends on the number of paths with the same mMd, which increases by manipulating
the number of dimensions of the network and the nodes’ coordinates as shown in the following.
Let us consider a 3D cubic lattice (4D torus) of side ln=2 ·sl, and (nin , no ) of coordinates [(x0 , y0 , z0 ); (x1 , y1 , z1 )].
We can formulate :{| x0 − x1 |; if : |x0 − x1 | ≤ sl
X =ln − |x0 − x1 | ; if : |x0 − x1 | > sl

(6)6

S. Marcello, Y. Shunra and M. Ruggero

Physica A 618 (2023) 128627

Fig. 4. Depiction of a network in 2D and 3D. The comparison between the two networks shows how by increasing the dimension, the number ofpaths with the same mMd increases accordingly. The presynaptic activity, increases as well in the network at a higher dimension, rising the chanceof the SGE effect.
In panel(a), two-dimensional grid, there are three different paths from nin to no with mMd = 3. In this situation we may have at most two presynaptic
signals s1 , s2 reaching no at almost the same time and triggering the SGE effect.
In panel(b), the three-dimensional grid induces 6 different paths with the same mMd = 3 between nin and no . In this case, the larger number ofpaths may let converge at most 3 presynaptic signals at around the same time, s1 , s2 , s3 , to no , increasing the probability for the SGE effect tobecome relevant.

{
Y =

{
Z =

|y0 − y1 |; if : |y0 − y1 | ≤ slln − |y0 − y1 |; if : |y0 − y1 | > sl

(7)|z0 − z1 |; if : |z0 − z1 | ≤ slln − |z0 − z1 | ; if : |z0 − z1 | > sl

(8)

Here X, Y, Z represent the least number of steps on each dimensional axis to reach no from nin , which coincides withthe definition of the minimum Manhattan distance mMd. Therefore, we can write mMd = X + Y + Z . For example, let
us assume that after our calculation we have (X, Y, Z ) = (3, 4, 5). It follows mMd = 3+4+5 = 12. It means that the least
number of steps from nin to no is 12. Starting from nin , we may move firstly 3 steps over the x-axis, 2 steps over the y-axis,
5 steps on the z-axis, and finally the last 2 steps again on the y-axis. That would be an example of the path from nin to
no with mMd = 12. We are now interested to know the number of all the paths with mMd = 12. Intuitively, the solution
may come by permuting the partitions of X, Y, and Z, and precisely is given by the multinomial coefficient [24]:P =(X + Y + Z )!

(9)

X !Y !Z !

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 thesame mMd to increase accordingly (as in Fig. 4).
2.3. Process time τ
We investigated so far how the SGE can be boosted by the number of paths between two neurons. We investigate now
the existence and relationship of this phenomenon upon the definition of the process time τ .
Specifically, we defined:7

S. Marcello, Y. Shunra and M. Ruggero

Physica A 618 (2023) 128627

Fig. 5. 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 increasingi,j
the randomness results in a wider range of variability of µf ,fp from 0 ms to 50 ms, with a central value of 25 ms. In this case, µ can become smalleri,jthan the processing time τfp , and positively influence the phase space of the neuron fs dynamic, provoking a faster postsynaptic response. In panel
4, for a small randomness level, the variability of µ never gets below the value of τ and the SGE effect will not take action.

• Process time τnif,j : the interval from the arrival of a pre-synaptic spike to neuron nf , the processing of this input andthe resulting emission of a post-synaptic spike. Here, i and j represent the loops of the central delay (cd)j and noise
delay (nd)i .
• The axonal delay µfi,,jfp : is the time for a spike to go from neuron nf to one of its adjacent neuron nfp , explained asin Eq. (4).
• The total transmission delay is: dfi,,fjp :
i,j

i,j

i,j

df ,fp = µf ,fp + τf

(10)

In Fig. 5 we represent schematically how the SGE effect works. The processing time τ (red line) and the transmission
time µ (blue line). We theorize how the existence of the process time τ allows the effect to take place. While thetransmission time µ, for any given central delay, represents the time for crossing the distance between two neurons,τ represents the time the neuron takes to emit the post-synaptic spike once the pre-synaptic stimulation has arrived. Thevery existence of this mechanism allows hypothesizing that while a previous signal is being processed, (PSPs)B in nf in
the plot, in the interval τ , another spike can reach up from nf, prompting a faster emission and causing the SGE effect totake place.
2.4. Single spike propagation
To better understand the characteristics of process time τ , we simulate a single impulse spike in a 3D grid of RS neurons.
The plot in Fig. 6 represents a simulation of 1000 ms in a 3D lattice grid network of (7 × 7 × 7) RS Izhikevich neuronsof a single spike crossing the fastest path from the input nin :12 to the output no :155 (the neuron number represents its
position in the grid). In this simulation, the initiator neuron emits a single spike at t = t0 which stimulates its neighbor’s
neurons. This action produces post-synaptic spikes that become pre-synaptic signals toward adjacent neurons and theprocess repeats forward across the grid.
Considering as spikes all impulses in Fig. 6 of positive membrane potential, we notice that the second spike (occurringat the first post-synaptic neuron 61) happens approximately at time 100 ms. The next adjacent neuron 154 spikes aroundtime 110 ms. The first spike does not fully propagate to the six adjacent neurons. Instead, its potential fragments insmaller PSPs along the adjoining connected neurons spread along the system. A single neuron may be reached by at most
six PSPs at approximately the same time and, eventually, provoke a spike (green dot). Noticeably, a single impulse spike
can propagate and create the remarkable action potentials across the whole network in a 1000 ms time frame.
8S. Marcello, Y. Shunra and M. Ruggero

Physica A 618 (2023) 128627

Fig. 6. The membrane potential (vertical axis) of 6 adjacent neurons representing the propagation of a single spike in a 2D network of regularspiking (RS) Izhikevich neurons. The central delay is 10 ms with a 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 duration near the average axonal delay value (cd = 20 ms in the simulation shown).

2.5. Spiking activity: RS vs heterogeneous model
Here, we analyze and compare the relationship between neuron-type randomness and the average number of spikes
occurring on the whole network over a period of time. We study the homogeneous and heterogeneous models on the 3Dmodel (4D torus) first. We run several simulations as the heterogeneity level varies from 0 to 1. Level 0 corresponds tohaving all neurons of the same type (RS Izhikevich), while, level 1 is a fully heterogeneous network of Izhikevich neurons
as defined in Ref. [10]. For all simulations, the input (initiator) nin :12 is triggered by a constant current of 10 mA. Knowingthat by increasing the axonal delay heterogeneity the spikes propagate faster, a logical prediction is to expect the numberof spikes to increase as the heterogeneity increases. Nevertheless, by interpreting Fig. 7, we encounter a counter-intuitiveeffect as the number of spike decrease as variability rises. In the plot, we have the results for 15 different simulations
9S. Marcello, Y. Shunra and M. Ruggero

Physica A 618 (2023) 128627

Fig. 7. The number of spikes against axonal delay randomness (range 0–20 ms), when the central delay is 22 ms, for 6 values of neuronal
heterogeneity and for the three different network dimensions. The number of spikes and the standard deviation is calculated over hundreds ofneurons, depending on the size of the network. The blue dot (1) indicates full heterogeneity the yellow dot (0.0) zero heterogeneity (all RS neurons).
In all three networks, the number of spikes for H = 1.0 is almost double compared to H = 0.0. Noticeably, the number of spikes decreases asspatial variability increases and the effect is stronger for models at higher levels of heterogeneity. The panels on the right are the linear regression
coefficient (slope, top) for the six heterogeneity levels (from 0.0 to 1.0). The mean value of the number of spikes for each heterogeneity level (spikes,
center) and the standard deviation (SD, bottom) of the respective mean.

10S. Marcello, Y. Shunra and M. Ruggero

Physica A 618 (2023) 128627

Fig. 8. The panel compares the variations in the number of spikes for each central delay and for six levels of heterogeneity. For each value of central
delay cd, the randomness level (nd)i varies in 15 steps as in Fig. 7. The number of spikes is calculated across all the neurons of the network forthe duration of the simulation. Since there are 15 simulations of different randomness levels, each point is calculated by averaging up all those 15
simulations. We range cd from 1 to 71 ms. The six models have the same constant, current I = 10 mA and synapse strength w = 18. For maximum
neuronal heterogeneity, the model shows the most robust spiking activity, we call this phenomenon stochastic neuron enhancement (SNE).

of 1000 ms each for six levels of heterogeneity. For each simulation, the central delay is fixed at the value of 22 ms,
while δi varies from 0 to 20 ms. The first simulation reports the average number of spikes of the network for a central
delay of 22 ms and maximum variability of 0 ms. This means that all the connections between neurons have the sametransition time µ: 22 ms. Additionally, being all neurons of RS type, it implies that they require approximately the sametime τ to process action potential signals. In Fig. 6, by considering µ and τ to be constant, we can discern a quite smoothsynchronization of the signals in the network. Specifically, we can expect a situation where, for instance, the signals
s1 , s2 , s3 of Fig. 4(b), will reach no at the same time, since µ and τ are constant. The synchronization leads the PSPs tocluster together into no approximately at the same time and, eventually inducing a spike. On the other hand, by increasing
the axonal delay randomness, the transmission time µ changes accordingly to Eq. (4), ergo, the signals loose synchrony
and arrive at different times to no lowering the possibility for the PSPS to add up and forming a spike. By increasing thenoise, synchronization decreases and so does the number of spikes.
Interestingly, we notice while the behavior of the network’s spiking activity maintains the same shape (Fig. 7) as
central delay changes in value, the distribution rises by increasing the heterogeneity (Fig. 8). All the networks deliver
a very different amount of spikes by varying the neuronal heterogeneity for each simulation. Specifically, while theheterogeneous axonal delay influences the distance between neurons, the neural heterogeneity is intrinsically related
to the neuron model, acting also on the variability of τ , namely, how different neurons would process a pre-synaptic
signal with different process time.
3. Results and discussion
We now compare the evidence of the SGE effect on three different networks and confirm the mathematical formulation
expressed in Eq. (9).
To verify the SGE effect from one dimension to another, the choice of nin and no must be done accordingly. As pictured
in Fig. 4, in panel (b) we chose no to be in a different plane compared to its counterpart in panel (a), where, even if themMd is 3 in both cases, the number of paths duplicates in the 3D grid model.
11S. Marcello, Y. Shunra and M. Ruggero

Physica A 618 (2023) 128627

Fig. 9. Comparison of regression’s slope for the three-dimensional toroidal networks 3D (blue), 4D (red) and 5D (green). The simulations show
the variation of the central delay (cd)j from 1 to 71 ms at 2 ms steps. For each delay value, we execute 40 simulations. By changing the spatial
randomness of the grid, we measure the time of arrival of the first spike. These 40 values are linearly interpolated, and the negative value of theslope reveals the SGE effect; in other words, it shows that the interval between the first spike of nin and the first of the output no is decreasing
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 have synaptic strength w = 18. we used the Izhikevich model as implemented in pyNest. For each network grid (2D, 3D and 4D), we chose(nin ; no ) to have the same Manhattan distance mMd. According to Eq. (9), by maintaining the same mMd and increasing the network dimension,
the number of paths with the same mMd will increase and strengthens the SGE effect.

We consider the three network grids (2D, 3D, and 4D) of regular spiking (RS) neurons. First, we examine thesimulation’s data for the spike timing. Depending on the dimension of the model, the heterogeneous spatial structureof the grid reduces the time of propagation of the spiking information.
We changed the intrinsic axonal delay (cd)j from 1 up to 71 ms and tested twenty different levels of noise delay (nd)ifor each (cd)j . In Fig. 9 we compare the regression 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 value
for the central delay, while the (nd)i varies from 0 to cd.
The negative value indicates how the signal propagates faster. The three networks have the same initiator neuron nin inposition 12, driven by a constant current of 10 mA. Across the whole network all neurons are connected with a constantsynaptic 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(9). Respectively for the network of 2D (Torus 3D), 3D (Torus 4D), and 4D (Torus 5D), no is set up to be in position 116,
155, and 296. This is resulting, accordingly, in 6, 20, and 180 number of paths with the same mMd = 6, for each model.
Noticeably, the average value of the slope decreases as the dimension rises, from −1.09 (3D Torus) to −2.05 (4D Torus)and −3.42 (5D Torus). This proves that the information propagates faster in networks of higher dimensions in agreementwith our discussion leading to Eq. (9).
In other words, increasing the dimension enhances the number of paths from nin to no , and increases the number of
spikes arriving at no , resulting in a higher probability for the SGE effect to occur.
4. Conclusions
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 on borderlessareas. In the first part of the manuscript, we generalized what was found in previous studies related to the influence ofheterogeneous axonal delay at various dimensions, so as to secure the effect of facilitating the transmission of information
within a model neural 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 signal to the adjacent nodes on the grid. The diffusionof the signal varies according to the grid structure, based on a von Neumann neighborhood model in 2D, 3D, and 4D
dimensions. In all the simulations undertaken in this study, for all the three different dimensions of the network, the12

S. Marcello, Y. Shunra and M. Ruggero

Physica A 618 (2023) 128627

Fig. 10. A comparison of the SGE effect over two models of neural heterogeneity (H = 1, blue) and regular spiking model (H = 0, R.S., red). According
to Fig. 8, the H = 1 line shows the greatest amount of spiking activity, while the R.S. is the minimum. Nevertheless, the effect acts in a similar
measure on both networks. All the models have the same constant current input of 10 mA, a synaptic strength of 18, and the same number ofneurons. Noticeable, the heterogeneity offers the most suitable model regarding robustness (highest spiking activity) while keeping the benefit of
improving spike propagation times (SGE effect).

negative value of the regression line’s slope m, Fig. 9, indicates how heterogeneous axonal delay (noise) induces a fasterpropagation 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 called the stochastic grid enhancement (SGE) effect. We extended
the properties of this counter-intuitive result, generalizing the concept for any network dimension, showing that the
information contained in spiking activities travels faster under higher heterogeneous axonal delay randomness. Throughan analytical demonstration and several simulations, we have shown that at any dimension the increase of 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 neural type heterogeneity for the same set of spiking neural networks bymanipulating the parameters c, d of the Izhikevich model throughout a uniformly distributed random variable. The level ofneural heterogeneity ranges from 0 to 1, with 0 representing a perfectly homogeneous regular spiking (RS) neural model
network and 1 a heterogeneous one as defined in Eqs. (3). In all the cases, we set synaptic weights at w = 18, and theinitiator neuron nin input current at 10 mA.
Recognizably as shown in Figs. 8 and 10, by keeping constant the initial value of current, synaptic strength, and thenumber of neurons of the model, introducing cell type heterogeneity results in a large increase in spiking activity aphenomenon that we called here stochastic neuron enhancement, SNE. The information created from the spiking activity
was abundant in heterogeneous networks since the networks were able to increase spiking activity in the same period oftime and with equal beginning values for the number of neurons, current, and synaptic strength. As a result, we conclude
that the combination of neuronal heterogeneity and heterogeneous axonal delay 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 neural cell heterogeneity is therefore a biologically and computationallyeffective strategy for the brain to improve signal robustness.
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 andheterogeneous axonal delay appear to favor computational processes enhancing perception and other performance inhuman subjects. While in this study we kept the parameter manipulation to a minimum to better clarify the effectof heterogeneity on the information transmission, further statistical analysis such as synchronization factors, discharge
13

S. Marcello, Y. Shunra and M. Ruggero

Physica A 618 (2023) 128627

frequency and other possible statistical analysis as suggested to deepen the understanding of spiking informationtransmission.
CRediT authorship contribution statement
Salustri Marcello: Conceptualization, Methodology, Software, Formal analysis, Writing – original draft. Yoshida
Shunra: Data curation, Software. Micheletto Ruggero: Validation, Supervision.
Declaration of competing interest
The authors declare the following financial interests/personal relationships which may be considered as potential
competing interests: Marcello Salustri reports financial support was provided by Yokohama City University.
Ethical approval: This article does not contain any studies with human participants performed by any of the authors.
Data availability
Data will be made available on request.
References[1] N. Perez-Nieves, V.C.H. Leung, P.L. Dragotti, D.F.M. Goodman, Neural heterogeneity promotes robust learning, Nat. Commun. 12 (5791) (2021)
http://dx.doi.org/10.1038/s41467-021-26022-3.
[2] C. Koch, G. Laurent, Complexity and the nervous system, Science 284 (1999) 96–98, http://dx.doi.org/10.1126/science.284.5411.96.
[3] P. Dayan, L.F. Abbott, Theoretical neuroscience: Computational and mathematical modeling of neural systems, J. Cogn. Neurosci. 15 (1) (2003)
154–155, http://dx.doi.org/10.1162/089892903321107891.
[4] B. Gutkin, D. Pinto, B. Ermentrout, Mathematical neuroscience: From neurons to circuits to systems, J. Physiol.-Paris 97 (2) (2003) 209–219,
http://dx.doi.org/10.1016/j.jphysparis.2003.09.005.
[5] S.G. Tewari, M. Gottipati, V. Parpura, Mathematical modeling in neuroscience: Neuronal activity and its modulation by astrocytes, Front. Integr.
Neurosci. 10 (2016) http://dx.doi.org/10.3389/fnint.2016.00003.
[6] A.A. Faisal, L.P.J. Selen, D. M. Wolpert, Noise in the nervous system, Nat. Rev. Neurosci. 9 (2008) http://dx.doi.org/10.1038/nrn2258.
[7] Z. Budrikis, Forty years of stochastic resonance, Nat. Rev. Phys. 3 (771) (2021) http://dx.doi.org/10.1038/s42254-021-00401-7.
[8] M. Salustri, R. Micheletto, Heterogeneous axonal delay improves the spiking activity propagation on a toroidal network, Cogn. Comput. (2022)
http://dx.doi.org/10.1007/s12559-022-10034-2.
[9] W. Stacey, D. Durand, Stochastic resonance improves signal detection in hippocampal CA1 neurons, J. Neurophysiol. 83 (3) (2000) http://dx.doi.org/10.1152/jn.2000.83.3.1394.
[10] E. Izhikevich, Simple model of spiking neurons, IEEE Trans. Neural Netw. 14 (6) (2003) 1569–1572, http://dx.doi.org/10.1109/TNN.2003.820440.
[11] B. Vazquez-Rodriguez, A. Avena-Koenigsberger, O. Sporns, Stochastic resonance at criticality in a network model of the human cortex, Sci. Rep.
7 (2017) http://dx.doi.org/10.1038/s41598-017-13400-5.
[12] A. Tozzi, F.P. James, Towards a fourth spatial dimension of brain activity, Cogn. Neurodyn. 10 (3) (2016) http://dx.doi.org/10.1007/s11571-0169379-z.
[13] T. Suter, A. Jaworski, Cell migration and axon guidance at the border between central and peripheral nervous system, Science 365 (6456)(2019) http://dx.doi.org/10.1126/science.aaw8231.
[14] B. Sonal, G. Mahendra, Design of mesh and Torus topologies for network-on-chip application, Int. J. Reconfigurable Embed. Syst. (IJRES) 2 (2013)
http://dx.doi.org/10.11591/ijres.v2i2.3434.
[15] C. Wu, Y. Li, S. Chai, Design and simulation of a Torus topology for network on chip, J. Syst. Eng. Electron. 19 (4) (2008) 694–701,
http://dx.doi.org/10.1016/S1004-4132(08)60141-3.
[16] D.A. Zaitsev, A generalized neighborhood for cellular automata, Theoret. Comput. Sci. 666 (2017) 21–35, http://dx.doi.org/10.1016/j.tcs.2016.11.
002.
[17] J. von Neumann, Probabilistic logics and the synthesis of reliable organisms from unreliable components, Automata Stud. 34 (1956) 43–99.
[18] J. Eppler, M. Helias, E. Muller, M. Diesmann, M.-O. Gewaltig, PyNEST: A convenient interface to the NEST simulator, Front. Neuroinform. 2 (3)(2009) http://dx.doi.org/10.3389/neuro.11.012.2008.
[19] S. Chaturvedi, M. Boudjelal, A. Khurshid, Comparison of LIF and izhikevich spiking neural models for recognition of uppercase and lowercase
English characters, CiiT Int. J. Digit. Image Process. 6 (2014).
[20] S. Nobukawa, H. Nishimura, T. Yamanishi, The importance of neighborhood scheme selection in agent-based tumor growth modeling, Sci. Rep.
8 (2018) http://dx.doi.org/10.1038/s41598-017-18783-z.
[21] M. Coli, P. Palazzari, R. Rughi, The toroidal neural networks, in: 2000 IEEE International Symposium on Circuits and Systems, Vol. 4, ISCAS,
2000, pp. 137–140, http://dx.doi.org/10.1109/ISCAS.2000.858707.
[22] M. Madadi Asl, A. Valizadeh, P. Tass, Dendritic and axonal propagation delays determine emergent structures of neuronal networks with plastic
synapses, Sci. Rep. 7 (2017) http://dx.doi.org/10.1038/srep39682.
[23] R. Suwanda, Z. Syahputra, E.M. Zamzami, Analysis of Euclidean distance and Manhattan distance in the K-means algorithm for variations
number of centroid K., J. Phys. 1566 (2019) 696.
[24] M. Sato, T. Sado, Lattice paths restricted by two parallel hyperplanes, Bull. Inform. Cybern. 21 (3/4) (1985) 97–105, http://dx.doi.org/10.5109/
13371.

14