Bifurcations in a new two-cell spiking map: a numerical and experimental study

In this paper, a new nonlinear discrete-time map is presented. The map is based on a second-order dynamics that, despite the limited number of parameters, is able to produce a rich dynamical behavior, including the onset of spiking trends. This latter case will be particularly emphasized, since it allows to consider the introduced system as a novel discrete-time model for spiking neurons. The study is performed by using a numerical bifurcation approach. Moreover, the possibility to obtain a spiking behavior using noise is also shown. The implementation of the map using advanced microcontroller units and the obtained experimental results are discussed.


Introduction
The concept of discrete-time map has a long history in chaos studies [1].The logistic map [2] represents a paradigmatic example of this class of systems as it is able to generate a wide range of dynamical behavior, including bifurcations and chaos, on the basis of a single parameter.Recently, the interest in this area has been renewed [3,4] with the definition of novel discrete maps for engineering applications, such as image encryption [5,6] and image processing [7], or by including the discrete-time counterparts of innovative devices, such as memristors [8].Methods for constructing discrete-time memristive chaotic maps have been recently studied in [9,10].In the latter case, the pinched hysteresis characteristic of a discrete-time memristor induces, in a rather mathematically simple map, various complex behaviors, including a transition from chaos to hyperchaos.The area of research on this topic is appealing and new efforts are also done considering maps that are inspired by quantum game theory [11].
Within the class of nonlinear discrete maps, a particular relevance are assuming spiking maps [12] since they are reliable models of digital neurons.With respect to neuron models in the continuous-time domain [13,14], discrete-time neuron models [15] are, in fact, prone to digital implementations and practical applications.This topic is, therefore, of interest for several areas of research and application, such as neurocomputing and bioengineering.Recently, the efforts in realizing discrete-time spiking neurons have been reviewed in [16], where spiking maps derived by the neuron models, with emphasis on the Hodgkin-Huxley model, are presented.Indeed, simple models possessing a spiking behavior consistent with biological neuron observations have been introduced in [17], where an integrateand-fire mechanism has been discovered in a simple map.The interest in this area of research has been also devoted to discover discrete-time models for particular neurons, like the olivo-cerebellar system [18].Spiking maps are further useful both in modeling [19] and in studying networks behavior [20].An interesting example of discrete-time neuron model is discussed in [19] where it is shown that a discrete-time memristor model provides the key ingredient for modeling magnetic effects on the bursting and firing behavior of biological neurons [21,22].From the point of view of studying nonlinear dynamics in discrete-time models of spiking neurons, the work [23] emphasizes bifurcations, chaos and the richness of the dynamical behavior emerging from spiking maps and in [24] the possibility to have a transition from quasi-period spiking to an hyperchaotic bursting is explored.
In this communication, we present a new secondorder discrete-time map that presents a rich dynamical behavior.Its genesis is related to the two-cell continuous-time system studied in [25].Moreover, the introduction of two novel parameters allows to sensibly increase the complexity embedded in its behavior and to propose a simple hardware implementation of its dynamics.
The paper is organized as follows.In Sect.2, some preliminary details related to the two-cell continuoustime system are reported.In Sect.3, the new discretetime nonlinear map is introduced, discussing the properties of its equilibrium points, while Sect. 4 proposes an intense numerical bifurcation analysis in order to explain into details the emergence of the rich dynamics, also in terms of firing behavior.In Sect.5, the spiking dynamics of the proposed system is discussed providing a strategy based on the injection of noise to trigger the spiking behavior in the stability regions of the parameter space.Moreover, Sect.6 is dedicated to the experimental results derived from the simple implementation of the map obtained by using a lowcost microcontroller and Sect.7 draws the concluding remarks.

The two-cell nonlinear system
Let us consider the following set of nonlinear differential equations where It represents a second-order continuous-time nonlinear system originally introduced in [26] with the aim of extending the paradigm of reaction-diffusion cellular nonlinear networks (RD-CNN).It originates from the so-called Nossek cell [27] with the inclusion of the two external currents i 1 and i 2 whose effect impacts on the existence and stability properties of the equilibria [25].The system in Eq. ( 1) finds important practical applications in control as it is used as the core of central pattern generators to control the gait of hexapod robots [28].
The conditions under which the system in Eq. ( 1) oscillates have been detailed in [25].In particular, the concurrent presence of virtual and real equilibria, imposed by the external currents, ensures a limit cycle oscillation for 0 < μ < s and |i i | < 1.Without loss of generality, in the following we will fix the parameters of the two-cell system according to [25,28] as μ = 0.7, s = 1, i 1 = −i 2 = −0.3.For this set of parameters, the peculiar slow-fast dynamics appears, as shown in Fig. 1, as a consequence of the slow evolution toward the virtual equilibrium, followed by a fast evolution away from the virtual equilibrium, thus generating an homoclinic orbit during which the dynamic flow undergoes sudden variations in the rate of the state-variables change.
Let us now consider, instead of the piece-wise linear function (2), a continuous and differentiable nonlinearity such as This leads to the introduction of a further parameter α which provides the system with an interesting feature.The slow-fast dynamics in fact, yet being preserved for α ∈ [0.65, 1.66], can be modulated by tuning α, as shown in Fig. 2.This introduces an integrate-and-fire like behavior, where one variable assumes the role of action potential and the other assumes the role of a recovery variable.Thus, the system in Eq. ( 1) with nonlinearity as in (3) can be considered a novel model for spiking behavior.The analogy with the Morris-Lecar neuron model should be noticed [29], where the same type of nonlinearity is considered.

The two-cell spiking map
Let us now consider the Euler discretization with t = kT of the dynamics in Eq. (1) with nonlinearity (3) as: The discrete-time system described in (4) maintains the periodic spiking behavior observed in (1) for small values of T , as shown in Fig. 3 where the map iterations and the trajectory in the x 1 − x 2 plane are reported for T = 0.1 and different values of α.Moreover, Eq. ( 4) can be considered a nonlinear discrete-time map with T acting as a further parameter.Varying T , in fact, the behavior of system (4) with nonlinearity (5) drifts from that of system (1) with nonlinearity (3), showing the birth of bifurcation scenarios leading to a complex, yet spiking and bursting, behavior.
In the following, we will keep constant the parameters of the original continuous-time two-cell system, i.e., μ = 0.7, s = 1, i 1 = −0.3 and i 2 = 0.3, and investigate the dynamical behavior emerging in the parameter space T − α.
Let us start from determining the location of the equilibria of the discrete-time map, by imposing the conditions This leads to the following system of nonlinear equations and thus the equilibria must satisfy the equations The first consideration is that the equilibria of system (4) are independent of T .Being the nonlinearities in (7) transcendental functions, the solutions can be derived more efficiently by using a graphical approach based on plotting the nullclines (7) and detecting the intersection points.The nullclines for α = 1, α = 1.666 and α = 1.8 are reported in Fig. 4 showing a transition from one to three and then to five distinct equilibria.The location of the equilibria with respect to α retrieved by adopting this method is shown in Fig. 5.
As concerns the stability analysis, it appears evident that this depends also on the parameter T .In fact, the Jacobian matrix of (4) reads as By inspecting the modulus of the eigenvalues of J evaluated on each equilibria, it is possible to determine their stability as a function of T and α.According to this criterion, in Fig. 6 the regions in the parameter space T − α in which the system (4) admits one, two or no stable equilibria are reported.

Numerical bifurcation analysis
The two-cell discrete-time map introduced in this paper displays a plethora of complex dynamical behavior when varying the two parameters α and T .In this section, we aim at investigating numerically the bifurcation routes to chaos and the peculiar behavior of the system.
The numerical investigation on the system behavior is performed by inspecting the bifurcation diagrams with respect to the two parameters and reporting the corresponding maximal Lyapunov exponent λ max .The Lyapunov spectrum has been calculated by considering the Jacobian matrix in (8) and adopting the algorithm described in [30] to avoid ill conditioning and consequent numerical problems [31].
Let us investigate at first the bifurcation scenario with respect to T .According to Fig. 6, we can focus on three values of α for which the transition between different stability regimes can be observed.For α = 0.5, the system moves from one stable equilibrium to no stable equilibria when increasing T .When α = 1.2, the system admits no stable equilibria for all T .Finally, for α = 1.8 the system undergoes a transition between two stable equilibria and no stable equilibria.
The bifurcation diagram and the maximal Lyapunov exponent with respect to T when α = 0.5 is reported in Fig. 7a.They show the birth of quasi-periodic toroidal oscillations [32], characterized by λ max = 0, occurring when the stable equilibrium disappears.Eventually, the tori collapse onto periodic oscillations, also displaying odd periodicity.
For α = 1.2, chaotic and periodic windows appear along the considered range of T , showing the complex bifurcations cascade of periodic oscillations with increasing periodicity leading to the onset of a chaotic regime, as reported in Fig. 7b.
The bifurcation route for α = 1.8 obtained varying T reflects in the bifurcation diagram reported in Fig. 7c.A period-doubling route sets on giving birth to a sustained chaotic behavior with high values of λ max .
The behavior of the system with respect to Fig. 6 can be further explored by fixing T and varying α.In this case, there are two possible scenarios to be investigated, since the variation of α leads to the transition from one to none and then to two stable equilibria for T < 2, while for T ≥ 2 the system admits one or none stable equilibria when α varies.The two complex bifurcation routes corresponding to these scenarios are reported in Fig. 8.The diagram for T = 1.6 (Fig. 8a) shows the transition between the single stable equilibrium to a toroidal regime, with λ max = 0, which eventually collapse onto a series of periodic oscillations with even and odd periodicity, followed by the transition to a wide window of chaos.The presence of two stable equilibria for α > 1.65 leads to a multistable behavior in which, along with the two stable equilibria, a period-4 cycle can be retrieved, depending on the initial conditions.When T = 2.3 (see Fig. 8b), the toroidal oscillations are soon replaced by the birth of a first chaotic attractor, that coexists with a period-4 oscillation, as reported, for α = 0.56, in Fig. 9.Further increase in α leads to a novel periodic window originating a chaotic behavior, followed by an odd periodic window creating a route to a further chaotic window.Therefore, two different routes to chaos appear with respect to same parameter α.It is interesting to note that the behavior of the map occurring for T = 2.3 and α = 1.7 is intermittent, i.e., alternates between two different chaotic behaviors, as reported in Fig. 10a.Intermittency is robust with respect to the initial conditions, as shown in Fig. 10b.The presence of intermittency, as it occurs in other discrete-time maps [3], is linked to the odd periodic cycles observed in the considered area of the parameter space.The system in Eq. ( 4) with nonlinearity ( 5) is able to generate a wide plethora of non-trivial dynamical behavior.Interestingly, the dynamics of the system make it suitable for modeling spiking neurons.Its behavior, in fact, can be opportunely tuned by acting on system parameters α and T to show different firing patterns.In Fig. 11, we report the transition from a quasi-periodic spiking to a chaotic bursting observable when T = 2.3 and increasing α.The quasiperiodic torus obtained for α = 0.45, in fact, displays a quasi-periodic spiking pattern that becomes periodic for α = 0.5.Chaotic spiking occurs for α = 1.4,whose firing behavior is modulated by further increasing α until the emergence of an intermittent chaotic spiking behavior for α = 1.7.Increasing α has now the effect of passing from a spiking firing pattern to a bursting firing pattern, either periodic or chaotic.More specifically, a periodic bursting is obtained when α = 2.15, while a chaotic bursting appears for α = 2.2.
In order to provide an insight on the system behavior when varying both parameters T and α, the maximal Lyapunov exponent has been evaluated in the parameter space T − α.The diagram reported in Fig. 12 confirms the occurrence of several regions of chaos.

Noise-induced spiking behavior
The spiking dynamics of neuron models is often characterized by the so-called interspike interval (ISI), defined as the time-span between two successive spikes.Since both the continuous-time system (1) with nonlinearity (3) and its discrete-time counterpart (4) with nonlinearity ( 5) have a spiking behavior depending on α, the ISI has been evaluated for the latter for different values of α, as reported in Fig. 13 where the direct effect of the introduction of this parameter on the modulation of the spiking dynamics is shown.The spiking behavior disappears for α > 1.66.Since the discrete-time system (4) is the Euler discretization of the continuous-time dynamics (1), even decreasing T the behavior of the map reaches an equilibrium when α > 1.66, as confirmed in Fig. 14a for α = 1.7.
In order to guarantee a spiking behavior in the system (4) with nonlinearity (5) in a wider range of the parameter α, we can introduce a random perturbation in Eq. ( 4) as: Fig. 13 Interspike interval, in samples, for the spiking dynamics of the system in Eq. ( 4) with nonlinearity ( 5) with respect to α, when T = 0.1 where ξ 1 and ξ 2 are random processes with uniform distribution in the range [−1, 1] and η is the noise level.The effect is to perturb the stability of the coexisting equilibria, thus triggering the slow-fast dynamics responsible for the spiking behavior.The effect of the injected noise is shown in Fig. 14b, in which the behavior of the system for T = 0.1 and α = 1.7 is perturbed with a noise level η = 0.5.
In order to determine the effectiveness of the noise perturbation, the ISI has been evaluated for three different values of the noise level, varying α.In Fig. 15, the average value of the ISI as a function of α is reported for η = 0, η = 0.3, and η = 1.The range of α triggering a spiking dynamics is larger, the higher is the level of noise.The cost, however, is also to increase the variability of the ISIs.
This erratic behavior of the interspike interval occurring in the two-cell map when subjected to the random process, led us to investigate the variability of the distribution of the ISIs, which is clearly affected by the presence of the noise, as a function of η.To this aim, the ISIs have been evaluated fixing T = 0.1 and α = 1.7, and varying η, as shown in the diagram reported in Fig. 16a, where the continuous line marks the average ISI for each value of the noise level.The spiking behavior is triggered for η > 0.3, and the distribution of ISIs becomes more regular when η increases.This latter consideration can be further verified by giving a measure to the variability of the ISIs.Hence, we investigated the coherence of ISIs distribution computing the coefficient of variation (C V ) for the ISIs series obtained where σ ISI and μ ISI are, respectively, the standard deviation and the mean value of the retrieved ISIs series.The C V is reported in Fig. 16b as a function of η.The increase of the noise level leads to a non-trivial reduction of the variability of the ISI distribution.This is a consequence of the nature of the two-cell map for α = 1.7 and T = 0.1, i.e., the coexistence of two stable equilibria.The noise, in fact, triggers the spiking behavior thanks to the interplay of the two equilibria an among them with an intrinsically periodic nature.4) with nonlinearity (5) with respect to α, when T = 0.1, and for η = 0 (average in continuous line), η = 0.5 (average in dashed line) and η = 1 (average in dash-dotted line)

Experimental results
In this section, the implementation of an electronic circuit aimed at realizing the two-cell discrete-time system introduced in this paper is presented.Discrete-time systems can be easily implemented in digital computing devices, with the only constraint of the digit precision.However, we want to implement a hybrid digital/analog circuit which produces an output voltage following the dynamics of the two-cell map in (4).
The circuit is based on the microcontroller board Arduino UNO.It is a multipurpose board equipped with an ATMEL 32bit ARM-based processor, interfaced with an analog-digital converter and a serial/USB interface.It is fully programmable via an integrated development environment in C/C++.
The microcontroller is programmed implementing the iterator (4) with nonlinearity (5), and the two variables, coded as 8-bit words, are sent through the digital output of the board.These digital signals are then converted to an analog voltage by means of a standard R − 2R digital-to-analog converter.Arduino pins from 2 to 9 are used to produce the eight digits of the 8-bit word.The R − 2R digital-to-analog converter is realized on a breadboard, following the schematic reported in Fig. 17.Wires connect the output pins of the Arduino board with the inputs of the converter in such a way that the most significant bit (corresponding to pin 2) is connected to the first resistor R, and so on.The complete experimental setup is shown in Fig. 18.The microcontroller is programmed with the code reported in the "Appendix A." The hybrid circuit reveals to be effective in reproducing the behavior of the two-cell map, even in the 8-bit implementation.In Fig. 19a, b, the periodic behavior of the map is reported obtaining the period-5 cycle, for α = 0.5 and T = 2.3, and period-12 cycle, for α = 1.2 and T = 1.4.The intermittency corresponding to T = 2.3, and α = 1.7, and the chaotic behavior produced by the circuit with T = 2.3, and α = 1.8 are reported in Fig. 20a, b, respectively.Comparing the obtained signals with the numerical simulation, a perfect match as regards the dynamic behavior can be observed.Differences in the signal shape can be noticed, due to the characteristics of the microcontroller implementation.The first macroscopic difference is that the values assumed by the iterator, and then converted to an analog voltage, are not instantaneous but are maintained for a finite amount of time.This time consists in the duration of each loop cycle that depends on the complexity of the executed code and on the delay introduced in the main code (see "Appendix A").The delay is introduced in order to ensure that the main clock of the microcontroller completes a cycle and the digital output is correctly updated.The second difference stands in the amplitude of the analog voltage, that is constrained in the range [0V, 5V].This means that when producing the digital output corresponding to the value assumed by the state variable x 1 at each step, it must be adequately scaled (see "Appendix A").These differences, however, do not impact on the capability of the hybrid circuit The circuit allows also to investigate the triggering of the spiking behavior as a consequence of the introduction of a random process.This needs a reformulation of the code used to program the microcontroller, as detailed in the "Appendix B." The spiking behavior obtained thanks to the injection of noise is reported in Fig. 21a when α = 1.7,T = 0.1 and η = 0.5.Moreover, the experimental bifurcation diagram related to the ISI with respect to η is reported in Fig. 21b, showing a perfect match with 21 Behavior of the hybric circuit implementing the twocell map: (a) spiking behavior obtained thanks to the injection of noise when α = 1.7,T = 0.1, and η = 0.5; (b) experimental bifurcation diagram of the ISI with respect to η the numerical results in Fig. 16a, despite the peculiarity of the output analog voltage outlined above.

Conclusions
Spiking models play a crucial role in the definition of future paradigms for intelligent learning and artificial intelligence.In order to ensure a straightforward implementation on smartdevices and microcontrollers, discrete-time models with spiking and bursting behavior are fundamental.The importance of discrete-time nonlinear maps is also related to their capability of aid-ing in modeling identification from experimental data, where sampled data can be directly classified.
In this paper, we proposed a second-order discretetime map originated from a well-established continuoustime model with a slow-fast spiking dynamics.Rather than being a mere discretization of a continuous dynamics, the proposed system springs a non-trivial quantity of complex dynamics and a generally richer dynamical behavior with respect to the continuous-time counterpart.The use of the discretization step T as a further system parameter allows the discrete-time case to easily achieve strange behavior, including chaos, intermittency, cycles with odd and even periodicity, and even multistability in spiking dynamics.A transition from quasi-periodic to periodic spiking, and then to chaotic spiking, intermittency and bursting behavior can be achieved by tuning a single parameter.
Another crucial aspect is that the introduced model presents several regions in the T − α parameter space in which the coexistence of different behavior occurs.This characteristic makes the neuron model capable of reacting differently in the presence of different initial conditions, thus performing a classification.As the original continuous-time model has been adopted to design robust central patter generators for bio-inspired mobile robots [33], the discrete-time spiking neuron model here introduced can be considered in the context of practical applications on a robotic platform.Under this perspective, the multistable behavior observed may be exploited to allow different robot reactions in the presence of different environmental stimuli suitably mapped onto the dynamics initial conditions.
The possibility of modulating the spiking dynamics through a single parameter and the efficacy of noise in triggering a spiking behavior make the proposed system a realistic neuron model that can be adopted in novel spiking neural network architectures.

Fig. 4
Fig. 4 Nullclines for different values of α: a α = 1, the single intersection corresponds to a unique equilibrium; b α = 1.666, three intersection corresponding to three equilibria; c α = 1.8, five intersection corresponding to five equilibria

Fig. 5 Fig. 6
Fig. 5 Coordinates of the equilibria with respect to α

Fig. 10
Fig. 10 Intermittency in the two-cell map for α = 1.7 and T = 2.3: a intermittent behavior between chaotic dynamics (left panel) and particular of the internal attractor (right panel); b intermittency for ten different initial conditions chosen randomly from a Gaussian distribution with variance sigma = 4 and mean value μ = 0

16
Inducing spiking behavior with noise.a Bifurcation diagram with respect to α on the ISI for T = 0.1: three levels of noise are considered η = 0 (in blue), η = 0.5 (in red), and η = 1 (in black).Lines indicate the average value of the ISI.b Bifurcation diagram with respect to η and c coefficient of variation C V with respect to η, when α = 1.7 and T = 0.1.(Color figure online)

Fig. 17 Fig. 18
Fig. 17 Digital-to-analog converter based on R − 2R configuration with output buffer: R = 10 k , TL-084 operational amplifiers powered with a dual voltage supply V s = ±9 V