Fourier analysis of a delayed Rulkov neuron network

TL;DR

This paper uses Fourier and wavelet analysis to study synchronization in a delayed, chaotic Rulkov neuron network, developing an algorithm to optimize delay for improved synchronization.

Contribution

It introduces a novel Fourier-based algorithm to optimize delay in a neuron network, enhancing synchronization in chaotic, small-world networks.

Findings

Algorithm improves synchronization with optimal delay

Fourier and wavelet transforms characterize neuron behavior

Robustness tested with noisy, non-homogeneous neurons

Abstract

We have analyzed the synchronization of a small-world network of chaotic Rulkov neurons with an electrical coupling that contains a delay. We have developed an algorithm to compute a certain delay whose result is to improve the synchronization of the network when it was slightly synchronized, or to get synchronized when it was desynchronized. Our general approach has been to use tools from signal analysis, such as Fourier and wavelet transforms. With these tools, we have characterized the behavior of the neurons for different parameters in frequency and time-frequency domains. Finally, the robustness of the algorithm has been tested by using non-homogeneous neurons affected with a parametric noise.

Equations

$$x_n=\frac{\alpha }{1+x_{n-1}^2}+y_{n-1},$$

$$y_n=y_{n-1}-\beta \left(x_{n-1}-\sigma \right),$$

$$\sum _j{\delta }_{ij}(x_{j,n-\tau }-x_{i,n-1}),$$

$$x_n=\frac{\alpha }{1+x_{n-1}^2}+y_{n-1}+$$

$$+\delta (A{x}_{n-\tau }-D{x}_{n-1}),$$

$$y_n=y_{n-1}-\beta (x_{n-1}-\sigma ),$$

$$\sum _j{\delta }_{ij}(x_{j,n-\tau }-x_{i,n-1})=$$

$$=\delta \left(\sum _j{a}_{ij}(x_{j,n-\tau }-x_{i.n-1})\right)=$$

$$=\delta \left(\sum _j(a_{ij}{x}_{j,n-\tau })-d_{ii}{x}_{i,n-1}\right).$$

Full text

Rulkov model; synchronization; neuron networks; signal analysis

I Introduction

Neuronal dynamics constitutes an important field of research in mathematical biology and biophysics. This research increases the understanding of mental illnesses, provides tools for their diagnosis, or simply helps to understand fundamental processes in the brain.

Many aspects of the dynamics of the brain can be described by using mathematical models, which basically can be divided in continuous (ODEs or PDEs) or discrete in time (Recurrence equations or maps). Each one of them possesses advantages and disadvantages. Continuous models have a higher computational cost, which is increased even more when we consider them forming a network. Some well-known examples are the Hodgkin–Huxley, the FitzHugh–Nagumo, the Hindmarsh–Rose and the Kuramoto models Hodgkin and Huxley (1952); Nagumo et al. (1962); Hindmarsh and Rose (1984); Kuramoto (1984). In spite of the possible limitations of discrete-time models, they possess a clear advantage from the computational point of view. Since they are basically constituted by recurrence relations, they can be easily computed, so that the computation time is much lower than in the case of continuous models. Furthermore, the mathematical analysis can be simpler. However, the interpretation of the results might not be that simple. Among the discrete-time models (also called map-based neuron models) we can mention the Izhikevich, the non-chaotic Rulkov, and the chaotic Rulkov models, among others Izhikevich (2003); Rulkov (2002, 2001); Ibarz et al. (2011).

In this paper, we want to analyze a more complex structure where neurons are connected between them forming a particular network. Consequently, the mathematical tools that represent a certain pattern of connections are graphs, where the nodes represent the neurons, and the edges the connections. They can be expressed in matrix notation by using an adjacency matrix, where the neurons $i$ and $j$ are connected if the entry $(i,j)$ is 1, while they are not connected if entry $(i,j)$ is 0.

Small-world networks have been used in order to reproduce the brain connectivity He et al. (2007). They are a kind of random graphs where the typical distance between nodes $L$ scales as the logarithm of the number of nodes $N$, i.e., $L\propto\log{N}$. Among other relevant properties we can mention that some nodes possess much more connections than the average node, and their nodes usually form clusters Watts and Strogatz (1998).

We are considering the case in which the interaction between two neurons is modeled mathematically with an electrical coupling Nordenfelt et al. (2013), that is, we are not considering for convenience chemical couplings. Physiologically, this represents the transportation of ions through the gap junction channel between very close neurons. Moreover, the electrical couplings have been implemented with a time delay representing the communication time between neurons.

Finally, we have also assumed that not all neurons are equal, i.e., they do not have the same values of their parameters. We model this non-homogeneity by adding a Gaussian white noise in a key parameter of the chaotic Rulkov neuron.

The synchronization of the neurons in a network is a very interesting emergent property. One of the practical interest of this property is its relationship with some mental diseases. In particular, in the Alzheimer disease or schizophrenia, the synchronization of neurons is lower than the expected. However, in epilepsy or Parkinson disease the synchronization of neurons is greater than expected Uhlhaas and Singer (2006); Rubchinsky et al. (2011).

Our approach is based on signal analysis, i.e., using Fourier and wavelet transform. The point of these tools is that the Fourier transform allows us to study the frequency domain of our data, while wavelet transforms do the same but in a time-frequency domain. With the Fourier transform we can study the frequencies in an asymptotic stable oscillation while by using wavelet transforms, we can study the evolution of these frequencies in time.

The main goal of this paper has been to develop an algorithm that improves the synchronization of neuron network. This algorithm computes a delay that included into the electrical coupling, helps to improve the synchronization of the neuron network. We tested this algorithm for a non-homogeneous neuron network, showing its robustness for low and medium noise intensities.

II The Model

The chaotic Rulkov model Rulkov (2001); Ibarz et al. (2011) is a well-known map-based neuron model defined as

[TABLE]

where $\alpha$, $\beta$ and $\sigma$ are constant parameters, the variable $x_{n}$ is the fast variable representing the transmembrane voltage of a single neuron, and $y_{n}$ is the slow variable standing for the slow gating process. The difference of the time scales between the two variables is determined by a sufficiently small value of the parameter $\beta$, i.e., $0<\beta\ll 1$. For different choices of the parameters, the model described in Eq. (1) can show different behaviors of the neurons in a qualitative manner, such as: rest, spikes, bursts, and chaotic spikes Wang and Cao (2014).

As we explained earlier in the Introduction, the main goal of this paper is to study the global behavior of some coupled neurons, as it happens in the real brain. For this purpose, we need to connect these neurons mathematically using a network, in our case we will use a small-world network.

If we want to build a small-world network, we need to order the $N$ nodes, and then connect each node with its $k$ first neighbors in a cyclical way. For instance, the first node should be connected with the $k$ latest nodes. Finally, each node has a probability $p$ to rewire one of these edges with any other node. Thus, this kind of graph has three parameters: (i) $N$, the number of nodes; (ii) $k$, the number of first neighbors connected with one neuron; and (iii) $p$, the probability of rewiring an edge. This last parameter indicates the random nature of the graph, i.e., as $p$ increases the graph will be more random.

The randomness of the graph implies that, for the same choice of parameters, there might be different degrees of synchronization that will depend on the graph itself. This phenomenon makes sense in the real brain, since all brains have the same macroscopic structure but from a microscopic point of view the connection patterns are different for every person.

The node connections given by the graph represent electrical couplings between the neurons. Mathematically, the coupling is modeled with the addition of a term in our model equation. Furthermore, since the distance between neurons is not zero, we have a non-zero electrical transmission time, what implies that we need to introduce a delay $\tau$ in the coupling. For the neuron $i$, the mathematical representation of these couplings have the form:

[TABLE]

where $i$ is the evaluated neuron, $j$ the remaining neurons, and $\delta_{ij}$ a coupling constant between neurons $i$ and $j$. Since the electrical connection between two neurons is bidirectional, the coupling strength $\delta_{ij}$ is equal to $\delta_{ji}$. This coupling would be equal to zero if and only if $i$ and $j$ are not connected.

Summarizing all this previous information, we proceed to build the complete model of the network. In this case we will also fix the same value for the delay and the coupling strength for all neurons, so that the model can be written, in matrix notation, as

[TABLE]

where $x$ and $y$ are vectors, and the entry $i$ is the value of $x$ or $y$ for the neuron $i$; $A$ is the adjacency matrix of the graph, $\tau$ the delay, and $\delta$ is a coupling constant defined as $\delta=1/(3(k+1))$. The operator $D$ is a diagonal matrix where each entry $d_{ii}=\sum_{j}a_{ij}$, i.e., it is the number of connections of the neuron $i$. The operations $x^{2}$ and $\alpha/(1+x^{2})$ are point-wise for each element of the array. The coupling term is the same as in Eq. (2) but in a matrix notation:

[TABLE]

The first term is the matrix vector product $Ax_{n-\tau}$ and the second term is $Dx_{n-1}$.

Some properties that we can study about the global behavior of the network are periodic; in fact, we will focus on them. These properties can be easily studied using spectral methods like the Fourier or the wavelet transform. So, once the time signals for every neuron have been obtained, we use the FFT (Fast Fourier Transform) and the CWT (Continuous Wavelet Transform) algorithms of MATLAB for their analysis.

The FFT gives us information about the phase and the frequency of the neurons while the CWT performs a time-frequency analysis, where we can see how the frequency varies in function of time. We use the Morse wavelet with symmetry parameter $\gamma=3$ and time-bandwidth product $P^{2}=60$ Olhede and Walden (2002).

III Frequency Analysis for a Single Neuron

In this section we will analyze the behavior of the Rulkov neurons (Eq. (1)) for different values of the parameter $\alpha$. In our case, we will analyze not only the time series but its frequency domain, because it will give us a better understanding of the behavior of the system.

We present some of the more interesting results in Fig. 1. On the left panels we represent the time series corresponding to different values of $\alpha$, while in the right panels we plot their corresponding Fourier transform. We can observe that the time series present the typical spiking behavior for low values of $\alpha$, where their corresponding Fourier transforms have a maximum for their fundamental frequency and its harmonics, as shown in Figs. 1(a-b). As $\alpha$ increases, the previous spiking regime starts to present a decreasing oscillation after having reached its maximum. In this case, we obtain a transform quite similar to the one shown for lower values of $\alpha$. Both cases present the characteristic spectrum of periodic but not smooth functions like the sawtooth function, as shown in Fig. 1(c-d).

If we increase a bit more the value of $\alpha$, we reach values where the behavior changes dramatically, presenting a mixture of spikes interspersed with a variable number of bursts. For this case, the frequency domain shows different maxima corresponding to each kind of bursts-spike group, as shown in Fig. 1e-f.

Finally, when the value of $\alpha$ is higher the system presents a plain bursting dynamics, as can be observed in Fig. 1g. Its corresponding spectrum is centered in the principal firing frequency that can be seen in Fig. 1h. The lower maxima in frequencies represents the internal oscillations inside the bursts. Finally, if we increase a bit more the value of $\alpha$, we have chaotic bursts characterized by a noisy spectrum, where no maxima are observed in the frequency domain, as illustrated in Fig. 1(i-j).

In the cases $\alpha=4$ and $\alpha=4.75$, we have observed that the neurons did not stabilize in a specific periodic oscillation; so we have used the wavelet transform for a better understanding of these neurons in their time-frequency behavior. In the case $\alpha=4$ there might be a periodic pattern of bursts-spike groups. And in the case $\alpha=4.75$ there might be some kind of order in the chaotic oscillations. Our numerical simulations shown in Fig. 2a indicate that there is a kind of noisy alternation of the frequencies implying a change of the number of bursts for each kind of bursts-spike group, as we expected. At the same time, no fundamental frequency appears at any time as shown in Fig. 2b. This kind of behavior is expected from a chaotic oscillation.

IV Neurons in a Small-World Network

When we have a group of interacting neurons, a very interesting global property is their possible synchronization. In the rest of the paper we will focus on the study of this property. We use a use a small-world network with parameters $N=50$, $p=0.2$ and $k=2$ to model the connections of the neurons.

The numerical simulations of the dynamics of this network are plotted in Fig. 3. As we did previously for the case of a single neuron, we will start our analysis with a low value of the parameter $\alpha$, and will study the changes in the global behavior as $\alpha$ increases.

In Figs. 3(a-b), we plot the results for $\alpha=2$. In this case, we observe how all neurons are spiking in a fully synchronized state. This behavior is similar to the one observed in the case with only one neuron (Figs. 1(a-b)). When we increase the value of the parameter $\alpha$ $(\alpha=3.75)$, we observe that the neurons are not synchronized. This is shown in Fig. 3c, where the time series corresponding to different neurons are plotted in different colors. This desynchronization manifests itself in the corresponding frequency domain (Fig. 3d), where a certain dispersion can be observed around the largest maximum.

Finally, if we increase the value of the parameter $\alpha$ to $4.75$, we observe that the system possesses a bursting dynamics, as shown in Fig. 3e. If we represent the corresponding transform, Fig. 3f, we see two maxima where the largest one corresponds to the firing frequency. Now, we can appreciate a very interesting change in the behavior of the network. The bursting dynamics that happens for $\alpha=4.25$ (see Fig. 3f) in the case of the single neuron, appears here for $\alpha=4.75$.

The behavior of the neurons connected in a network depends of the network itself. Thus, their behavior or their synchronization will strongly depend on the corresponding network, even though their global behavior might be similar. The graph used in this section is plotted in Fig. 4, where each color represents a different fundamental frequency for the neuron in a simulation with $\alpha=3.75$. We can see that neurons with the same fundamental frequencies are usually clustered.

V Delay Algorithm

Analyzing the results obtained in the previous section, we observe that for some values of the parameter $\alpha$ not all neurons oscillate at the same frequency. One option for improving synchronization between the neurons is adding a delay to the model, so our goal here is to focus on developing an algorithm to calculate the delay $\tau$ for which this kind of systems can be synchronized showing its performance, as well.

The steps of the algorithm that we have developed to calculate the delay $\tau$ used in Eq. (3) are:

1. (i)

Compute the fundamental frequency (non-zero frequency with the largest amplitude) of each neuron. 2. (ii)

Calculate the period that corresponds to the fundamental frequency of the previous step. 3. (iii)

Use this period as the delay and compute again Eq. (3).

We use the delay $\tau$ obtained with this algorithm to solve the previous network for $\alpha=3.75$, using the same network parameters: $N=50$, $p=0.2$ and $k=2$. The results of this process are represented in Figs. 5, 6b and 7b, in contrast with Figs. 3c and 3d, and Figs.6a and 7a, where no delay was used in the simulation.

In the time series of Fig. 5 we plot the same neurons as represented in Fig. 3c. Now, we can observe that after an initial desynchronized state, the system starts synchronizing as time passes until its final state of synchronization. This synchronization is fast in frequency as we see in Fig. 7 but slow in phase, i.e., 30000 iterations are needed to achieve it. In the frequency domain plot, the maxima are harmonics of the fundamental frequency as in Fig. 1b or Fig. 3b, in the same manner as what happens in a spiking dynamics. Space-time plots are very convenient to compare the synchronized and non-synchronized regimes, because they allow to show the state of all neurons for a certain time interval, helping to visualize whether the system is synchronized or not.

In Fig. 6, we draw the space-time plot of the latest iterations of the simulations. In the right panel where the delay has been used, the synchronization is clearly visible, while in the left panel, where there is no delay, there is no synchronization. In Fig. 7, we have plotted the wavelet transforms of Fig. 3c and Fig. 5a. We can observe in the left panel, where we did not use any delay, that there is no steady frequency at all. However, we observe a steady fundamental frequency in the right panel where there is a delay. These results confirm all that we have already pointed out previously with the time series analysis.

VI Effect of noise on $\alpha$ parameter

In order to simulate a more realistic neuron network, we are going to introduce some non-uniformity in the neurons. As a matter of fact, we are interested to analyze how this non-uniformity affects the global behavior of the network and the robustness of our delay algorithm.

We implement this scheme with the addition of an internal noise in parameter $\alpha$. Specifically, we have used a Gaussian white noise $\xi$ with zero mean, variance one and uncorrelated in space, i.e., $\left<\xi_{i},\xi_{j}\right>=2\delta(i-j)$, which is applied to the parameter $\alpha$ for each neuron. So that now we have $\alpha=\left[\alpha_{1},\alpha_{2},...,\alpha_{n}\right]^{T}$ in Eq. (3), where $\alpha_{i}=\alpha_{0}+D\xi_{i}$ and $D$ is the noise intensity. We fix $\alpha_{0}=3.75$ as in the previous section, and we use the same value of $D$ for all neurons.

We have carried out several numerical simulations for different noise intensities, and without delay. For a low noise intensity, the system is slightly more desynchronized and the neurons present more fundamental frequencies than in the noiseless case (see Figs. 8(a-b)). Now, we proceed to repeat the simulation in presence of noise and with the delay. Our results show that the system gets synchronized due to the delay, as shown in Figs. 8(c-d). To compare those situations we have used their space-time plots, Fig. 9, where we can see the synchronizing effect of our delay.

We are considering now a high noise intensity, where $D=0.75$. We use this value of $D$ because is the lower intensity for which the system cannot be synchronized using our delay in any of our simulations. With a high noise intensity, the system seems less desynchronized than the case with a low noise intensity (compare Fig.8a with Fig.10a, and Fig.9a with Fig.11a). This is counter-intuitive, since we would expect a more disordered state with a higher noise intensity.

Even though the situation now seems more ordered, our algorithm cannot synchronize neurons completely. Yet it improves the synchronization, as it can be seen in Figs. 10(c-d) and 11b. The reason is that the bursts do not need to rest before starting firing again. The fundamental frequency is the firing frequency, but here we have lower frequency maxima representing bursts that do not rest between firing.

In the left panel of the wavelet transform plot, Fig. 12a, there is no fixed frequency over time. Nevertheless, in the right panel, Fig. 12b, we can see a steady fundamental frequency but not in its harmonics. This could be an effect of the algorithm failure.

VII Conclusion

We have characterized the behavior in frequency domain of chaotic Rulkov neurons. This provides useful information on the particular dynamics of the neurons. The same analysis has also been done on a small-world network of neurons.

We have been mainly focused in what happens with the synchronization in phase and frequency of the network by using tools derived from signal analysis. We have been interested in the effect of adding a delay to the network couplings. Our numerical simulations have proved that with the use of a particular delay, the synchronization of the system is improved when it was slightly synchronized or gets synchronized when it was desynchronized. We have developed an algorithm to compute this particular delay, which is related to the fundamental frequency of some neurons.

A more realistic approach to a real neuron network is to assume that not all neurons are equal, i.e., they do not have the same values of their parameters. To model this non-homogeneity with an internal white noise applied to parameter $\alpha$, which is one of the key parameters of the Rulkov neuron. We have tested our algorithm to compute the delay that synchronizes the network in this case. Our algorithm is robust against low and medium noise intensities, but it fails for higher intensities.

Acknowledgments

This work was supported by the Spanish State Research Agency (AEI) and the European Regional Development Fund (FEDER) under Project No. FIS2016-76883-P.