Distribution of spiking and bursting in Rulkov’s neuron model

Large-scale brain simulations require the investigation of large networks of realistic neuron models, usually represented by sets of differential equations. Here we report a detailed fine-scale study of the dynamical response over extended parameter ranges of a computationally inexpensive model, the two-dimensional Rulkov map, which reproduces well the spiking and spiking-bursting activity of real biological neurons. In addition, we provide evidence of the existence of nested arithmetic progressions among periodic pulsing and bursting phases of Rulkov’s neuron. We find that specific remarkably complex nested sequences of periodic neural oscillations can be expressed as simple linear combinations of pairs of certain basal periodicities. Moreover, such nested progressions are robust and can be observed abundantly in diverse control parameter planes which are described in detail. We believe such findings to add significantly to the knowledge of Rulkov neuron dynamics and to be potentially helpful in large-scale simulations of the brain and other complex neuron networks.


Introduction
Despite the significant progress in neuroscience during the last years, there are still several unveiled questions and mysteries about the brain functionality. The approximately 8.6×10 10 neurons are the basic building blocks of the complex neural network underlying the human brain [1][2][3][4]. There exist several types of neurons and, since the seminal work of Hodgkin and Huxley (HH) [5], numerous models were proposed to describe how neurons function.
Discrete-time single-neuron models producing realistic response properties of real-life neurons were recently a e-mail: jgallas@pks.mpg.de (corresponding author) used to study spatiotemporal patterns of activity observed in sensory processing, memory formation and other cognitive tasks [1][2][3][4]. The popularity of using maps as neuron models is the fact that maps are computationally very inexpensive while providing realistic representations of biological neurons. Many applications of single-neuron models already exist, but it seems that a systematic study of their response as a function of the control parameters lack. Map-based models are easier to treat and analyze and have been used to describe specific situations, especially when considering coupled elements characterized by chemical [28,29], electrical [30,31] or both aspects [32][33][34][35]. The topology of neural networks is also an essential player in their dynamical behavior, as reported for, e.g., scalefree [32], global [36], mean-field [37], small world [38], and Apollonian [39,40] networks. The specific coupling is an essential issue in what concerns both synchronization of neural networks [22,[41][42][43][44], and other dynamical aspects [45].
The purpose of this paper is to explore the control parameter space of one of the most popular and successful discrete model for neurons, namely a map proposed by Rulkov [21]. This model reproduces well many neural behaviors such as emergent bursting from nonbursting cells [19], spiking-bursting [21] and the origin of chaos [46], the birth of self-sustained subthreshold oscillations [47], patterns of bursting [48], spatiotemporal chaos ordering [49], or synchronization features in two neurons [28,50], in large one-dimensional neural networks [42] and two-dimensional lattices [22], and also in scale-free [35,43,51] and in small-world networks [52,53]. The great interest in Rulkov's model lies in the fact that it can reproduce quite well the spiking and spiking-bursting behaviors observed in real-life neurons [20,21].
Complementing results presently available, we report a detailed and systematic fine-scale characterization of the shape, size, and organization of the stability phases observed in the parameter spaces of Rulkov's neuron model. This is done in two ways, by characterizing stability phases with standard Lyapunov exponents, and by using the periodicities computed in discrete time which constitutes a similar technique to that of the socalled isospike diagrams [54][55][56][57][58][59][60][61][62], namely by constructing phase diagrams based on the number of spikes of periodic oscillations associated with periodic phases. As described below, we find Rulkov's neuron model to display stability phases organized in an exceedingly complex way, where periodic regions are abundantly observed in distinct control parameter planes. Despite the overall complexity, it is possible to recognize extended stability phases exhibiting an ordered organization that obeys simple rules.
As it is known, Lyapunov stability diagrams merely divide the control parameter space into just two phases: chaotic and non-chaotic. In sharp contrast, periodicity and isospike diagrams discriminate the nature of all individual phase, thereby providing a really detailed control parameter plane characterization of the oscillations supported by the equations of motion. As illustrated below, isospike diagrams are compelling to characterize stability phases and reveal exciting features that are not visible with Lyapunov exponents. For detailed surveys comparing the Lyapunov and isospike characterization of stability phases see, e.g. Refs. [63,64] and references therein.
Section 2 introduces Rulkov's model briefly. Section 3 presents high-resolution stability diagrams for different slices of the parameter space of the model, focusing in the stability regions which exhibit different patterns of periodicity, i. e., nested families of arithmetic progressions of oscillatory phases. Finally, Sect. 4 presents our conclusions and some perspectives opened up by our findings.

Model and methods
The two-dimensional map of Rulkov is governed by the following equations [21]: where In the above expressions, x and y are the fast and slow dynamical variables, respectively. The map is controlled by three parameters, namely σ, α and μ.
Several earlier works reported parameter planes for the Rulkov model considering μ = 0.001, depicting bifurcation diagrams for the plane (σ, α). Such diagrams record the boundaries between regions of silence, spiking and bursting, and spiking regimes [21]. From such diagrams it is also possible to identify some chaotic points [46], or to delimit a region in which irregular (chaotic) spiking and spiking-bursting are possible to determine by the lines σ th = 2 − α/(1 − μ) and L ts , obtained numerically [42]. Another issue is to determinate well-defined chaotic phases [22], as shown here in Fig. 1a, or to demarcate the some of the above regions using the concept of Neimar-Sacker bifurcation [24].
The aforementioned aspects are shown in Fig. 1 for μ = 0.001. Thus, the delimitation of the regions of regular and chaotic behavior is made with Lyapunov exponents, in Fig. 1b, or with the number of spikes per burst, in Fig. 1c. In this paper, we use the computation of periodicities that also permits us to describe the dynamical behavior of the Rulkov map with the advantage to distinguish the dynamical behavior of both variables x and y. In Fig. 1d, the parameter plane is described in terms of such periodicities corresponding to the fast variable x. This type of diagram was previously used in Refs. [54][55][56]65,66].
We start by focusing on the spiking behavior, either as a burst of spikes or as tonic spiking. To study bursts of spikes, we consider a point in a region of the α × σ plane where such bursts are typical. To this end we fix α = 18.7 and σ = 0.1, and compute the periodicities, say p x and p y , in both fast and slow variables, as well as the number n s of spikes per burst. Using time series, Fig. 2a-d illustrate typical bursting behaviors recorded in the fast variable x, and the periodic features of both variables of the neuron when varying μ. We also show the corresponding phase planes, Fig. 2e-h, for representative values of μ. It is clear that the periodicities p x and p y have the same value, which plummets with the increase of μ. Similar behavior is observed in the number n s of spikes per burst. From Fig. 2 is easy to compute the periodicity of x and y by subtracting two consecutive numbers corresponding to minima for the burst of spikes (x-variable), or in the maxima of the slow variable y. It is noticeable that the periodicities for both variables coincide as they should. In Fig. 2a Fig. 3a-b for such dependences). Consequently, the relationship between n s and p x(y) is linear as depicted in Fig. 3c. It is noteworthy mentioning that the literature, e.g. Refs. [21,47], systematically considers the limit μ → 0, where bursts with many spikes are abundant.
One of the main points of this work is getting the parameter diagrams to study the dynamical behavior of a neuron model. As shown in Fig. 2a-d, the time series of both variables x and y have the same periodicity. Thus, the phase diagrams α vs. σ built using the periodicities of each one of the variables are similar. Henceforth, for simplicity, we focus on the periodicities related to the slow variable y. Indeed, as it is explained in Sect. 3, the phase diagrams of Fig. 4 constitute the heart of the paper because the main results are extracted for the diagrams shown therein. The details of the coloring are explained afterward.
We report a new bifurcation scenario, denoted as a nested period-increasing scenario where each period can be obtained from a linear combination of a basis of two periods indicated by p 1,0 and p 0,1 . The abovementioned basal periodicities expressions constitute a shorthand notation explained in detail in Sect. 3.

Results and discussion
As mentioned, for regular behaviors the periodicities of x and y coincide. Thus, one may use indistinctly x or y to study neuron periodicities. In Fig. 1 we presented phase diagrams for μ = 0.001, the original value considered by Rulkov and in subsequent literature. Here, we consider the dynamics on the α × σ control parameter plane, focusing in the periodic phases observed for μ = 0.25, a value for which the neuron displays particularly a plethora of intricate dynamical behaviors. Figure 4 shows stability diagrams obtained by recording over extended regions of the control parameter plane σ ×α. Panel (a) shows a stability diagram recording p y , the periods observed in the slow variable y, while panel (b) shows the corresponding representation using standard Lyapunov exponents. Manifestly, the periodicity diagram displays much more information than the Lyapunov diagram, which is only able to discriminate periodicity from lack thereof. Figure 4a contains sub-windows with rich patterns, for instance the four regions inside boxes A, B, C, and D, which are shown magnified in Fig. 4c-f.
In order to explain the meaning of our representation describing the periodicities found in the parameter diagrams, we use henceforth the color code located in the lower part of Fig. 4, where each color represents the maximum prime number of the periodicity factored. For instance, the color representing p = 35 = 5 · 7 is 7, but this color might also represent periodicities of 14 This convention guarantees that an individual color will represent the multiples of each prime number from 2 to 127. Chaos is represented in blue in this color code and stationarity (fixed point) in white.
Using the above convention, we computed several phase diagrams, not shown here, for both periodicities p x and p y and μ = 0.25. We found that these phase diagram pairs are very similar in almost all regions, with the same periodicities. Consequently, we arbitrarily selected the slow variable y to study the period-icities in parameter regions where nested periodicities sequences are abundant, as in Fig. 4a.
In Fig. 4a there exist interesting regions that are framed and zoomed-in Fig. 4c-f showing in greater detail representative sequences of nested periodicities. A projection of a line of Fig. 4c for a fixed value of the parameter α allows us to perform a detailed analysis of the distribution of periodicities and detect the importance of what we called primary periodicities. The period-adding behavior is illustrated by using an αinterval of Fig. 4d when σ = −1.2. Another interesting aspect related to the existence of accumulation points is observed in Fig. 4e, f; these points might be determined when σ and α are suitably varied simultaneously, as explained later on. However, the most relevant aspect of all these regions is that they permit describing the sequences of periodicities constituting arithmetic progressions following a single rule explained afterward. In what follows, we analyze each of the phase diagrams that give essential information regarding the dynamical behavior of the neuron. For example, for α = 25.0375 indicated by a green line in Fig. 4c, we obtain the bifurcation features shown in Fig. 5: (a) a bifurcation diagram exhibiting the main periodicities where we observe that the sequences of period cascades are pretty different from the typical period-doubling cascades, (b) a bifurcation diagram in terms of Lyapunov exponents where one sees the regular (negative exponents) and chaotic (positive exponents) behaviors, and (c) a diagram of periodicities p x that illustrates the details of a rich dynamical behavior where it is possible to identify some regularities on the period sequences.
Since there are many possible sequences of arithmetic progressions, we start by considering the region labeled 1 in Fig. 5c and amplified in Fig. 6a. This region exhibits many sequences of arithmetic progressions, seven of them are represented with different colors in Fig. 6b, employing the multipliers k and l accompanying to the basic periodicities p 1,0 = 15 and p 0,1 = 16, respectively. When terms of individual sequences overlap, they can be distinguished by the color related to the Panel c contains five σ-intervals with U -shaped patterns that will serve to study the periodic sequences in more detail Table 1 Main terms and recurrence equations for the sequences indicated in Fig. 6b Sequence classification of the sequence to which the terms belong.
To avoid the apparent messiness in Fig. 6b, we represent each arithmetic progressions separately in Fig. 6ci, adopting the same colors assigned in Fig. 6b. The basis of our analysis of arithmetic progressions is the fact that every sequence consists of specific linear combinations of two elementary periodicities, p 1,0 and p 0,1 , the smallest ones found in a particular region exhibiting the U -shape patterns. Consequently, the n-th term of every arithmetic progression might be written as where k = k(n) and l = l(n). This simple equation may be used to summarize all arithmetic progressions seen in Fig. 6c-i. The corresponding expressions are collected in Table 1. One last property is the fact that in Eq. (4), and for every σ-interval considered, we always find p 0,1 = p 1,0 + 1. Consequently, we have: p n (k(n), l(n)) = (k + l)p 1,0 + l .
Thus, there is only a single fundamental periodicity, namely p 1,0 that can be denoted simply by p. This situation is manifested strongly in the case considered in Figs. 8 and 9. The result expressed by Eq. (4) might be visualized considering a small interval of α along the line σ = −1.2, as indicated on the rightmost region of Fig. 4d. In this α-interval, three periodicities are present: p x1 = 14 (red), p x2 = 27 (brown), and p x3 = 13 (cyan). From Fig. 4d, note that there is an obvious relation-ship p x2 = p x1 + p x3 among the period-addings. To illustrating this relation, in the top row of Fig. 7 we present return maps for some representative values of α, namely, 19.57, 19.77, and 19.97 when σ = −1.2. For each one of these values, the periodicities correspond to 14, 27, and 13 as stated above. We note the shape of the return map with periodicity 27 is a sort of "sum" of the return maps with periodicities 14 and 13. A similar situation occurs for the time-series shown in the bottom row of Fig. 7.
It is interesting to see the ordering sequences of the structures that have p 1,0 = 11 as their main periodicity manifested in the right part of Fig. 8a with an arithmetic progression given by p n = np 1,0 . As illustrated in Fig. 8b, there is also a "secondary" sequence of periodic regions described by the progression p n = (2n + 1)p 1,0 . On the other hand, in the left part of Fig. 8a, the sequence of structures follows an arithmetic progression given by p n = np 1,0 + p 0,1 , where p 0,1 = 12. A remarkable feature observed is the existence of a myriad of accumulation points whose trends are enhanced when viewed in Fig. 9, a phase diagram constructed with Lyapunov exponents. In this figure, white regions correspond to periodic phases with main periodicity p 1,0 = 11, while colors represent chaotic phases. The arrows indicate the trends towards some of the accumulation points in the diagram.
Finally, Fig. 4f shows several sequences of periodicity regions which follow arithmetic progressions and are dominated by a primary periodicity p 0,1 = 13 such as those located in the central part of the phase diagram: among others, where p 1,0 = p = 12. As before, we find the existence of several accumulation points when suitably varying both control parameters simultaneously.

Conclusions and outlook
This work reports a detailed fine-scale investigation of the dynamical response over extended parameter ranges of a Rulkov neuron, a computationally inexpensive model which is able to reproduce the spiking and spiking-bursting activity of real biological neurons. We keep the neuron parameters σ and α fixed and study the variation of μ over wide range of values. Highresolution stability diagrams are reported for both the standard μ = 0.001 value, and for μ = 0.25, a rich fountain of complex neuron oscillations. For small values of μ, we find oscillations with a considerable number of spikes per burst, as well as a dramatic fall of the number of spikes, following a power-law, these quantities when μ increases. We find that μ values different from μ = 0.001, the standard value frequently considered in the literature, allow the system to exhibit exuberant dynamical behavior. We report an in-depth analysis for μ = 0.25. In particular, we report the discovery of sequences of remarkably nested arithmetic progressions among the periodic pulsing and bursting phases of the neuron which, can nevertheless be expressed as simple linear combinations of pairs of certain basal periodicities. Nested progressions are robust and can be observed abundantly.
The dynamical landscape seems to be very complicated and with intricate periodic regions, at least in a first view of the periodicities distribution. Nevertheless, an analysis starting from the one-parameter bifurcation diagram based on periodicities exhibiting certain U -shaped distributions which allowed us to find elementary expressions of the arithmetic progressions describing the periodic regions. We found a general equation describing such arithmetic progressions based on two primary periodicities to express the recurrence equation linked to the sequences. The latter is related to the period-increasing scenario, a feature that seems to be usual in bifurcation diagrams of dynamical systems. Remarkably, the fact that basic periodicities are found to be consecutive numbers allows one to express progressions with just one primary, or basal, period. We also found certain "primary periodicities" which in some cases determine accumulation points of a sequence. Some progressions are found to depend solely on a prime number, with the fundamental periodicity related to such prime. In any case, other sequences appear when the other primary periodicity manifests itself.
From a mathematical viewpoint, the Rulkov map shows a rich dynamical behavior that deserves further study, both analytical and numerical. For instance, in the control plane μ×σ, periodicity regions have a multi slices-shape which seem to convergence towards fixed points. Our findings might contribute significantly to the better understanding of Rulkov neuron dynamics and be potentially valuable for large-scale simulations of the brain and other complex neuron networks.
Acknowledgements This work was initiated during a visit of the authors to the Max-Planck Institute for the Physics of Complex Systems, Dresden, gratefully supported by an Advanced Study Group on Forecasting with Lyapunov vectors. GMRA was supported by the German Academic Exchange Service, Re-invitation Programme. IMJ was supported by the Hungarian National Research, Development and Innovation Office under Grant K-125171. JACG was supported by CNPq, Brazil, Grant PQ-305305/2020-4.
Funding Open Access funding enabled and organized by the Max-Planck Society, Projekt DEAL.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: All data are included in the paper.].
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.