Cognitive Computation
https://doi.org/10.1007/s12559-022-10034-2

Heterogeneous Axonal Delay Improves the Spiking Activity
Propagation on a Toroidal Network
Marcello Salustri1 · Ruggero Micheletto1
Received: 3 December 2021 / Accepted: 26 May 2022
© The Author(s), under exclusive licence to Springer Science+Business Media, LLC, part of Springer Nature 2022

Abstract
Several studies have looked into how noise affects neural networks and actual brains as evidenced by transcranial random
noise stimulation, which improves cognitive performance. This research aims to broaden this understanding by concentrating on the network structural heterogeneity realized by adding noise to a neural model network’s axonal propagation delay.
We utilized the pyNEST neural network simulator to model a network of 400 artificial Izhikevich neurons connected by afolded von Neumann neighborhood to form a toroidal shape where axonal propagation noise simulates a variable spatial
spacing between neurons. In this network only one neuron is regularly spiking at first because it is specifically stimulated by a 10mA
external current, while all the other neurons have no external input and are stimulated solely by the activity of their neighbors. The
forward propagation of the spiking wave from the original neuron to its neighbors, and then to distant nodes on the toroidal
network, was investigated. For each simulation, we recorded the activity of all the network changing several parameters toverify differences of spike activity in different positions on the torus. By manipulating heterogeneity, we discovered that
adding noise helps the signal reach distant neurons in 20% less time, compared to when there is no heterogeneity. We demonstrated for the first time that structural heterogeneity in a neural network can favor the propagation of spiking waves. This
result is in line with other findings that suggest that a certain level of noise is good for the brain, extending this concept to
the network physical structure.
Keywords Noise delay · Spike · Stochastic resonanceIntroduction
Modern neuroscience copes with the considerable development of mathematical models describing the biological functioning of dynamical systems resembling the brain. To thelarge amount of brain-related data accumulated over time,
there has been a proportional engagement in mathematicalcomputing simulation to repeat and match the experimental
results [1]. Furthermore, the counterintuitive effect of stochastic resonance plays a role in understanding how noise
improves the performance of a system instead of reducing
it [2].
Supported by the Yokohama City University* Ruggero Michelettoruggero@yokohama-cu.ac.jp
1

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

In a larger perspective, the brain is a complex system
whose signals transmissions take place in a noisy heterogeneous environment [3]. Several studies show how the presence of stochastic resonance provokes substantial improvements in signal detection [4]. In this paper, hence, we wantto contribute to this field, studying the spiking activity propagation in a model network and finding a simple mathematical model explaining the mechanism for which the presenceof heterogeneity [5] can improve the propagation of a signal
within neurons. Noticeable examples are studies on transcranial random noise stimulation (tRNS), where subjects arestimulated by large electrodes with weak random electric
stimuli that actually improve their different motor, sensory,and even cognitive tasks [6].
In this paper, we focus on another type of noise: the heterogeneity of inter-neuronal distance. In a real brain, neurons
are not placed in a regular lattice, like in a crystalline structure. Instead their position is affected by a certain degree ofrandomness. Since in a model network, the distance between
neurons is represented by the propagation delay, in this study

13

Vol.:(0123456789)

Cognitive Computation

we identify how the increasing of noise in axonal propagation time decreases the time delay of signal transmission
through a simple network model of 400 neurons. Particularly, using an Izhikevich neuron model implemented on a
toroidal framed network [7, 8], we investigate the different
time delays of a signal traveling across the network while
the system is undertaken to different heterogeneity levels
through axonal delay manipulation. The heterogeneity has
been implemented by using a uniformly distributed noise
parameter. Neurons in the brain deal with information in
confined areas, which generally are named and classified
depending on their function. These areas, however, do nothave defined ends or borders. In this sense, the choice of atorus instead of a confined layer or plane is natural becauseof its border-less nature. The utilization of toroidal topology in previous studies [9], together with the mathematical topological nature of the torus, which we might want to
exploit for future research, led us to frame our network in
the toroidal shape.
We find that the increase of axonal delay noise in the
network corresponds to a decrease in the total propagation
time of the signal from an input defined as initiator nin , toan arbitrary test neuron called output no. Communication ofneurons takes place through action potentials (spikes), which
travel to neighbors neurons inducing other spikes. Overall
this creates a wave of spikes that propagates along the network. The Hodgkin–Huxley model, introduced in 1952 andawarded with Nobel Prize in 1963, mathematically describes
all the plausible action potential patterns properties of biological membranes with a large set of nonlinear differential
equations [10].
Over the years, the models have been extensively analyzed for developing methods that can be used for more
complex systems. An extraordinary compendium of the
Hodgkin–Huxley model was presented by Eugene M.
Izhikevich [11], who enhanced a model that accurately
resembles the spiking and bursting behavior of corticalneurons, combining the biological outcomes of their modeland the mathematical interpretation of integrate-and-fireneurons with a two-dimensional system of ordinary differential equations. Implemented by a PC program, Izhikevich
could reproduce around twenty different responses whichcharacterize real neurons spiking behavior.
As mentioned above, our model network lies on a toroidal
structure to emulate confined neuronal domains. Moreover,the symmetric properties of this topology give us the advantage to have always four neighbors for each node and thatmakes the network a perfect n-dimensional grid.
Specifically, each node is connected to its 4 nearest neighbors and corresponding nodes on opposite edges are connected in a cyclic fashion. The communication of each node
can take place in 4 directions, +x or east, -x or west, +y or
north, and -y or south. A torus network can be defined as a

13graph G = (N, C), where N(G) and C(G) are, respectively,
the nodes and connections of G [12]. The total nodes of the
2D torus are n2 , and the two-dimensional torus structureleads naturally to a von Neumann-shaped neighborhood[13], used to define a set of cells surrounding a given cellprovided with four inputs.
We initiate the simulation by choosing a neuron (corresponding to a node in our toroidal network) that acts asa stimulus and that we name “initiator.” This is connected
to a constant external current of I = 10mA to induce it tospike regularly for the total duration of the simulation (inmost of our tests 1000 milliseconds). Each single node isconceived to possess an intrinsic axonal delay d = cd + 𝛿 ,where cd is called from now on “central delay” and𝛿 = nd ⋅ 𝛼 . The parameter nd is called “noise delay” whereas
𝛼 = [(2x) − 1]; (0 ≤ x ≤ 1) is a random variable with uniform distribution used to implement the axonal propagation
delay heterogeneity. We, hence, analyze the propagation of
the signal within the network up to the chosen output. To
gather more interesting data, we chose the two input-output
neurons to have the maximum distance between them in our
structure.
Given the symmetric toroidal shape of the network andconsidering the model to be framed over a 20×20 grid of
400 neurons, the maximum distance from two neurons is a
Manhattan distance of 20, to which corresponds 10 nodes
along the diagonal [14].
At the beginning of the simulation, only the initiator neuron is spiking, then the spike activity propagates up to a
far away test neuron that we name the “output” node. The
simulation calculates the time it takes for the first spike toreach the chosen output node. We run several simulations,in each of which we increase the noise delay using the randomly uniform distribution of the parameter x of 𝛼 . Since
the dis-homogeneity is random and uniform, we should not
expect any advantage due to the noise, instead, we will provethat it favors a faster propagation of spiking activity and wewill give a simple theoretical framework to understand why.

Methods
To implement our simulation, we used Python and the NEST
library version 2.2 [15], which enabled us to portray a neuraldynamic experiment by creating and connecting neurons tosimulate a spiking network in continuous time. With the NEST
library, we were able to implement specific functions to large
sets of nodes whose connections have a configurable delay/
weight as well as parameters and state variables. The libraryfunctionality to implement loops and integrate the Izhikevich
differential equations, simplifies the process of software development [16].

Cognitive Computation

In our structure, we used the neuron model developedby Eugene Izhikevich. Among several models for studying spiking neural networks (SNN), we use the model proposed by Izhikevich because of its high plausible biologicaldynamics, fundamental to reproducing well the non-linear
phenomena dynamics of the neural network. Without realistic neuronal dynamics, we will not be able to explain thefaster spiking activity propagation caused by the increased
spatial heterogeneity in the network. Furthermore, by properly setting a few parameters, the Izhikevich neuron model
reproduces most of all spiking activities, offering a solid
accuracy to reproduce spike patterns [17]. Concerning the
propagation delay between neurons, in NEST, it is implemented by a pipeline/buffer procedure. At each time step, aRunge-Kutta engine processes the Izhikevich model differential equations in which the input of the presynaptic neuron is taken into account by a weighted sum (the weights
represent the synaptic connection strength parameter).
However, if a neuron is delayed by a certain time, its signal
circulates in a memory register, a pipe, long as the delay
is measured in time steps number. In other words, at eachtime step, the presynaptic membrane potential 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 timedelay between presynaptic and post-ynaptic neurons in aseamless manner without interfering with the theoreticalmodels of the delays in this study.
The 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 characterize both the spiking and the bursting
behavior of the neurons [11]{ �
v = 0.04v2 + 5v + 140 − u + I
u� = a(bv − u)with after the spike, a reset process is given by:{v←cif v ≥ 30mV →
u←u+d

(1)

, and u� = du.
where v� = dv
dt
dt
The variables v and u and the parameters a, b, c and d
are all dimensionless [11], t is the time. In the mathematical model, v represents the membrane potential and u amembrane recovery variable, providing the negative feedback to v. After each spike, according to Eq. (1) the membrane voltage and the recovery variable are reset down tothe parameter c and to u + d , respectively. I is the variable representing the external current stimulus, while theadjustment 0.04v2 + 5v + 140 allows the membrane potential v and time t to be scaled, respectively, to mV and ms.

Fig. 1  There are several shortest paths to go from A to B, but the
Manhattan distance, which in this case 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

The parameters a, b, c and d, according to their values,
enable the equation to depict various firing patterns model.
In particular:a ∶ (0.02 msec−1 ) represents the recovery time for u; toan increment of a, resulting in a quicker recovery for u.
b : (0.20)[dimensionless] underlines the strength ofu to the sub-threshold fluctuations of the membranepotential v. An increment of b is translated into astronger interrelation between u and v, meaning possible sub-threshold oscillations and low-threshold spiking
dynamics.
c ∶ (−65mV) represents the after-spike reset value of themembrane potential v.
d : (6mV) describes the after-spike reset of the recovery
variable u.
The model of Izhikevich empowers us to define each neuron with an accurate and biologically feasible mathematicalstructure, leaving us with the only need to establish a plausible connection structure for each node/neuron of our network.
A reasonable choice to design the dynamic network,
find its natural place in the two-dimensional lattice von
Neumann neighborhood (4-neighborhood model as inFig. 2), where each cell is connected with its four adjacent cells located on its north, south, east and western
position [13]. In machine learning, the Manhattan distance(or taxicab geometry) is used to calculate the distance (or“closeness”) between two points/nodes. Contrary to the
Euclidean distance, there may be several different paths
with the same Manhattan distance between two points, as
shown in Fig. 1.

13

Cognitive Computation

Fig. 2  Toroidal network of 16 neurons plotted in 4 × 4 square matrices. To create a biologically plausible continuous system, which is to
avoid the presence of boundaries or borders in the structure, the firstrow is connected to the last row and the first column to the last one

Once we have defined the new metric, we generate a twodimension matrix of 400 neurons connected according to the
von Neumann neighborhood model and consider the firstrow to be up-connected with the last row (and vice-versa),
and the first column to be left-arrow connected with the lastcolumn and vice-versa (see Fig. 2 for a schematic exampleof the structure). In this way, by theoretically bending theedges, we reshape the two-dimensional matrix in a treedimension torus, on the surface of which lays our dynamic
network that represents a model of a neural aggregate where
information is confined into. This leads to a resembling affinity with the toroidal neural networks (TNN) [18], defined by
2D toroidal topology. The toroidal structure addresses thestudy to several appealing geometric interpretations.
For instance, as we can see in Fig. 3, the developmentof spike rate on a toroidal surface tends to approach thegeodesic line, a locally length-minimizing curve. From a
topological perspective, the torus shows important homeomorphic properties, empowering the research to deepen
further studies in computational topology for improving theunderstanding of brain functionality [19, 20].
Once we framed the network in a toroidal model anddefined the neurons according to Izhikevich’s archetype,we run the simulations. In the first part of the program, webuilt a matrix whose nodes represent the neurons, connected
in a pattern to form the grid of toroidal structure with 400
neurons (20 × 20 arrays). To maintain the system in its simplest form all neurons are considered identical, there are
no inhibitory synapses and no thalamic currents. We want
only to focus on the effect of dis-uniform axonal propagation
delay [21], so we study the model at its prime conditions. In
this, we considered RS (regular spiking) neurons, cause they
address the most typical neuron class in the cortex.
We use the following terminology to identify the delay andstrength connection for a neuron. Each of them is here characterized by its synaptic strength connection w, which is constant
all along with the network and will not be considered in our
calculation.

13

As introduced above, the intrinsic axonal delay for eachneuron is, in general, defined by the variables d = cd + 𝛿
where cd represents the time it takes for the spike to propagate
along the axon and 𝛿 it is a random variation around this value.
We have two main loops. The external one runs the integer jthat varies the neuron central delay cd, whereas the internalloop runs i that changes its stochastic noise delay part 𝛿. We
initially assigned a value to the variable ns (for noise steps)and (cd)max , representing respectfully the number of i-loop
and j-loop such that:1 ≤ i ≤ ns , and 1 ≤ j ≤ (cd)max , i, j ∈ ℕ
When i reaches ns, after ns i-loops, j increases of 1. For simplicity, we consider the value of the j-th central delay (cd)jto be equal to j, so it assumes the integers values (1, 2, 3, 4...
(cd)max ). Therefore, this structure realizes a heterogeneous
intrinsic delay di,j for each neuron. The variable parameter
is formulated as follows:

di,j = (cd)j + 𝛿i,j = (cd)j + (nd)i,j 𝛼i,j

(2)

where
• 𝛼i,j = [(2xi,j ) − 1]
• (cd)j = central delay; 1 ≤ (cd)j ≤ (cd)max , j ∈ ℕ• 𝛿i,j = (nd)i,j ⋅ 𝛼i,j ; (𝛿i,j ≤ (cd)j )• for each j: (nd)i,j = i(cd)j, 1 ≤ i ≤ ns, i ∈ ℕ where0 ≤ xi,j ≤ 1 is an aleatory variable with uniform distribution assigned to each neuron at the i-loop and j-loop. In thisway 𝛼i,j is always within the −1,+1 interval.
ns

The simulation starts with j=1 (cd=1) which, at each i-loop,
assigns to each neuron an intrinsic delay di,j , whose value,according to Eq. (2), ranges as follows:(cd)j − (nd)i,j ≤ di,j ≤ (cd)j + (nd)i,j
For any single (cd)j , we run ns number of i-loops (ns-simulations) and get the value for (nd)i,j as follows:(nd)i,j = i

(cd)j
ns

, 1 ≤ i ≤ ns, i ∈ ℕ(3)This framework realizes a network of 200 neurons, eachof them having a different (heterogeneous) intrinsic axonal
delay di,j . In the whole network, only the special neuron nin is
stimulated by an external current I = 10mA. Once this starts
to spike regularly, its activity propagates to its neighborswith a certain delay. Then its neighbors once activated, willtransmit their activity to others and so on. To measure thespike activity propagation speed between two defined neurons, a neuron distant from nin is used to test the arrival timeof the spike information coming from nin itself. This neuron

Cognitive Computation

Fig. 3  In A1 and B1 we run three simulations on the toroidal network
of Fig. 2. Neuron 3 is considered the input of the network, being theonly one connected to the external input and to initially spike. The
three plots in A1 show the spike rate calculated in three different simulations of 1000 mseconds. In B1 other three simulations are shown,
but this time the input neuron is in position 1. On B1 the path 1-59-13-14-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 high spike frequency tends tolay on the path 3-7-11-15. Alternative paths with the same Manhattan
distance are 3-2-14-15 and 3-4-16-15. In the three-dimensional torus,
if we consider the Euclidean distance, the sequences in red in A2 andB2 (1-5-9-13) tend to lay on the longitudinal geodesic, which, eventually, would be interesting in terms of optimization of a network

is indicated with no . The subscripts in and o stand looselyfor input or initiator and output, respectively.
Consequently, we define:Accordingly, we can infer the spike propagation delay fromnin to no as:

•
•

i,j
fout as the time of the first spike of no.
i,jfin as the time of the first spike of nin.

i,j

i,j

Δfgi,j = fout − fin

(4)

13

Cognitive Computation

Fig. 4  The spike activity propagation delay on the torus. Here thefirst spike happens at time t = 0 at the node located at index 30, thedelay is calculated when the neuron at the node in the position indicated on the horizontal axis has its first spike. The tested neurons arethose located at the geodesic of the toroidal grid. There is a linear
dependence between nodes distances and propagation m = 2.5 msec/
neuron represents the inclination of the linear regression. In thesetests the intrinsic delay and delay heterogeneity are kept constant at
cd = 50 msec and 𝛿 = 12.5 msec, respectively

where g represents the minimum number of steps from ni
to reach no , that is, the minimum Manhattan distance [18]between the two neurons.
i,jFor each i, therefore, we identify a specific Δfg , and fora whole i-loop we gather ns values for it.
i,j
All the Δfg of one i-loop have the same central delay(cd)j = j . When i reaches ns, the central delay (cd)j is incremented of 1 and starts another cycle of i-loops.
As already stated, when j reaches (cd)max , the simulai,jtion stops, leaving us with ns × (cd)max values of Δfg to beanalyzed.
From the data, we can infer if and how the increase ofnoise delay (nd)i,j has affected the transmission of the spikesignal. More accurately, we want to ascertain if the increasei,jof (nd)i,j decreases Δfg . For this purpose, in the next section
we consider to compare the regression line for each i-loop.

spike can be also called time of arrival in the sense that it
represents the time in which the information from nin arrives
at another neuron. As expected, we have a linear proportionality between the Manhattan distance between the test
neuron and the initiator. In our specific case, the network
is a grid of 20 × 20 neurons. Thus, identifying the neurons
with an index that counts each of them along the grid, andsetting the initiator in the center of the second raw, gives theinitiator the index of 30 (20 neurons for the first raw, plus
10 neurons for the middle position in the second raw). With
this indexing convention, we test the time of arrival of the
spike-wave, for each neuron along the geodesic toward theoutput neuron 230 (placed as well on the geodesic line, on
the central raw).
As noticeable from the simulation results of Fig. 4, thepropagation delay in milliseconds is growing linearly with
distance as intuitively we expect.
So, considering the time of spike arrival for each neuron
in the network, we could elucidate if and how noise improves
the circulation of spike activity in a robust and reproducible
way. In Fig. 5 are shown two maps representing the toroidal
network at two different noise levels. Each block represents
one of the neurons and the color represents the time of thefirst spike in milliseconds. Since the only neuron that is
stimulated externally is the initiator located in position x =10, y = 2 in the grid (its index is then 30), all the neurons
that are spiking, do so because of spike axonal propagation
along the network. When the noise level is lower (Fig. 5A,
ndi,j = 5msec) the number of neurons never reached by the
spiking activity is larger (deep purple color in the map).
Instead when the noise level is higher (Fig. 5B, ndi,j =31msec), the opposite phenomena is observed. Overall
speaking, this leads to conclude that increasing the noisei,jwill cause the total propagation delay df ,f of the network topdecrease. If we consider Eq. (4), we can understand that byincreasing the noise parameter (nd)i,j (that is when i is big),i,je,jthe probability to have Δfg < Δfg for most simulations
increases. Namely, the regression line will be negative, thatis, m < 0.
To describe the effect of noise on the delay, we reckon thegeneral equation for the regression line [22]:ȳ = m� x̄ + 𝛼�Results and Discussion
To explain the effect of noise on the delay, firstly we observe
how the spike activity propagation time evolves along with
the network when it has a constant value of delay heterogeneity. In Fig. 4, with we pose (cd) = 50msec, with(nd) = (cd)∕4 = 12msec. Again, only neuron nin is spiking
initially, and in the plot we show the “time of first spike”
of all other neurons. Since the spike activity is propagatingfrom nin in all directions in the network, this time of the first

13

where
∑n

m
�=

(xi − x̄ )(yi − ȳ ) 𝜎x,y
= 2
∑n
𝜎x
(x − x̄ )2
i=1 i

i=1

Here we introduce• x̄ and ȳ as the average of the xi and yi.
• 𝜎x and 𝜎y as the standard deviations of x and y.

(5)

Cognitive Computation

Fig. 5  We show the spike delay for each neuron of the 20 × 20 toroidal network. For each neuron in the network, we report the first spike
occurrence’s time, representing the time of arrival of the spike information of the neuron that initially is spiking. This is the initiator neuron nin of index 30, located in x = 10 and y = 2 in the map, which isspiking from t = 0 due to a constant input current I0 = 10mA. The
color bar shows the time of the first spike for each neuron in milliseconds. The dark blue regions are where the neurons never spike,• 𝜎x2 and 𝜎x,y as the variance and covariance.

i.e., the spiking activity from the initiator nin doesn’t arrive in timein those locations (our test simulation is of 500msec). Noticeably, inpanel 1 the noise level is small (ndi,j = 5 msec) a good portion of
the network is never reached by the spike activity, whereas when
the noise level is higher (panel 2, ndi,j = 31 msec) the circulation of
spikes is faster with more neurons reached by the activity of the initiator

𝜎x2 =In our model (xi,j ;yi,j ) = ((nd)i,j ;Δf i,j ).
We now try to give a mathematical formulation for m andprove 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, wheni,jthe noise (𝛿)i increases, the delay Δfg decreases. From Eq.
(5), we calculate the slope for the linear regression1, where:∑ns(xi − x̄ )(xi − x̄ )2(6)𝜎x = i=1
ns

∑ns

i=1

(xi − x̄ )2ns

𝜎x,y =(x − x̄ )(yi − ȳ )
i=1 i
ns

(7)

our, respectively, values are:
•
•

(cd)
(cd)j (ns+1)j
xi = (nd)i,j = i ns j ; x̄ =
2ns
∑ns
i,j
Δf
j
yi = Δf i,j ; ȳ = i=1
ns

To simplify the symbolism from now on, we reckon
cdj = cd, and Δf i, j = Δf i . We omit the index j since cdj isconstant during the calculation of a single m while only iranges from 1 to ns. By substituting xi , x̄ , yi , ȳ respectfullyin Eqs. (6) and (7), we get:

1(cd)2 (ns2 − 1)
12
ns2

and
ns
ns
(ns + 1) ∑ i(cd) ∑iΔf ]𝜎x,y = 2 [ (iΔf ) −2
ns i=1i=1which we can write as:
ns−1

𝜎x,y = (cd)ns2

2
∑

i=1

∑ns

=

|(i −

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

From (5), we can formulate:m=∑ ns−12i=1

�(i − (ns+1))�(Δf (ns+1−i) − Δf i )22(cd) ⋅ (ns12−1)

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

(8)

The Eq. (8) demonstrated how, by increasing the index
i (i.e., noise), the delay of a signal is reduced comparedto 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.

It is just reported a part of the demonstration.

13

Cognitive Computation

Fig. 6  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 ofthe initiator nin. These times represent the times of arrival of the spiking activity from an initiator neuron placed at the position 30 to a testor output neuron placed far away in the toroidal network (in position230). We performed 40 simulations with increased noise levels, fromcd/40 msec up to a maximum of cd = 100 milliseconds. The resulti,jing histogram is centered around a delay difference Δfg of about 800
milliseconds. In the central panel, we plot the same arrival times differences in function of the noise (nd)i,j . We have 40 increasing valuesof (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, positioning itself
i,jon values greater than 800. Eventually, by increasing the noise, Δfg
decreases accordingly, reaching values below 800 msec. The dotted line represents the linear interpolation with an inclination ofm = −2.15, as indicated in the figure’s title. In the bottom panel, weplot all the values of m obtained by previous simulations, where eachone has a different neuronal intrinsic delay (cd)j (ranging from 5 to
100 msec at 1 msec steps). While the top and central panel refer only
to the value of (cd)j = 93.07, the bottom panel reports all the previous
values of m up to the current (cd)j . A negative m indicates that thevalue of delay is decreasing

In our study, we changed the intrinsic axonal delay (cd)jfrom 1 up to 100 milliseconds. For each of these values, wetested forty different levels of noises, ranging from (cd)j ∕40
up to (cd)max = 100msec. The top panels of Figs. 6 and 7
show the distribution of the spike arrival times in the output neuron no. In Fig. 6, the central intrinsic delay is large,
93msec, so only 20 arrival times are reported. Our simulation lasts 1000 msec, and depending on the variability of thestochastic 𝛿i,j factor, in some cases, the neuron no is nevertriggered to spike (the spike information from nin doesn’t
arrive in time). Instead, as in Fig. 7, the histogram is complete when the central delay ranges around 31msec. It has 40
arrival times, with the distribution concentrated on a shorter
delay (approximately 350msec). The central panel of thesefigures shows the time of the first spike at no for increasing(nd)i,j . The dark dotted line in the linear regression depicts

the significant characteristic of how such time decreases as
the noise level rises. Interestingly, examining the bottom
panel of both figures, the slope of the regression m is always
negative, with a minimum outcome around cd = 35 msec
(the y axis is reversed).
The bottom panel of Fig. 6 reports the value of m for each
cd, from the first to the last simulation (i.e., 2 ≤ cd ≤ 93.07).
We were able to reproduce the minimum value of m around
cd = 35 msec for repeated identical simulations with different random seeds. The peak of m was always located arounda global minimum between 30 ≤ cd ≤ 35.
The slope m, which indicates how noise reduces the
propagation speed within the network, remains negative inall simulations we’ve conducted. The fact that this occurs
in a repeatable manner demonstrates the robustness of thisfeature in our model neural network.

13

Cognitive Computation

Fig. 7  Same as Fig. 6. 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 testing for intrinsic delay values in the range 31 ≤ cd ≤ 35. When the
intrinsic delay is centered at 35 msec, the signal from nin propagates
fasterModel
We‘ll now try to give a theoretical interpretation of thenegativity of m.
To do so, we must investigate the neuron’s dynamic system characteristics, which necessitates an understanding ofFig. 8  Panels 1 and 2: a stable equilibrium. A small presynaptic signal (PSP) goes from neuron nf to nf1. The membrane potential willfluctuate around the equilibrium point before returning to its restingstate (black dot, the equilibrium point). Panels 3 and 4 depict the
responses of two different presynaptic stimuli. The small PSP (A)the membrane potential/recovery system’s phase portrait.
Conditions such as the equilibrium, limit circle, attractor [11], and other characteristics in the phase space help
to clarify the spiking behavior from a dynamical system
point of view.
We analyze here two phases, starting with the situation of a neuron in an equilibrium state, whose membranepotential is resting [23].
In Fig. 8(1), input currents produce a small amount ofpresynaptic potential (PSP) [24], in the phase diagram ofFig. 8(2) its membrane potential moves a bit around theequilibrium (black spot, attractor-resting potential) and,after few milliseconds, it will go back to its resting position.
Instead, in Fig. 8 panel 3, we show the effect of twopresynaptic stimuli. A little signal PSPs (A) causes a slightchange in the equilibrium, whereas a bigger signal PSPs
(B) causes the neuron’s intrinsic dynamics to spike after ashort time. The period between the arrival of the stimulus
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 postsynaptic neuron while that same neuron isstill firing as a result of a prior PSP signal B. Consider thecase where a spike leaves a neuron nf and travels in the
direction of neuron nf 1.
Firstly we define:• Process time 𝜏f i,j : that represent the time for neuron nf toprocess a presynaptic spike and deliver its own postsyn-

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

13

Cognitive Computation

Fig. 9  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 the axonali,jtransmission 𝜇f ,f is wide (from zero to 50ms in this example wherep
the central value of the delay is 25ms). Since this value can becomei,jsmaller than the processing time 𝜏f , this will positively influence thephase space of the neuron’s dynamic, resulting in a faster postsynaptic response. In panel 2 it is shown the opposite situation in which thevariability of the axonal delay is smaller and never realizes this condii,jtion. Since a longer than average, 𝜇f ,f doesn’t increase postsynapticp
process times, whereas a shorter than average does, this indicates how
a degree of randomness helps the overall circulation of spike activity
on the toroidal network

aptic spike, the two indexes i and j represents the current
central delay (cd)j and noise delay (nd)i,j.
• Axonal transmission delay 𝜇fi,j,f : the time for a spike to go

According to Eq. (11), the spike produced by nf will reachnf 1 soon, while the postsynaptic is still processing B′ . The
summation of these two membrane potentials, like in thephase diagram of Fig. 8(4), will push the differential equation solution to higher membrane potential, facilitating thei,jspike generation thus reducing the processing time 𝜏f .

ppfrom neuron nf to one of its four adjacent neuron
nfp,(p=1,2,3,4), explicated as:
i,j

i,j

𝜇f ,f = (cd)j + ndf ,f = (cd)j + i
p

p

(cd)j
ns

p

i,j

⋅ (2xf ,f − 1)p

(9)

• Total transmission delay dfi,j
:
,f
p

i,j
i,j
i,j
df ,f = 𝜇f ,f + 𝜏f
p
p

(10)

We need four indexes because in our model, each changing of the noise delay (nd)i,j , corresponds to an i-simulation,
which associates to any single neuron. The random variable
0 ≤ xf ,fp ≤ 1, characterizes the delay transmission from nf toi,j
nfp , and a process time 𝜏f . Therefore, each i-simulationi,jdefines df ,f .
p

Once everything is in place, the network operates as
expected: a spike reaches a specific input neuron, say nf .
After a process time tauf , it departs nf and reaches four adjacent neurons nfp,(p=1,2,3,4) with four different axonal transmisi, jsion delays 𝜇f , f defined in Eq. (9).
p

Assume that some additional earlier signals A ′ and B ′ ,respectively, excite the presynaptic neuron f and the postsynaptic neuron fp at the same time (notice that the simultaneous
stimulation of f and fp has been chosen for clarity, the actualinterval between the two stimuli A ′ and B ′ is arbitrary, andwill not influence our general conclusion). We now want tofocus on the case when :
i,j

i,j

𝜇f ,f < 𝜏f
p

p

13

(11)

When this event happens, the total time interval between theinitial presynaptic stimulus caused by A ′ and the final posti,jsynaptic spike df ,f , will be lower than when the inequality1
Eq. (11) is not satisfied. In other words, because of the rani,jdomness factor within the axonal transmission delay 𝜇f ,f ,pthis delay time will be randomly smaller or higher than average. However, only when this is smaller, the arrival of asecond spike in the postsynaptic neuron will induce a loweri,jprocess times in it. Instead in the other case, when 𝜇f ,f isphigher, the process times 𝜏f will not be affected. This asymi,jp

metry creates the effect that we observe in the simulations.
The increasing degree of noise in the network increases theprobability of this condition happening in multiple neurons,
inducing faster spiking activity propagation and improving
the overall circulation of information along with the network, see Fig. 9 for a sketch representing the phenomenon.

Conclusions
We constructed a model network of 400 neurons placedon a 20x20 toroidal grid, neurons on the border are connected to those on the other extreme of the grid realizing
the toroidal structure to emulate real brain patches where
information exchange is confined on localized domains.
We wanted to study the influence of structural heterogeneity in the network and see how it may help the diffusion

Cognitive Computationof information in the system. To do that, we made a single
source of spikes, an initiator neuron nin that is permanently
connected to an external current I that makes it regularly
spiking for the duration of the tests. The spiking activity ofthis neuron represents the information and it is transmitted
to the adjacent nodes on the grid, distributed accordingly
to a von Neumann neighborhood.
Our study shows that in this simplified model network,
the time it takes for the first spike from the initiator neuron nin to reach and induce spiking to distant neurons isdecreasing with the increased heterogeneity in the axonal
propagation delay. The heterogeneity is expressed by anoise factor in the axonal delay term, manipulating thisparameter we found that noise favors a faster spread of
information along with the network.
We devised a theoretical interpretation of this counterintuitive result, which is involving the nonlinear dynamical nature of the timing and propagation of neuronal spike
activity. In our treatment, the presynaptic signal induces aspike with two different time delays, the first is considered
a processing time that is intrinsically related to the differential equation representing the neuronal electrochemical
dynamics and a transmission time that is inherent to thetransfer along the axon. Our theoretical model shows thatthe random variations of transmission times due to the heterogeneity, favor faster processing reducing the overall spiking activity transfer. We have shown rigorously that thisis caused by an increased probability for a second PSP tointervene when the neuron is still elaborating a previous
one, effectively enhancing its action and accelerating the
spike processing.
Since in biological reality, the axonal propagation delay
is a proxy of the axon length or the inter-neuron distance,
this shows that the randomness by which neurons arelocated in the brain could be interpreted as a strategy to
improve and optimize the transmission of information in
biological systems.
These results contribute to the understanding of the fundamental properties of information transfer in the brain,in particular, it gives a model for an explanation of otherparadoxical phenomena in which noise appears to favor
complex computational processes like for example transcranial random noise stimulation (tRNS) that has shown
to be the most effective in enhancing neural excitability
and improving perception and other cognitive tasks inhuman subjects.

Declarations
Ethical Approval This article does not contain any studies with humanparticipants performed by any of the authors.

Conflicts of Interest Marcello Salustri declares that he has no conflictof interest. Ruggero Micheletto declares that he has no conflict of interest.

References
1. Winfree AT. The geometry of biological time. New York:
Springer. 2001;12:2971–4. ISSN 978-1-4757-3484-3. https://​doi.​
org/​10.​1038/​35065​725.
2. Shannon CE. A mathematical theory of communication. Bell Syst
Tech J. 1948;27:373–423. 623–656.
3. Koch C, Laurent G. Complexity and the nervous system. Science.
1999;284:96–8. https://​doi.​org/​10.​1126/​scien​ce.​284.​5411.​96.
4. Stacey WC, Durand DM. Stochastic resonance improves signal detection in hippocampal ca1 neurons. J Neurophysiol.
2000;83(3). https://​doi.​org/​10.​1152/​jn.​2000.​83.3.​1394.
5. Perez-Nieves N, Leung VCH, Dragotti PL, et al. Neural heterogeneity promotes robust learning. Nat Commun. 2021;12. https://​doi.​org/​10.​1038/​s41467-​021-​26022-3.
6. Moret B, Donato R, Nucci M, et al. Transcranial random noise
stimulation (TRNS): a wide range of frequencies is needed forincreasing cortical excitability. Sci Report. 2019;9. https://​doi.​
org/​10.​1038/​s41598-​019-​51553-7.
7. Faisal AA, Selen LPJ, Wolpert DM. Noise in the nervous system.
Nat Rev Neurosci. 2008;9. https://​doi.​org/​10.​1038/​nrn22​58.
8. von Neumann J. Probabilistic logics and the synthesis of reliable organisms from unreliable components. Automata Studies.
1956;34:45–99.
9. Gardner RJ, Hermansen E, Pachitariu M, Burak Y, Baas N,
Dunn BA, Moser M, Moser EI. Toroidal topology of population
activity in grid cells. Nature. 2022;602. https://​doi.​org/​10.​1038/​
s41586-​021-​04268-7.
10. Schwiening CJ. A brief historical perspective: Hodgkin and Huxley. J
Physiol. 2012;590(11). https://​doi.​org/​10.​1113/​jphys​iol.​2012.​230458.
11. Izhikevich EM. Simple model of spiking neurons. IEEE Trans
Neural Netw. 2003;14(6).
12. Wu C, Li Y, Chai S. Design and simulation of a torus topology for
network on chip. J Syst Eng Electron. 2008;19:694.
13. Tzedakis G, Tzamali E, Marias K, Sakkalis V. Routes to chaos induced
by a discontinuous resetting process in a hybrid spiking neuronmodel. Cancer Informat. 2015;14. https://​doi.​org/​10.​4137/​CIN.​S19343.
14. Suwanda R, Syahputra Z, Zamzami EM. Analysis of euclidean
distance and manhattan distance in the k-means algorithm forvariations number of centroid k. J Phys. 2019;1566:696.
15. Eppler J, Helias M, Muller E, Diesmann M, Gewaltig
M-O. Pynest: a convenient interface to the nest simulator. Front
Neuroinform. 2009;2(12). ISSN 1662-5196. https://​doi.​org/​10.​3389/​neuro.​11.​012.​2008.
16. Gewaltig MO, Morrison A, Plesser HE. Nest by example: An
introduction to the neural simulation tool nest. Computational
Systems Neurobiology. 2012;533–558. https://​doi.​org/​10.​1007/​
978-​94-​007-​3858-4_​18.
17. Nobukawa S, Nishimura H, Yamanishi T. The importanceof neighborhood scheme selection in agent-based tumor
growth modeling. Sci Report 2018;8. https://​doi.​org/​10.​1038/​
s41598-​017-​18783-z.
18. Coli M, Palazzari P, Rughi R. The toroidal neural networks. 2000
IEEE International Symposium on Circuits and Systems (ISCAS).
2000;4:137–140 vol. 4. https://​doi.​org/​10.​1109/​ISCAS.​2000.​858707.
19. Expert P, Lord LD, Kringelbach ML, Petri G. Editorial: Topological neuroscience. Network Neuroscience. 2019;3(3). https://​doi.​
org/​10.​1162/​netn_e_​00096.

13

Cognitive Computation
20. Tozzi A, James FP. Towards a fourth spatial dimension of brain
activity. Cogn Neurodyn. 2016;10(3). https://​doi.​org/​10.​1007/​
s11571-​016-​9379-z.
21. Madadi Asl M, Valizadeh A, Tass P. Dendritic and axonal propagation
delays determine emergent structures of neuronal networks with plastic synapses. Sci Report. 2017;7. https://​doi.​org/​10.​1038/​srep3​9682.
22. Schneider A, Hommel G, Blettner M. Linear regression analysis: part
14 of a series on evaluation of scientific publications. Dtsch Arztebl
Int. 2010;107(44). https://​doi.​org/​10.​3238/​arzte​bl.​2010.​0776.

13

23. Chrysafides SM, Bordes S, Sharma S. Physiology, resting potential. NCBI National Center for Biotechnology. 2021.
24. Aldrich SB. The use of multiple neurotransmitters at synapses.
Synaptic Transmission. 2019;449–480. https://​doi.​org/​10.​1016/​
B978-0-​12-​815320-​8.​00021-1.
Publisher’s Note Springer Nature remains neutral with regard tojurisdictional claims in published maps and institutional affiliations.