Routes to Chaos Induced by a Discontinuous Resetting Process in a Hybrid Spiking Neuron Model

Several hybrid spiking neuron models combining continuous spike generation mechanisms and discontinuous resetting processes following spiking have been proposed. The Izhikevich neuron model, for example, can reproduce many spiking patterns. This model clearly possesses various types of bifurcations and routes to chaos under the effect of a state-dependent jump in the resetting process. In this study, we focus further on the relation between chaotic behaviour and the state-dependent jump, approaching the subject by comparing spiking neuron model versions with and without the resetting process. We first adopt a continuous two-dimensional spiking neuron model in which the orbit in the spiking state does not exhibit divergent behaviour. We then insert the resetting process into the model. An evaluation using the Lyapunov exponent with a saltation matrix and a characteristic multiplier of the Poincar’e map reveals that two types of chaotic behaviour (i.e. bursting chaotic spikes and near-period-two chaotic spikes) are induced by the resetting process. In addition, we confirm that this chaotic bursting state is generated from the periodic spiking state because of the slow- and fast-scale dynamics that arise when jumping to the hyperpolarization and depolarization regions, respectively.

the piecewise-constant and piecewise-linear systems as types of piecewise-smooth systems [19][20][21][22][23] . Tsubone et al. 23 proposed a systematic method using an analytical approach to predict the parameter regions for the chaotic states in piecewise-constant systems. Meanwhile, Mitsubori and Saito 19 and Nakano and Saito 20 developed a systematic method for piecewise-linear systems. However, the analysis of chaos and bifurcation in piecewise-smooth systems such as the Izhikevich neuron model, which generally include non-linear terms, requires the evaluation of the Lyapunov exponents and characteristic multipliers against exhaustive parameter sets. Several indices considering the effect of the resetting process using the saltation matrix have been proposed 16,[24][25][26][27] . Coombes et al. 16 utilised this approach and conducts an analysis of the planar non-linear integrate-and-fire model in a large parameter region. They revealed that the parameter region for chaotic states is located at the boundary between burst firing and fast spiking. The Izhikevich neuron model clearly features various types of bifurcation and routes to chaos under the effect of the state-dependent jump in the resetting process [27][28][29][30][31] . However, neither the relation between the chaotic behaviours and the state-dependent jump nor the mechanism for inducing the chaotic states using the resetting process has been revealed.
Revealing this mechanism requires evaluating the influence of the state-dependent jump on the trajectory in a continuous system through a comparison of two systems: versions with and without the resetting process. The aforementioned Izhikevich neuron model cannot be applied for this purpose because of its divergent behaviour in the spiking state when the resetting process is removed. We previously introduced a preliminary approach to this issue based on a comparison of spiking neuron models such as the FitzHugh-Nagumo neuron model with and without the resetting process 32,33 . In actual neural systems, the system state transits from resting to spiking through both Hopf and saddle-node bifurcations 1,34 . However, the FitzHugh-Nagumo neuron model permits spiking activity only via Hopf bifurcation. A non-linear equation for the membrane recovery variable with a sigmoidal function has been reported to generate spiking activity via both types of bifurcation 34 .
In our approach, we first adopt a continuous two-dimensional (2D) spiking neuron model with a non-linear equation for the recovery variable in which the orbit in the spiking state does not exhibit divergent behaviour. We then add the resetting process to the model. Utilising a rigorous method to analyse the bifurcation and chaos in hybrid systems (i.e., using a characteristic multiplier of the Poincaré map and Lyapunov exponents with a saltation matrix), we evaluate several routes to chaos, which cannot be achieved in a hybridised FitzHugh-Nagumo neuron model 33 , by changing the resetting process parameters and comparing the structure of the attractor between the system versions with and without the resetting process.

Model and Method
Spiking Neuron Model with the Resetting Process. The FitzHugh-Nagumo neuron model 10,11 is driven by 2D ordinary differential equations with the following form: where v and u represent the membrane potential of a neuron and the membrane recovery variable, respectively. The parameter a determines the shape of a v-nullcline =  v ( 0), while b/c and c exhibit the sensitivity of u and its time constant, respectively. In real-world neural systems, the system state transits from resting to spiking through both Hopf and saddle-node bifurcations 1,34 . However, the FitzHugh-Nagumo model only permits spiking activity via Hopf bifurcation. Spiking neuron models such as the Morris-Lecar neuron model 35 can be used to generate spiking activity via both types of bifurcation. Such models apply the following non-linear equation for  u with a sigmoidal function for v 34,36 : where α is the time constant for u and parameters β and ε determine the shape of the sigmoidal function. Instead of the linear equation given by Eq.
(2), we use Eq. (3) with the parameters set to (a, α, ε) = (0.1, 0.1, 0.05) as the equation for  u. For the continuous spiking neuron model above, we implement the resetting process given by the following equation: where v r and d represent the after-spike reset values of the membrane potential v and the recovery variable u, respectively. In other words, we consider it to be a spike event if v reaches v peak . The state-dependent jump induced by Eq. (4) converges to a continuous trajectory under the condition v r → v peak and d → 0 in the case which v peak is set to a maximum value of v for the orbit of the continuous spiking neuron model given by Eqs (1) and (3).
We numerically analysed this model in SUNDIALS, our non-linear differential/algebraic equation solver simulator, using the backward-differentiation formula method with Newton's iteration 37 . As this function can detect intersection points between a trajectory and the user-defined manifolds, we utilise it to detect the points for v ≥ v peak . The time step for numerical integration is related to the time precision needed to detect intersection points. As the fixed width in this solver is not user-tunable, relative and absolute tolerances are set instead; we set these to sufficiently small values (10 −14 ) to achieve sufficient numerical precision for the evaluation of a chaotic spiking activity and the detection of the intersection. to quantify the chaotic activity in the version of the spiking neuron model with the resetting process 25,26,28 . The evolution of the orthogonal vectors of perturbation l j (j = 1, 2) for Eqs (1) and (3) in a system with a continuous trajectory in spike intervals between the i-th and (i + 1)-th times (t i ≤ t ≤ t i+1 ) is described as follows: where Λ i+1 is the matrix (l 1 , l 2 ). The Jacobian matrix J for Eqs (1) and (3) is given as follows: The saltation matrix at t = t i is given by the following equation: represent the values of (v, u) before and after spiking, respectively. Λ k (T k+1 , T k ) (k = 0, 1, …, N−1) can be expressed as follows in case spikes arise in the range [T k : T k+1 ] 28 : The initial perturbation Λ 0 (T 0 , T 0 ) and the evolution period τ = T k+1 − T k are set to unit matrix E and 0.1, respectively.
The Lyapunov spectrum λ j is calculated as follows based on the norm of l j where the evolution period for λ j is set to Nτ = 10 5 (N = 10 8 , τ = 10 −3 ).
Poincaré map. We set a Poincaré section Ψ(v = v peak ) to conduct a bifurcation analysis in a system with a state-dependent jump. The dynamics of the system behaviour on Ψ are given by the Poincaré map as follows: i l l i We evaluate here the profile of ψ l on a return map of u i − u i+l . In practice, we solve Eqs (1), (3) and (4) against the initial values of (v r , u 0 ), and obtain the u values that pass through Ψ at times l as u l . The values (u 0 , u l ) are plotted on a return map of u i − u i+l .
Characteristic multiplier of the Poincaré map. The variational equations of Eqs (1) and (3) in a system with a continuous trajectory in spike intervals between the i-th and (i + 1)-th times (t i ≤ t ≤ t i+1 ) are defined as follows: where Φ indicates the state transition matrix. Φ(t, t 0 ) can be expressed as follows for the case in which spikes arise in the range [t : t 0 ]: The initial state transition Φ(t 0 , t 0 ) is set here as the unit matrix E. We calculate the characteristic multiplier of the solution for the l period (u 0 = ψ l (u 0 )) as follows: l l 0 where u 0 indicates the initial value of orbit x 0 = (v r , u 0 ) at t = t 0 . The projection P and embedding P −1 to and from the Poincaré section are defined as follows by local sections Π 0 and Π 1 , respectively: where w indicates the local coordinate on Ψ. According to the literature 27 , local sections Π 0 and Π 1 are set as follows to solve Eq. (17): where q 0 (v, u) and q 1 (v, u) are the scalar functions used to determine the local sections. Eq. (17) can be developed as follows to utilise the aforementioned sets (see the detailed derivation from Eq. (22) to Eq. (23) in ref. 27 ): |μ l < 1|, μ l = −1 and μ l = 1 represent the stable condition, period-doubling bifurcation and tangent bifurcation, respectively.

Results
Bifurcation in a Continuous 2D Spiking Neuron Model. First, we demonstrate the system behaviour of the continuous spiking neuron model given by Eqs (1) and (3) for the cases of β = 0.5 and β = 0.3, which are hereinafter called regions #1 and #2, respectively. For region #1, Fig. 1 and vector field of   v u ( , ) in the case I = 0; these are denoted by a dotted line, a dashed line, and arrows, respectively. In this region, the fixed points (i.e. the points at which the v-nullcline intersects with the u-nullcline) are located at (v, u) ≈ (0, 0), (0.10, 0), (0.35, 0.06). The system trajectories are attracted to the stable fixed point ((v, u) ≈ (0, 0)). The other intersection points are unstable fixed points. In the I = 0.004 case (upper panel of Fig. 1(c)), the orbit indicated by a solid line exhibits a limit cycle along the vector field because of the effects of the destruction of the pair of stable and unstable fixed points. Next, as shown in the lower panel of Fig. 1(c), a spiking activity appears in the time series of v(t). In region #2, the fixed point is located at (v, u) ≈ (0, 0), and, as the point is stable, the system trajectories are attracted to it ( Fig. 1(b)). Meanwhile, in the I = 0.04 case, a limit cycle emerges as a result of the destabilisation of the fixed point ( Fig. 1(d)).
We evaluate the eigenvalues m j (j = 1, 2) of J around the stable fixed point ((v, u) ≈ (0, 0)) to specify the bifurcation against I. Figure 2 Fig. 2(b). For region #2, the right panel of Fig. 2(a) shows the dependence of m max Re( ) j j on I. In this region, the fixed point is stable in  − . ≤ . I 0 005 0 0193 but becomes unstable when  . I 0 0193 because the complex values of m j , m 1 and m 2 are complex conjugates and pass through a real axis at I ≈ 0.0193 (right panel of Fig. 2(b)). Therefore, this bifurcation is classified as Hopf bifurcation.  Figure 4 shows the dependence of the bifurcation diagrams of u i and λ j (j = 1, 2) on the v r parameter in regions #1 and #2 with the value of parameter d fixed at d = 0.01 in Fig. 3. This result shows the chaotic states, which exhibit an irregular behaviour of u i , and the Lyapunov exponents λ 1 > 0, λ 2 = 0 in region #1 observed in the parameter set for (   .
The upper panel of Fig. 5(a) shows the time series of v(t) as an example of the chaotic spiking pattern in region #1 (v r = 0.33, d = 0.01). In this spiking pattern, v(t) exhibits two types of behaviours after the resetting process (spike). In one case, v(t) enters the hyperpolarization mode and decreases to v(t) ≈ 0.25. In the other case, v(t) increases to v peak but does not do so through hyperpolarization. This spike repeats several times. In the former behaviour, (v, u) jumps in the region of <  v 0 through the v-nullcline in the v-u phase plane, as shown in the lower panel in Fig. 5(a). Meanwhile, in the latter behaviour, (v, u) jumps in the region of >  v 0 but not do so through the v-nullcline via the resetting process. This spiking pattern is called a burst and is observed in actual neural systems. Figure 6 shows the bifurcation diagram after the value of u i + d is reset (a) and typical orbits at v r = 0.25, 0.3, 0.32 (b) for the transition scheme from spike to burst against changing v r . This result implies that (v,u) jumps in the region of < Fig. 6(b)). Meanwhile, at  . v 0 31 r , (v,u) jumps in the region of >  v 0 as well as in the region of <  v 0 (red circle in the v r = 0.32 graph in Fig. 6(b)). The upper and lower panels of Fig. 5(b) show examples of the chaotic time series of v (t) and the behaviour of (v, u) in the v-u phase plane in region #2 (v r = 0.14, d = 0.01). In this case, (v, u) always jumps in the region >  v 0 and the orbit exhibits a near-period-2 chaotic behaviour. We can also confirm the aforementioned mechanisms for generating a chaotic bursting behaviour and a near-period-2 chaotic behaviour in the Izhikevich neuron model 29, 31 .
Next, we investigate the dependence of the return map on v r . For region #1, Fig. 7(a) shows ψ(u i ) on the u i+1 − u i return map in cases with the resetting process, in which v r = 0.15, v r = 0.33 (corresponding to the value of v r in Fig. 5(a)) and v r = 0.395, and in the case without the resetting process. In the case without the resetting process, which is indicated by the dotted black line, ψ(u i ) exhibits a nearly constant value (≈0.01). Meanwhile, in  Fig. 5(a) can be considered to be induced by this structure. The stretching and folding structure then fades as v r further decreases, as in the v r = 0.15 case indicated by the solid red line. For region #2, Fig. 7(b) shows ψ (u i ) on the u i+1 − u i return map in the cases with the resetting process, in which v r = 0.12, v r = 0.14 (corresponding to the value of v r in Fig. 5(b)), and v r = 0.20, and in the case without the resetting process. ψ (u i ) is a nearly constant value in the case without the resetting process (≈0.03), which is indicated by the dotted black line, as well as for region #1. However, the stretching and folding structure arises as an effect of the resetting process (v r = 0.20 case, indicated by the solid blue line). Meanwhile, the frequency of the stretching and folding structure increases in the case in which the jump distance  is larger than v r = 0.20 (v r = 0.14 case indicated by the solid green line. The near-period-2 chaotic spiking activity can be interpreted as being produced by this structure. This frequency decreases as v r further decreases (v r = 0.12 case indicated by the solid red line).
Fixing v r at 0.33 in region #1 and v r at 0.14 in region #2, we evaluate the dependence of ψ(u i ) on the u i+1 − u i return map on parameter ε (Fig. 8). The stretching and folding structure with piecewise nearly linear maps appears, as shown in ε = 0.053, in the case in which ε increases in region #1 (corresponding to the shape of the u-nullcline approaching the step function). A stretching and folding structure with non-linearity emerges at u i ≈ 0.04 when the ε values decrease, such as at ε = 0.05, 0.03. In region #2, nine folding structures appear at ε = 0.05. The frequency of folding in this case decreases as the value of ε decreases (see ε = 0.045, 0.04).

Discussion and Conclusion
In this paper, we applied the resetting process to a continuous 2D spiking neuron model with a sigmoidal nullcline structure to reveal the mechanisms for the emergence of chaotic states in a hybrid spiking neuron model. We also evaluated the bifurcation and routes to chaos against two types of spikes generated by the parameter sets for saddle-node bifurcation (region #1) and spikes for Hopf bifurcation (region #2) by changing the value of the resetting parameter. Through an evaluation using the Lyapunov exponent with a saltation matrix and the index for the fixed-point stability on a Poincaré section, we demonstrated that two types of chaotic behaviour are induced by the resetting process.
A chaotic state with bursting characteristics emerged in region #1 through tangent bifurcation, with a jump distance that increased as v r decreased. This chaotic state moved to the periodic state through period-doubling bifurcation as the distance increased even further. We investigated the dependence of the return map on the  resetting parameter and compared the models with and without the resetting process. Consequently, we found that the non-linear stretching and folding structure of the attractor is induced by the resetting process. This structure is also affected by the shape of the sigmoidal function of the u-nullcline. We also confirmed that the bursting chaotic states emerged according to the adjustment of this structure's non-linearity.
Homoclinic orbits coexisting with fast and slow dynamics are needed to reproduce bursting 38 . For example, using the Hindmarsh-Rose neuron model as a continuous three-dimensional spiking neuron model can reproduce the bursting by fast dynamics of the membrane potential and recovery variable and the slow dynamics of the bursting variable 12 . Meanwhile, in hybrid spiking neuron models (Fig. 6) bursting can be reproduced using 2D systems because the hyperpolarization after resetting to <  v 0 and the depolarization in the inter-burst term after resetting to >  v 0 can play the roles of slow and fast dynamics, respectively. The analysis of the Hindmarsh-Rose neuron model indicated that a unimodal peak on the return map of the Poincaré section existed in chaotic bursting [39][40][41] . This structure, which was also confirmed in region #1 of our model (Fig. 7(a)), contributes to the generation of chaotic bursting. Investigations of this chaotic bursting by inter-spike intervals (ISI) have previously been described in the literature 42,43 . Gu 43 reported that the return map of the ISI exhibited a unimodal peak in chaotic bursting, with a peak value that decreased with the approach of chaotic spiking in a physiological experiment involving a neural pacemaker. Meanwhile, Innocenti et al. 42 demonstrated that this return map has two peaks in chaotic bursting and that these peaks merged with the approach of chaotic spiking in the Hindmarsh-Rose neuron model. The patterns of the ISI prior to the last-minute hyperpolarization might be interpreted as reflecting the number of peaks in the return map. In our case, we confirmed that the trend in our result is consistent with that of Gu (Additional Information) 43 . A chaotic state with a near-period-2 behaviour emerged in region #2 through tangent bifurcation as the jump distance increased. As the distance was increased even further, this chaotic state moved to a periodic state through tangent bifurcation. Evaluation of the return map showed that this chaotic state emerged from the non-linear stretching and folding structure of the attractor induced by the resetting process. The folding frequency increased as the jump distance increased; also, this frequency decreased as the value of parameter ε decreased.
In conclusion, the resetting process provides and enhances non-linear effects in attractors, a result that cannot be achieved in continuous spiking neuron models of less than two dimensions. Chaotic states tend to arise when a state-dependent jump exists with an appropriate distance. The two types of chaotic behaviour and bifurcation mentioned above can also be observed in the widely used Izhikevich neuron model 14,29,31 . Therefore, the effects induced by the resetting process revealed in this study might be utilised to generate various chaotic spiking patterns.
Further research based on this study should be undertaken to classify the types of transition from chaotic bursting to spiking and to evaluate the bifurcation and chaos in neural networks composed of hybrid spiking neurons.   Figure 9 shows the inter-spike intervals (ISI) ISI i = t i+1 − t i , where t i indicates the spike time (i = 0, 1, 2,…) in the cases with the parameter settings for chaotic activity (v r = 0.325, 0.33, 0.34) in region #1. As a result, the return map of the ISI exhibits a unimodal peak in chaotic bursting, and its peak value decreases with the approach of the parameter region of periodic bursting and spiking (Fig. 6).