Chaotic Resonance in Typical Routes to Chaos in the Izhikevich Neuron Model

Chaotic resonance (CR), in which a system responds to a weak signal through the effects of chaotic activities, is a known function of chaos in neural systems. The current belief suggests that chaotic states are induced by different routes to chaos in spiking neural systems. However, few studies have compared the efficiency of signal responses in CR across the different chaotic states in spiking neural systems. We focused herein on the Izhikevich neuron model, comparing the characteristics of CR in the chaotic states arising through the period-doubling or tangent bifurcation routes. We found that the signal response in CR had a unimodal maximum with respect to the stability of chaotic orbits in the tested chaotic states. Furthermore, the efficiency of signal responses at the edge of chaos became especially high as a result of synchronization between the input signal and the periodic component in chaotic spiking activity.

A considerable number of studies have been conducted on chaos and bifurcation in spiking neural systems, generating model systems that include the Hodgkin-Huxley, FitzHugh-Nagumo, and Hindmarsh-Rose models 35 . In particular, the Izhikevich neuron model, as a hybrid spiking neuron model, combines a continuous spike generation mechanism and a discontinuous after-spike resetting process; thus, the model can induce many kinds of bifurcations, and reproduce almost all spiking activities observed in actual neural systems simply by tuning a few parameters 36 . In addition, the variety of reproduced spiking patterns is high in comparison with other spiking neuron models 37 .
The hybrid spiking neuron model is one of the piecewise-smooth dynamical systems, in which dynamics are switched according to the state of the system 38 . Saito and colleagues have conducted chaos/bifurcation analysis and circuit implementation against piecewise-constant dynamical systems, and piecewise-linear dynamical systems, as simplified versions of the piecewise-smooth dynamical system [39][40][41] . In particular, Tsubone et al. proposed a systematic method to predict parameter regions for chaotic states using an analytical approach in the piecewise-constant dynamical system 41 . While in general, piecewise-smooth dynamical systems include non-linear terms similar to those seen in the Izhikevich neuron model, an approach for evaluating Lyapunov exponents and characteristic multipliers that considers the saltation matrix 38 through simulations against exhaustive parameter sets is needed. On considering this approach, it is clear that this model has various kinds of bifurcations and routes to chaos when under the effect of the state-dependent jump in the resetting process 34,[42][43][44] . However, the signal responses of CR have not been evaluated in chaotic states produced through different routes.
In our preliminary work, we confirmed the presence of CR in chaotic states induced by different routes (i.e., the periodic-doubling bifurcation route and intermittency route to chaos) in the Izhikevich neuron model 45 . In this paper, we build on our previous work and evaluate the signal responses in CR, and compare the characteristics across these chaotic states through two methods. We first examine the dependence of the signal response on the maximum Lyapunov exponent; then we identify the resonant zone in the parameter space of the applied signal frequency/amplitude.

Materials and Methods
Izhikevich neuron model. The Izhikevich neuron model 36, 37 is a two-dimensional ordinary differential equation of the form and with auxiliary after-spike resetting Here, v and u represent the membrane potential of a neuron and the membrane recovery variable, respectively. We extended Eq. (1) using a weak periodic signal I ext (t) as follows: 2 ext in which we adopted I ext (t) = Asin(2πf 0 t). Note that the sinusoidal signal was utilized merely as a typical example of a signal in a neural system.

Evaluation indices. Indices for evaluation of chaos and bifurcation.
To quantify the chaotic activity in the Izhikevich neuron model, the Lyapunov exponent with a saltation matrix is utilized. On a system with a continuous trajectory between the i-th and the (i + 1)-th spiking times, (t i ≤ t ≤ t i+1 ), the variational equations (1) and (2) are defined as follows: where Φ, J, and E indicate the state transition matrix, the Jacobian matrix, and a unit matrix, respectively. At t = t i , the saltation matrix is given by i In the above, (v − , u − ) and (v + , u + ) represent the values of (v, u) before and after spiking, respectively. In case spikes arise in the range [T k :T k+1 ] [ms], Φ k (T k+1 , T k ) (k = 0, 1, …, N − 1) 43 can be expressed as In our simulation, we set T k+1 − T k as the time required for 20 spikes (i = 20). We set 1000 [ms] as the maximum value in case T k+1 − T k lasts for 1000 [ms] before 20 spikes occur.
In order to conduct bifurcation analysis in the system with a state-dependent jump, we set a Poincaré section Ψ(v = 30). The dynamics of system behavior on Ψ are indicated by the Poincaré map u i+1 = ψ(u i ) where u i is the value of u on Ψ. In the literature 42 , the stability of a fixed point Here, u 0 = (v 0 , u 0 ) indicates the initial value of orbit u = (v, u) at t = t 0 . |μ l < 1|, μ l = −1, and μ l = 1 represent the stable condition, period-doubling bifurcation, and tangent bifurcation, respectively.
Indices for evaluation of signal response. We calculated the timing of the spikes against signal I ext (t) by using a cycle histogram F t ( ) 33 For example, for T 0 = 10, if the spike times were t k = 2, 6, 12, 16, 26, the values of t k mod (T 0 ) were 2, 6, 2, 6, 6. The cycle histogram then became F(2) = 2 and F(6) = 3.
To quantify the signal response, we used the following index of Eqs (11) and (14). The mutual correlation C(τ) between the cycle histogram F t ( ) of the neuron spikes and the signal Ĩ t ( ) ext is given by Scientific RepoRts | 7: 1331 | DOI:10.1038/s41598-017-01511-y FF 2 For the time delay factor τ, we checked max τ C(τ), i.e., the largest C(τ) between 0 ≤ τ ≤ T 0 .

Results
Parameter region for evaluating signal responses. Initially, we determined the parameter regions where the chaotic state is produced. The left panels of Fig. 1(a) and (b) show the dependences of the maximum Lyapunov exponent λ 1 on parameters c and d in the region around parameter sets for the spiking patterns of RS, IB, and CH (see the right part of Fig. 1(a)) and the region around the parameter set proposed by Izhikevich for chaotic spiking (see the right part of Fig. 1(b)), respectively. The chaotic states (λ 1 > 0) exist in −59 ≲ c ≲ −40, d ≈ 1.0 in the former case, and d ≲ −13 in the latter case. As the parameter regions for evaluating CR, we chose 0.82 ≤ d ≤ 0.92 in the former region (called region #1 below), and −15.5 ≤ d ≤ −11 in the latter region (called region #2 below). Figure 2 depicts the bifurcation diagram of u i on Poincaré section Ψ (black dot) and Lyapunov exponents (red dotted (j = 1) and green dashed (j = 2) lines) as a function of parameter d in region #1 case (a) and region #2 case (b). In Fig. 2(a), the period-doubling bifurcation (μ l = −1) arises at d ≈ 0.8348, 0.8828, 0.8916, 0.894, and the chaotic state (λ 1 > 0, λ 2 = 0) appears d ≳ 0.894. Hence, the period-doubling bifurcation route to  chaos exists in this region. However, as shown in Fig. 2(b), the tangent bifurcation (μ l = 1) arises at d ≈ −11.9 and the chaotic state d ≲ −11.9 (λ 1 > 0, λ 2 = 0) appears. This chaotic state produced by tangent bifurcation indicates the alternating laminar and turbulent modes of intermittency chaos in a general way 46 ; this dynamic was demonstrated in our previous work 34 . That is, the intermittency route to chaos exists in this region.
Signal response in chaotic resonance. In the above mentioned chaotic parameter regions #1 and #2, we evaluated the response against a weak signal (A = 10 −2 , f 0 = 0.1). To begin with, we compared the cycle histograms F t ( ) between periodic and chaotic states. As shown in Fig. 3, in the cases of both region #1 (a) and region #2 (b), F t ( ) in the periodic state (solid line) does not fit Ĩ t ( ) ext (dotted line) because the periodic response against Ĩ t ( ) ext induces growth in its values at specific bins. On the other hand, F t ( ) in the chaotic state fits Ĩ t ( ) ext according to a chaotic response with scatter timing against I ext (t). This tendency can also be observed in their C(τ) as shown in Fig. 3(c) and (d). That is, C(τ) becomes approximately 0 in the periodic state, but C(τ) exhibits a sinusoidal shape in the chaotic state. In the following evaluations, we use max τ C(τ) to characterize the signal response, because the sinusoidal shape of C(τ) with period T 0 can be identified by amplitude and lag corresponding to τ τ C max ( ) and its τ value. Furthermore, this signal response is evaluated using τ τ C max ( ) and λ 1 . Figure 3  Scientific RepoRts | 7: 1331 | DOI:10.1038/s41598-017-01511-y dependence of τ τ C max ( ) (upper) and λ j (j = 1, 2) (lower) on parameter d in regions #1 and #2, respectively. In region #1 (Fig. 3(a)), the neuron exhibits the periodic spiking (λ 1 ≈ 0, λ 2 < 0) in 0.82 ≲ d ≲ 0.88, and the chaotic spiking (λ 1 > 0, λ 2 ≈ 0) in 0.88 ≲ d ≲ 0.92. In the periodic spiking state, the value of τ τ C max ( ) is less than 0.1; whereas in the chaotic spiking state, the value of τ τ C max ( ) is higher in comparison with the periodic spiking state. In particular, at the d ≈ 0.89 location around the bifurcation to chaos, called the edge of chaos 47 below, τ τ C max ( ) has a peak value (≈0.8). This can be interpreted as CR arising in the chaotic region (0.88 ≲ d ≲ 0.92). In region #2 (Fig. 3(b)), the chaotic spiking state (λ 1 > 0, λ 2 ≈ 0) arises in −15.5 ≲ d ≲ −12 and τ τ C max ( ), and is a high value due to the effect of this chaotic spiking state. Also, the value of τ τ C max ( ) indicates a similar tendency for region #1 (Fig. 3(a)), i.e., at the d ≈ −12.3 location around the bifurcation to chaos, τ τ C max ( ) has a peak value (≈0.9).
Furthermore, Fig. 4(a) and (b) show the scatter plots between τ τ C max ( ) and λ 1 obtained in Fig. 3 in the cases of region #1 and region #2, respectively. The red dotted line indicates the mean value of τ τ C max ( ) in bin λ 1 with window Δλ 1 = 0.001. From these results, in both regions, τ τ C max ( ) peaks at the appropriate value of max τ C(τ) ( τ ≈ . τ C max ( ) 0 7 at λ 1 ≈ 0.03 in region #1 and τ ≈ . τ C max ( ) 0 9 at λ 1 ≈ 0.04 in region #2). The points for this appropriate value for λ 1 correspond to the points representing the edge of chaos in Fig. 3. That is, the signal response in CR has a unimodal maximum with respect to the stability for chaotic orbits represented by λ 1 , and this peak is localized at the edge of chaos.
Signal sensitivity in the edge of chaos region. In the edge of chaos, i.e., the chaotic state near the bifurcation point, the power spectrum for the time series of system behavior has several peaks. In the periodic bifurcation route to chaos, the trajectory is restricted to the narrow space around the multiple-periodic trajectory before the points at the bifurcation to chaos. Therefore, the power spectrum of the chaotic state inherits the peaks from the power spectrum of the multiple-periodic state, while in the intermittency route to chaos, the laminar state dominates in the time series of system behavior. Hence, the power spectrum has peaks near the frequency components of the laminar state. The upper panels of Fig. 5 show the power spectrum of v(t) under the signal-free condition in the edge of chaos in region #1 (a) (d = 0.896) and #2 (b) (d = −12.0). For the reasons described above, the power spectrum has several peaks. Furthermore, as shown in the lower panels of Fig. 5, resonant frequency/amplitude zones and points (max τ C(τ) > 0.5) indicated by the red line and black points, respectively, in the signal-adapted condition. Here, its frequency f 0 corresponds to the horizontal line of the upper panels of Fig. 5. From this result, the resonant zones have a tendency to distribute near the peaks for the power spectrum in the signal-free condition. This is especially significant with the weaker signal amplitude regions.

Discussion and Conclusions
We showed herein two distinct routes to chaos, the period-doubling bifurcation route and the intermittency route, by using the Lyapunov exponent with a saltation matrix and index for stability of a fixed point on the Poincaré section. Furthermore, under the condition of receiving input from a weak periodic signal, the enhancement of the signal response by the effect of chaotic spikes (chaotic resonance) was confirmed in the chaotic regions induced by the above routes to chaos. Specifically, in both chaotic regions, the signal response in CR had a unimodal maximum with respect to the stability for chaotic orbits represented by λ 1 . Thus, it can be interpreted that the instability of the chaotic orbit in CR plays a role of the noise strength in SR.
Furthermore, we have confirmed that the peak of the signal response was located in the edge of chaos. There, we identified the periodic components in chaotic spiking activity as shown in Fig. 5. In the case of a relatively large signal strength, we found broadening of the signal frequency region in which the efficient signal response was high. On the other hand, in the case of a small signal strength, the region of high efficiency was restricted to the immediate neighborhoods of frequencies for the periodic components in chaotic spiking. This characteristic of signal response in relation to signal strength and frequency, called Arnold's tongue, is widely observed in synchronization phenomena 48 . Therefore, the high efficiency of signal responses in the edge of chaos could be interpreted as synchronization between the input signal and the periodic component in chaotic spiking activity.
In future work, we intend to evaluate the signal response in CR in large-sized spiking neural networks.