Spatial evolution of Hindmarsh–Rose neural network with time delays

Spatial relations between neurons in the network with time delays play a crucial role in determining dynamics of the system. During the development of the nervous system, different types of neurons group together to enable specific functions of the network. Right spatial distances, thus right time delays between cells are crucial for an appropriate functioning of the system. To model the process of neural migration, we proposed simple but effective model of network spatial evolution based on Hindmarsh–Rose neurons and Metropolis–Hastings Monte Carlo algorithm. Under the specific assumptions and using appropriate parameters of the neural evolution, the network can converge to the desirable state giving the opportunity of achieving large variety of spectra. We show that there is a specific range of network size in space which allows it to generate assumed output. A network or generally speaking a system with time delays (corresponding to the arrangement in the physical space) of specific output properties has a specific spatial dimension that allows it to function properly.


Introduction
Complex nervous systems, such as human brain consisting of about 10 billion neurons [1], are the systems that most of important questions related to, remain unanswered.Neurons integrate signals and encode information and although activity of a single nerve cell can be well-modeled today, the phenomena occurring in complex biological neural networks still need to be studied to understand the mechanisms behind various brain functions.
In recent decades, many works focused on the synchronization between neurons [1][2][3][4][5] as this mechanism was believed to be vital for several cognitive functions.The concept of small-world networks by Watts and Strogatz was also included [1,6].The recent innovation is the introduction of time delays to network models [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23] which allows to study the influence of the limited speed of information processing on the network dynamics.This speed limitation is indeed present in neuronal communication as the action potential propagates with the speed of tens of meters per second which is a significant aspect if the physical size of nerve tissue taken into consideration.There are various sources of information time delay between neurons, e.g.limited speed of transmitting action potential through the axon, different types of synapses (chemical, electrical), the release of neurotransmitter, the condition of myelin sheath.Time delays are the important parameters influencing the formation and transitions between neuronal firing patterns [12][13][14].It has been shown that time delay drives the collective behaviors of neurons in the network [15].Applying time delay to the network can enhance its synchrony or stop it from synchronizing [12][13][14][18][19][20].It also changes the bifurcation image of neural activity.Appropriate choice of time delays and coupling strengths allows to observe intermittent behavior of neurons [19] or spiral waves propagating in the network [22].
Recent papers on modeling dynamical systems including experimental reference to biologic neural networks also deserve interest.In [24] authors compare their numerical simulations to recordings of hippocampal place cells of rats.Simultaneously, using artificial neural networks for information processing it is possible to predict behavior of dynamic biologic system straightforward, not through modeling the system itself [25].On the other hand, it is known nowadays that during early stages of development of the nervous system the process of neuronal migration occurs.As the positioning of nerve cells constrains local neuronal signal, the migration leads to grouping of different classes of neurons together providing them appropriate spatial relationships and thus ability of the appropriate interaction [26].Recently, it has been shown that mechanism not covered in the field before are also important factors when neuronal behavior considered, e.g.autapse connections and even physical effects as electromagnetic induction and radiation [17].
What we propose in our paper is to combine two main concepts: time delays corresponding to positioning neurons in the physical space and their space-evolution in time based on the Metropolis-Hastings Monte Carlo algorithm [27].In opposition to our predecessors we do not set constant distribution of time delays but allow them to vary as the nodes move on the surface due to the widely known and simple Markov chain Monte Carlo rules.We use Hindmarsh-Rose neuron model with connection topology of a regular ring lattice which is in general independent from the arrangement of neurons in the physical space.We show that taking a specific target function and appropriate simulation parameters it is possible for the system to evolve from one specified state to another.As a specified state we consider the power spectrum of the output from the network.The networks evolving in our study generate power spectra consisting of the chaotic and periodic parts.In particular, the network can evolve from the state of entirely chaotic output signal to the state of entirely periodic signal, generating a large variety of spectra.We also study physical extent of evolving network to answer the question whether there is a specific physical size that allows it to function in the pre-assumed way.
The paper is organized as follows.Neural network model with time delay coupling is introduced in Sec. 2. In Sec. 3 we describe the model of the spatial evolution and the target function.Details of numerical studies and results are given in Sec. 4. In the last Sec.5 we briefly present conclusions and possible extensions to this study.

Neural network model
A single Hindmarsh-Rose neuron model can be represented with three differential equations, which are as following [28]: where  is the membrane potential of the neuron and , recovery variable, is related to the fast ion current (as Na + or K + ).Adaptation current is represented by  (e.g.Ca + ).The equations for  and  control the fast dynamics and the equation for  controls the slow dynamics of the model. is an external stimulus which drives the activity of the single neuron.The model is able to produce a large variety of signals observed in biological nerve cells, as chaotic spiking, periodic spiking and bursting discharges.As many other authors [2,3] for the numerical simulations we use:  = 1.0,  = 3.0,  = 1.0,  = 5.0,  = 0.006,  = 4.0,  = −1.56.Time series of variable  for different driving current regimes and the bifurcation diagram of a single neuron are presented in Fig. 1.Neural network can be represented as a set of  connected nodes, each of them containing the model (1).The coupling between the nodes is realized by the differential components with time delay present that corresponds to the arrangement of the network in the physical space, e.g. on the plane surface.The neural network model is then described by the following set of differential equations [10,28]: where the effect on  of coupling from  is given by: where ,  = 1, 2, … ,  and  is the coupling strength; (  ) × is the symmetrical (  =   ) connectivity matrix, so if a link between  and  exists then   = 1 and if not then   = 0.The neuron  receives the signal from the neuron  after the time of   , i.e.   is the time delay that only depends on the spatial distance between the nodes.The time delay length is given by: where int[•] stands for integer part, ∆ for time step,   is the distance between neurons  and ,  is a scaling factor (in this work we use  = 13) and (  ,   ), (  ,   ) are the coordinates of neurons.Time delays in the network we consider are selected in a way that corresponds to the arrangement of the network on the plane surface, so   denotes the Euclidean distance.

Evolution model
To study the spatial evolution of the network we worked out the following concept.We choose few of the neurons from the network to be output neurons and collect the signals they produce.The reason behind selecting few neurons instead of all of them is the observation that in the vast majority of complex neural networks (both applied information processing and biological models) not all of the nodes contact with external environment [29,30].Then, using Fast Fourier Transform the power spectrum density of the sum of the signals is computed.The power spectrum density function was chosen as it can be understood as a single output from the group of nerve cells of specified type and, simultaneously, it contains the information on different signals present in the output of this group.Due to the network arrangement on the surface, thus due to the various sets of time delays, the neurons can produce different signal patterns, therefore different power spectra can be observed.The examples of these are presented in Fig. 2. Generally speaking, in the model described above we observed that a neuron tends to produce a chaotic signal if most of its connections are short-distance (distance of ~1 unit in our model) or very-long-distance (roughly, distance of two orders of magnitude longer, ~100 units).A periodic signal can be observed if most of the connections are of the length between the two states mentioned above (long-distance, an order of magnitude longer than 1 unit).The relation between the distance from other neurons and the signal produced is shown in Fig. 3.It is worth to be emphasized that neurons that are not chosen as output neurons also have important impact on the output spectrum.Moreover, notice that the neurons in the network we study are in desynchronized state.
It is now possible to arrange a network where output neurons produce signals from different regimes (chaotic, periodic), so as the output of the network we have a specified (target) power spectrum.Then, we start with a new network, where all neurons are located in one, small-spanned group on the surface.To change the locations of the neurons and therefore time delays of their connections we use the Metropolis algorithm.During each iteration of the spatial evolution we randomly select a node and move it slightly on the surface.If the resulting position causes that current power spectrum fits the target spectrum better, we accept the new position; if not, we accept it with probability  = exp(−  ⁄ ), where  is the change of the target function value and  is the temperature.To measure the similarity of power spectra (target function) we use Pearson correlation coefficient.As the correlation coefficient takes values from the range of [−1, +1], we put the target function as () = 1 − (  (),   ), where   () and   stand for current power spectrum and target power spectrum, respectively.The evolution is stopped when stopping criteria are satisfied (e.g.low value of target function ).
In our simulations the search in physical space is enriched by varying the driving current   of the neuron  that is currently moving.Thus during each evolution iteration we randomly choose a spherical vector from the space (, , ) of two physical dimensions and one current dimension.This additional feature is introduced to obtain a larger variety of possible output states of the system.

Numerical studies
One of the main obstacles to study the system is computational load which increases with the size of the network as the coupling   () based on time delays is different for all of the neurons.Although we use highly optimized code, due to this issue we set a network of 10 neurons; as we will see later this number is enough for observing interesting phenomena.To check the robustness to the network size we also performed the simulation of network of 20 neurons.
The connectivity matrix so the connection topology is a regular ring lattice.Notice, that we do not specify time delays of the connections as they will be clearly defined by the initial positions of nodes, randomly selected on the planar surface later.The next vital preliminary question is the size of node neighborhood in the network.Fig. 4 shows target spectra for different values of neighborhood .For small values of  (e.g. = 2) the network is very sparse, thus, for node  the effect   () from other neurons is very weak and output neurons are driven mostly by their driving currents   .For high values of  (e.g. = 10) the network is dense, so the strong effect   () causes that dynamics of all the output neurons switches to chaotic regime.For further studies we apply the moderate value of  = 4 as it allows both chaotic and periodic output signals to co-exist in one network.As the output nodes we choose 3 of 10 neurons, specifically 1-st, 4-th and 7-th neuron.As  = 4, each -th node is connected to the nodes ( − 2)-th, ( − 1)-th, ( + 1)-th and ( + 2)-th, so the selection is to prevent output neurons to be linked directly what could disturb network output adaptability.
In following numerical studies we use the fourth-order Runge-Kutta algorithm with the time step of ∆ = 0.01.This time step is short enough preserve the neuron signal but also long enough to shorten the computational time of the spatial evolution to acceptable level.Variables   ,   ,   are initially set to random values.The total integration time of each simulation run is 20000.The signals from the last 13000 time units are used to produce the output spectrum of the network (as described in Sec.3.).Then first 1200 samples of spectrum (as it fully covers the region on interest) are smoothed with the moving mean with window of 48 and used to calculate the correlation coefficient.
The other parameters we use in the studies are as follows.The general coupling strength is  = 0.044.The spatial steps of the (, , )-space evolution are depended on  and ∆, ∆ ∈ [−4, +4], ∆ ∈ [−0.02, +0.02].The real step in each iteration of the evolution depends on the random spherical vector as mentioned in Sec. 3. The temperature during the evolution is  = 0.02 and after achieving the criterion  < 0.04 it becomes  = 0.005 to stabilize the system at the final state.Now, having set all parameters for Metropolis time evolution, we put the target spectrum of output neurons as shown in Fig. 4 (chaotic and two periodic peaks).
At the beginning of the simulation the nodes of the network are grouped on the plane surface with random locations (from the range of   ,   ∈ [0,1]) and random driving currents   ∈ [3.8, 4.6] (corresponding to periodic regime for a single neuron).For the beginning state of that kind only the chaotic peak is present in the spectrum, so the value of  is relatively high.Note that the connection topology (ring lattice) is constant during the evolution and the only parameters that evolve are space positions and driving currents of neurons.The exemplary initial state of the network to be evolved is presented in Fig. 5.During the evolution the network spreads on the surface allowing some nodes to move away from the others and therefore changing their time delays.As Metropolis algorithm mainly selects better states,  decreases in time (as shown in Fig. 6).The process often finishes after a finite time with the value of  relatively close to 0, which corresponds to convergence of the network spectrum and the target spectrum (as shown at Fig. 7).The process of evolution is presented in Fig. 8.The spatiotemporal patterns of neuron signals at various stages of evolution are shown in Fig. 9. Action potentials of one of the output neurons changing during the evolution are presented in Fig. 10.Fig. 5. Ring-lattice network connections with neighborhood  = 4 (left; the ring is only for better visualization of links) and the exemplary initial spatial state of neurons on the planar surface (right; the real spatial "look" of the network).The arrangement on the surface clearly defines the time delays of all of the links.Fig. 6.Target spectrum and the starting spectra of two equivalent but different realizations of evolution (the only difference is the initial spatial and current state of the network).The target spectrum is the spectrum that the network output (3 of 10 neurons from the network) fits to during the evolution.The starting spectra are the spectra of network output when the network is in its initial spatial state.All parameters are set as described in the text (Sec.4).

7.
Target spectrum and the spectra of the evolved networks (two separate runs of evolution).The resulting spectra of Run 1 and Run 2 are the output spectra of networks that fitted to the target spectrum, therefore the evolution finished with success.All parameters were set as described in the text (Sec.4).Fig. 8. Target functions () of the networks from Fig. 6 and 7 during the evolution.Horizontal axis is presented in logarithmic scale.The system reaches the spatio-current state of the network that fits the target function but the time needed for the system to evolve is up to 10000 iterations of Metropolis-algorithm evolution.All of the parameters were set as described in the text (Sec.4).Fig. 9. Spatiotemporal patterns generated by all of the neurons in the network at three stages of the evolution: at the beginning (left), after 500 iterations (center) and after 10000 iterations (right).Horizontal axis presents the number of the node, vertical axis is time of the single simulation (the last part).During the whole process of evolution the system stays in the desynchronized state.The output nodes are #1, #4 and #7.At the beginning of the all neurons present bursting chaotic behavior.When the system evolves, neurons #4 and #7 change their dynamics to regular periodic.Fig. 10.Action potentials of one of the output neurons (neuron #4) changing during the evolution.At the beginning of the evolution the neuron presents bursting chaotic behavior.As the evolution goes on, the dynamics of the neuron changes to roughly periodic (500 evolution iterations) and to regular periodic (here, plotted after 10000 evolution iterations) with the period corresponding to one of the target spectrum peaks.
Assuming the same target spectrum but starting from the widely-spanned network (e.g.  ,   ∈ [−10, +10]), we observed that the evolution converges to a positive result (Fig. 11) and resulting evolved network spatial property (mean length of the 3 shortest links) is similar to the one that started from the narrow-spanned group of neurons (as shown in Fig. 12).Taking into consideration Fig. 3, that result implies that the network with time delays of specific pre-assumed output properties has a specific spatial dimensions that allows it to function.To check whether our study is robust to the network size we performed the simulation covering network of 20 neurons with 6 output nodes constructed as in the previous case.One of the resulting power spectra is presented in Fig. 13 and indicates that increasing the network size generates consistent results and that during the spatio-current evolution described here any target spectrum can be achieved (i.e.unlimited number of peaks, unrestricted magnitude of both chaotic peak and periodic peaks) if only sufficient number of output neurons and suitable size of the network provided.

Discussion and conclusions
In this paper, we have investigated the phenomenon of spatial evolution of the Hindmarsh-Rose neural network with time delays depending on the arrangement of neurons on the physical surface.It has been shown that under the specific assumptions and using appropriate parameters the neural evolution can converge to the desirable state, despite the simplicity of the target function and evolution model, and a large variety of spectra can be achieved.This is an interesting phenomenon as it is known that periodic and chaotic stimuli can change behavior of the neurons that are stimulated and the "control" network generating a large variety of signals can be evolved using a simple physical process.We have clearly shown that there is a specific range of network size in space which allows it to generate pre-assumed output.If this size is not reached the output dynamics is always chaotic.In turn, if all of the links exceeds the threshold length then the network loses its ability to produce chaotic signals.If links are very long periodic dynamics of output neurons also degenerates.It can be concluded that a network (or a system, generally speaking) with time delays (corresponding to physical space) of specific pre-assumed output properties has a specific spatial dimension that allows it to function.It is in compliance with earlier results where different phenomena like intermittent behavior or spiral waves were observed due to the specific set of conditions including proper time delays [19,22] but here, in contrast, the network evolved itself fitting its output to the target.The successful use of Metropolis-Hastings Monte Carlo algorithm to evolve the system indicates that this approach can be used to search for specific states in the multidimensional parameter space of neural networks.
There are many possible extensions to our work.The networks we simulated were small due to the huge computational load.Varying parameters as driving current   or coupling strength   can provide further interesting information.It is also possible to apply another target function or more complex dynamics to the model.Ultimately, modeling evolution of larger networks with small world connectivity matrices with shortcuts and clusterization would be an issue of the highest interest.
We hope our study will contribute to the understanding of complex behavior in population of interacting neurons.

Fig. 2 .
Fig. 2. Examples of different power spectra.The wide peak on the left of each diagram (0 ÷ 200 a.u.) is a result of the chaotic output neurons and peaks on the right-hand (600 ÷ 1100 a.u.) are result of the periodic output neurons.

Fig. 3 .
Fig.3.Relation between the distance from other neurons and signal produced.When the output neuron is located close to the other neurons (all links of the length of 1) it produces a chaotic signal.When proper distance applied (all links of the length of e.g.10), then signal turns to periodic one defined by driving current   .If further increase applied to the distance the signal starts to degenerate.

Fig. 4 .
Fig. 4. The spectra for different values of neighborhood .Upper-left:  = 2, upper-right:  = 4, lower-left:  = 6, lower-right:  = 10.For further studies we apply the moderate value of  = 4 as it allows both chaotic and periodic output signals to co-exist in one network.

Fig. 11 .
Fig. 11.Spectra and () diagrams of the networks evolved from the widely-spanned initial spatial conditions (  ,   ∈ [−10, +10] -left diagrams;   ,   ∈ [−50, +50] -right diagrams).It can be noted that the chaotic peak on the right-hand spectrum diagram is not completely fulfilled.If initial conditions are more sparse the evolution is less likely to converge.All of the parameters were set as described in the text (Sec.4).

Fig. 12 .
Fig. 12.Average length of the 3 shortest links in the network, averaged over 5 evolution runs.Both, initially widely-spanned and initially narrow-spanned networks tend to have similar value of this parameter if evolved positively (the resulting spectrum matches the target spectrum).Error bars represent standard deviation.This result indicates that there is a specific minimal physical size of the network to realize the target output spectrum.

Fig. 13 .
Fig. 13.() diagram and the spectrum of the network of 20 neurons evolved in the same way as the previous examples with 10 neurons.The target spectrum consists of one chaotic and two periodic peaks.Each peak is the result of the output signal from 2 output neurons.The success of the evolution of the network with 20 nodes indicates that the network and evolution models are robust to the network size.