Analysis of Chaotic Resonance in Izhikevich Neuron Model

In stochastic resonance (SR), the presence of noise helps a nonlinear system amplify a weak (sub-threshold) signal. Chaotic resonance (CR) is a phenomenon similar to SR but without stochastic noise, which has been observed in neural systems. However, no study to date has investigated and compared the characteristics and performance of the signal responses of a spiking neural system in some chaotic states in CR. In this paper, we focus on the Izhikevich neuron model, which can reproduce major spike patterns that have been experimentally observed. We examine and classify the chaotic characteristics of this model by using Lyapunov exponents with a saltation matrix and Poincaré section methods in order to address the measurement challenge posed by the state-dependent jump in the resetting process. We found the existence of two distinctive states, a chaotic state involving primarily turbulent movement and an intermittent chaotic state. In order to assess the signal responses of CR in these classified states, we introduced an extended Izhikevich neuron model by considering weak periodic signals, and defined the cycle histogram of neuron spikes as well as the corresponding mutual correlation and information. Through computer simulations, we confirmed that both chaotic states in CR can sensitively respond to weak signals. Moreover, we found that the intermittent chaotic state exhibited a prompter response than the chaotic state with primarily turbulent movement.


Introduction
By virtue of recent developments in brain measurement technology, it is now recognized that information is transmitted among neurons according to not only their firing rate but also their spike timing. Therefore, spiking neuron models, which describe spike timing, have attracted considerable attention. The Hodgkin-Huxley (HH) approach [1] is known to be the most useful spiking neuron model for simulating neurodynamics, and does so by describing the capacitance of membranes and the characteristics of the resistance of ion channels. This model is represented by four equations involving several physiological parameters pertaining to membrane potential, the activation of the Sodium (Na) and Potassium (K) ion-currents, and the inactivation of the Na ion-current. This model can reproduce almost all spiking activity observed in neural systems by tuning the parameters. However, because of the complexity of the physiological parameters involved, researchers have proposed many neuron models that use a smaller number of parameters, and hence are simpler than the HH model, by focusing on membrane potential behavior of spiking modes and response characteristics during spiking activity, such as the integrate-and-fire neuron model and the FitzHugh-Nagumo neuron model [2]. Of such simplified models, the Izhikevich neuron model [3] combines continuous spike-generation mechanisms and a discontinuous resetting process following the spikes, and can reproduce major spike patterns that have been observed experimentally by tuning a few parameters. Furthermore, the variety of the spiking properties obtained through this model is greater than those obtained through other models [4]. The Izhikevich neuron model can simulate chaotic spiking activity under certain conditions regarding specific parameter sets and given adequate numerical precision by integrating small time steps and accurately detecting the spiking points [4,5]. Furthermore, it has been claimed that in order to determine this chaotic state, the Lyapunov exponents need to be calculated using a saltation matrix in order to account for the state-dependent jump in the resetting process [5]. Note that if the effect of the state-dependent jump is ignored, i.e., if the saltation matrix is not used, the maximum Lyapunov exponent has a positive value even in periodic states [5].
In the last few decades, researchers have realized that several processes for signal detection and transmission in neural systems are supported by the mechanism of stochastic resonance (SR) [6][7][8]. SR is a phenomenon whereby the presence of noise helps a nonlinear system amplify a weak (sub-threshold) signal. SR can be observed in many kinds of systems that have three ingredients: a kind of barrier/threshold, a source of noise, and a weak input signal. Chaotic processes also cause a phenomenon similar to SR called chaotic resonance (CR) with the ingredient of a deterministic fluctuating activity instead of a source of noise in stochastic processes. As regarding CR, two kinds of fluctuating activities have been considered, so far. One is the case whereby external additive chaotic signal is applied to the system instead of stochastic noise [9,10]. The other is the case whereby external additive chaotic signal is absent and alternatively intrinsic chaotic activities are utilized. Recently, there have been many studies of CR in the latter condition. The characteristic of this CR was initially investigated using simple models [11][12][13][14] such as a one-dimensional cubic map. By focusing on the neural system, researchers observed that chaos exists at several hierarchical levels, from the electrical response of a single neuron to the activity of the entire brain as an assembly of neurons [15][16][17][18]. CR has recently been studied in neural systems, such as by using the chaotic neural network [19,20] and the inferior olive (IO) neural system [21][22][23][24]. In [19,20], the superior signal response capability of CR to conventional SR was exhibited by using a chaotic neural network model based on the mean firing rate of neurons. Furthermore, the work in [21][22][23][24] suggested that CR plays a part in the function that allows IO neurons to transmit error signals containing large amounts of information for cerebellar learning in continuous spiking neuron models. In the CR phenomenon in spiking neural systems, chaotic behavior leads to the generation of spikes not at specific times, but at varying scatter times for each trial as input signals. Thus, the frequency distribution of these spike timings against the input signal becomes congruent with the shape of the input signal [21][22][23][24]. However, no study to date has investigated and compared the signal responses of some chaotic states in CR revealed by the bifurcation analysis of a spiking neural system with a resetting process, such as the Izhikevich neuron model. In past research, we have analyzed the bifurcation and the signal response of CR in the Izhikevich neuron model [25][26][27]. However, these studies involved issues related to the accuracy of the numerical simulations [25,26] and the quantitative bifurcation analysis [27].
In this paper, we examine and classify the chaotic characteristics of the Izhikevich neuron model by using Lyapunov exponents with a saltation matrix and Poincaré section methods in order to address the measurement problem caused by the state-dependent jump in the resetting process. Following the structure adopted in our past work [25][26][27], we then rigorously evaluate the signal response in CR for classified chaotic states through a numerical verification method [28] and a quantitative method to specify the bifurcation in the system with the resetting process [29].

Izhikevich neuron model
The Izhikevich neuron model [3,4] is a two-dimensional (2D) system of ordinary differential equations of the form v 0 ¼ 0:04v 2 þ 5v þ 140 À u þ I; ð1Þ In the above equation, v and u represent the membrane potential of a neuron and the membrane recovery variable, respectively.  the time evolution of v(t) and the system trajectory in a phase plane (v, u), respectively. Due to the resetting process, when v(t) exceeds 30 [mV], the system state (v, u) (discontinuously) jumps to point ((v, u) % (−65, 0.3)), as shown in Fig 1 (b). In our simulation, we numerically analyzed this model through non-linear differential/algebraic equation solvers by using the backward differentiation formula method [28] to achieve sufficient numerical precision in order to evaluate chaotic spiking activity. This method is more precise than Euler's method, which was adopted to reproduce only periodic spiking in [3].
To determine the uniformity of the neuron spikes, we adopt the coefficient of variation for inter-spike intervals [30]: T k is the k-th spike interval (T k = t k+1 − t k ), and Var(T k ) and < T k > are the variance and the mean of T k , respectively. CV here becomes 0 in the one-period state and positive in the nonperiodic state, including chaotic states.
As an index to check whether a state is chaotic, the maximum Lyapunov exponent is ordinarily used for systems with continuous trajectories, and is calculated by Here, jd k (t l = 0)j (k = 1, 2, Á Á Á, N) are N perturbed initial conditions applied to the system trajectory at t l = 0, and d k (t l = τ) represent their evolution in time for t l 2 [0: τ] [31]. Let us consider the case where a neuron fires at t l = t s , k = i. Since the time evolution of the system's trajectory is discontinuous in the resetting process, d i (t s ) receives the interruption, and jd i (τ)j is rendered irrelevant to the evolution of the system's trajectory. Due to this influence, λ 1 loses its accuracy in such situations. Therefore, trials for new measures of such a system are needed [29,32,33]. One such proposed measure is the insertion of the Lyapunov exponent into a saltation matrix in order to address the stability of the system's trajectory, including state-dependent jumps [5]. Therefore, we use the Lyapunov exponent with a saltation matrix to analyze chaotic states in the Izhikevich neuron model.

Fundamental properties of the model
The Izhikevich neuron model can reproduce major firing patterns, such as regular spiking, intrinsically bursting, chattering, and fast spiking [3,4]. Moreover, research has suggested that this model can simulate chaotic behavior with appropriate parameter values (a = 0 (1) and (2)  ). The dynamics of (u 1 , u 2 , Á Á Á, u N ), which is the evolution of u over time on C, is defined as a Poincaré mapping u i+1 = ψ(u i ). As shown in Fig 2(c), on the return map (u i , u i+1 ), the solution for u i+1 = ψ(u i ), the orbit of u i and u i+1 = u i are indicated by dotted, solid, and dashed lines, respectively. It has been observed that the orbit of u i exhibits chaotic behavior in the range −102 ≲ u i ≲ −90, and the shape of ψ displays a stretching and folding structure as a feature of the non-linear map.
To quantify this chaotic activity in the Izhikevich neuron model, we use the Lyapunov exponent with a saltation matrix. 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 Eqs (1) and (2) are defined as follows: where, F, 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 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 : can be expressed as By using the eigenvalues l k j ðj ¼ 1; 2Þ of F k (T k+1 ,T k ), the Lyapunov spectrum λ j is calculated by 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 takes 1000 [ms] before 20 spikes occur. We investigated the behavior . Furthermore, the system came to rest (non-firing) (λ 1 < 0, λ 2 < 0) for I ≲ −104.5, whereas periodic firing (λ 1 % 0, λ 2 < 0) was observed at I ≳ −94.5.
By fixing the value of I at −99, we investigated the bifurcation of this system by replacing I with d by using a bifurcation diagram consisting of (u 1 , u 2 , Á Á Á, u N ). Fig 4(a), (b), and (c) represent the bifurcation diagrams of u i , λ j , and CV, respectively, as functions of d. For d ≲ −11.9, the chaotic trajectory was distributed in the range −103 ≲ u i ≲ −80, and the system exhibited a chaotic state (λ 1 > 0) and irregular spiking activity (CV > 0), excluding periodic windows (λ 1 = 0), given periodic-1 (CV = 0) and multiple periodic (CV > 0) states. As the value of d increased, those of λ 1 and CV decreased. They converged at 0 for d ≳ −11.9, i.e., the system assumed periodic states and exhibited periodic spiking. In order to conduct bifurcation analysis in the system with a state-dependent jump, we used the evaluation method intended to assess the stability of a fixed point u 0 = ψ l (u 0 ) (l = 1,2, Á Á Á) on a Poincaré section. In the literature [29], this stability has been evaluated by In the above, u 0 = (v 0 , u 0 ) indicates the initial value of orbit u = (v, u) at t = t 0 . jμ < 1j, μ = −1, and μ = 1 represent the stable condition, period doubling bifurcation, and tangent bifurcation, respectively. In Fig 4, the tangent bifurcation at l = 2 arises at d % −11.9. Through this tangent bifurcation, the system transitions into chaos at d ≲ −11.9, as shown in Fig 4(a) and (b). Furthermore, in the region d % −11.9 as a bifurcation point, the trajectory (v, u) and the time series of v(t) were examined. At d = −11 in Fig 5(a), the time series of v(t) (left) and the trajectory (right) indicate periodic spiking and a one-period state, respectively. As the value of d decreases, the behavior of the system becomes irregular, as shown in Fig 5(b), (c), and (d). It can be observed that the durations of the apparently periodic instances of spiking seemed to have decreased during episodes of chaotic behavior in the system with a reduction in the value of d. . In order to focus on the oscillatory behavior of u i , we used the return map (u i , u i+2 ). Fig 6 (right) shows the orbit of u i (solid line), and the solutions to u i+2 = ψ 2 (u i ) (dotted line) and u i+2 = u i (dashed line). When d = −11 ((a)), the orbit of u i remained at the intersection (% (−98.5, −98.5)) of u i+2 = ψ 2 (u i ) and u i+2 = u i , and there were two unstable fixed points on both sides of this stable fixed point at u i % −101.5 and −91.5. As noted above, the tangent bifurcation arose at d % −11.9, i.e., a pair consisting of an unstable fixed point and a stable one destroyed each other. Through this tangent bifurcation, at d = −12 ((b)), the orbit of u i exhibited sluggish movement (called laminar state) in the region where the slope of ψ 2 was approximately 1.0 (−102 ≲ u i ≲ −94), and irregularly active movement (called turbulent or burst state) in regions with larger slopes () 1). Such chaotic, dynamic alternation between laminar and turbulent states is called intermittency chaos [34,35]. Note that the term of burst is not used in this paper to avoid confusion between the chaotic movement and the neural spike patterns in neurodynamics such as intrinsically bursting and chattering bursting. As the value of d decreased, the area of the region producing the laminar state, where the slope of ψ 2 was approximately 1.0, shrunk as well. The turbulent state was subsequently dominant in the dynamics, i.e., the system transitioned from intermittency chaos to chaos involving primarily turbulent movement, as shown in Fig 6 (b)-(d).
We thus confirmed that in the parameter region (a = 0.2, b = 2, c = −56, d = −16, I = −99) proposed by Izhikevich as the parameter set that produces chaotic behavior, the periodic state transitions into chaos through tangent bifurcation and intermittency chaos, i.e., the intermittent route to chaos exists in this region.

Response efficiency in chaotic resonance
In this section, in order to reveal the efficiency of signal response in a chaotic state, we investigated the response of the system to a weak periodic signal that could not be detected in a periodic state. First, we introduced an extended Izhikevich neuron model with a signal. Second, to

Analysis of CR in Izhikevich Neuron Model
To quantify the signal response, we used the following indices, i.e., Eqs (13) and (14) and Eqs (15)- (17). The mutual correlation C(τ) between the cycle histogram FðtÞ of the neuron spikes and the signal SðtÞ is given by C SF ðtÞ ¼< ðSðt þ tÞÀ < SðtÞ >ÞðFðtÞÀ < FðtÞ >Þ > : For the time delay factor τ, we checked max τ C(τ), i.e., the largest C(τ) between ÀT 0 2 t T 0 2 . As an extensively used index for evaluating information transmission, we used mutual information, which is information transmitted from input S to output F, MIðF; SÞ ¼ HðFÞ À HðFjSÞ: Here, S and F consisted of m s (s 1 , s 2 , Á Á Á, s m s ) and m f (f 1 , f 2 , Á Á Á, f m f ) event states, into which SðtÞ (−A * A) and FðtÞ (0 * its maximum value) were equally divided, respectively. H(F) and H(FjS) are given by HðFjSÞ ¼ À where P(s j ) and P(f j ) represent the occurrence probability of s j and f j , and P(f j js i ) is the conditional probability for the occurrence of f j and s i . In our simulation, we set m s = m f = 20. However, if the maximum value of FðtÞ was smaller than that of m f , m f was set to the maximum value because FðtÞ was the integer of firing counts.

Dependence on parameter d
This section concerns the response of the system represented by Eq (11). As mentioned in the section "Fundamental properties of the model, " the system featured periodic firing activity in the region (−12 ≲ d ≲ −5) and chaotic activity, which was generated through the intermittent route to chaos, in the region (−17 ≲ d ≲ −12). We now examine the behavior of the system in the Izhikevich neuron model with a weak sinusoidal signal S(t) (A = 0.3). Fig 7 shows the bifurcation diagram of u i ((a)), λ j ((b)), and CV ((c)) as a function of d. In the region −17 ≲ d ≲ −12, the behavior of u i exhibited chaotic activity (λ 1 > 0). However, as d increased, λ 1 and CV decreases; the system was entrained by S(t) (λ 1, 2 < 0), and exhibited a period-2 state at −12 ≲ d ≲ −11.5. In the region −11.5 ≲ d ≲ −5, the system exhibited a periodic state (λ 1 % 0, λ 2 < 0) and u i showed a slight motion when the range of u i was % 1.

Sensitivity of signal response in chaotic resonance
In the section "Fundamental properties of the model, " we observed that the chaotic state with primarily turbulent movement and the intermittent chaotic state coexisted in the region −17 ≲ d ≲ −12. In this section, we examine the sensitivity of signal response in the region of parameter d, including in the two chaotic states. Fig 11 shows    particular, the chaotic states along the boundary between the chaotic and periodic states, called the edge of chaos [19,36], exhibited the highest sensitivity and the promptest response among all chaotic states. Note that the distribution of points satisfied with the promptness was localized to a small boundary region in the region with high sensitivity.
We show the exact dependence between bifurcation and signal response given a weaker signal (A = 0.01) than the one used in the section "Dependence on parameter d, " which was set at approximately the strength of a signal to be detected almost on the edge of chaos. Fig 12 shows the bifurcation diagrams of u i ((a)), λ j (j = 1, 2) ((b)), CV ((c)), max τ C(τ) ((d)), and MI(F; S) ((e)) as functions of d. u i exhibited chaotic and irregular activity (λ 1 > 0, CV % 0.5), and periodic movement with slight motion in its range u i % 0.1 in the regions −17 ≲ d ≲ −12 and −12 ≲ d ≲ −5, respectively. The signal response in the region −17 ≲ d ≲ −14, which exhibited the impressive performance (max τ C(τ) % 0.9 and MI(F; S) % 1.7) at A = 0.3 (see Fig 10), degraded such as at max τ C(τ) ≲ 0.7 and MI(F; S) ≲ 1. Nonetheless, the signal response at the edge of chaos maintained satisfactory performance (max τ C(τ) % 0.9 at d % −12.19 and MI(F; S) % 1.6 at d % −12.5) and rapidness (jτj % 0.1 [ms]). Further, Fig 13 shows the relationship between max τ C(τ) and λ 1 in the region −13.5 d −11 of Fig 12. The red dotted line indicates the mean value of max τ C(τ) in the bin λ 1 with window Δλ 1 = 0.001. From this result, we see max τ C(τ) recorded a peak (% 1.0) at λ 1 % 0.04, i.e., we confirmed that signal response in CR has an unimodal maximum with respect to the degree of stability for chaotic orbits.
Dependence on signal frequency f 0 Finally, we evaluated the dependence of signal response on signal frequency in CR under the condition A = 0.01, d = −12.19. As shown in Fig 14, the dependence of max τ C(τ) ((a)) and MI (F; S) ((b)) on signal frequency f 0 recorded a peak at f 0 % 0.103 [kHz] with chaotic state (λ 1 > 0 ((c))). Thus, CR has a resonance frequency, as is the case with resonance phenomena in general.

Conclusion
In this paper, we examined the chaotic characteristics of the Izhikevich neuron model in detail by using Lyapunov exponents with a saltation matrix and Poincaré section methods, and discovered two distinctive states: a chaotic state with primarily turbulent movement and an intermittent chaotic state. In order to evaluate the signal response of CR in these classified states, we introduced an extended Izhikevich neuron model by considering a weak periodic signal, and defined a cycle histogram of neuron spikes, and the corresponding mutual correlation and information. Through computer simulations, we confirmed that both chaotic states in CR can sensitively respond to weak signals, and that the intermittent chaotic state exhibited a significantly prompter signal response than the chaotic state with primarily turbulent movement. In particular, the sensitivity and rapidity at the edge of the chaos, located along the border of the intermittency chaos and the periodic state, recorded the highest values of all other chaotic states. Furthermore, we confirmed that signal response in CR is dependent on the frequency of signal input, as is the case with resonance phenomena in general.
From the results obtained here, we expect that the Lyapunov exponent with a saltation matrix can be applied to other reset systems with a state jump, such as control systems and models in the social and financial sciences, as an index to determine whether they are chaotic. With regard to the CR phenomenon, we believe that CR plays an important role in information transmission in the nervous systems. From an engineering viewpoint, the high signal response efficiency in CR can be utilized for the development of devices to detect weak signals.
Another subject of research based suggested by this study is the evaluation of signal response in neuron assemblies consisting of Izhikevich neurons. We are currently exploring the coupling strength dependency of signal response in neuron assemblies.