DELAY-INDUCED MIXED-MODE OSCILLATIONS IN A 2D HINDMARSH-ROSE-TYPE MODEL

. In this study, we investigate a Hindmarsh-Rose-type model with the structure of recurrent neural feedback. The number of equilibria and their stability for the model with zero delay are reviewed ﬁrst. We derive conditions for the existence of a Hopf bifurcation in the model and derive equations for the direction and stability of the bifurcation with delay as the bifurcation parameter. The ranges of parameter values for the existence of a Hopf bifurca- tion and the system responses with various levels of delay are obtained. When a Hopf bifurcation due to delay occurs, canard-like mixed-mode oscillations (MMOs) are produced at the parameter value for which either the fold bifurcation of cycles or homoclinic bifurcation occurs in the system without delay. This behavior can be found in a planar system with delay but not in a planar system without delay. Therefore, the results of this study will be helpful for determining suitable parameters to represent MMOs with a simple system with delay.

1. Introduction. In many biological systems, the effects of delays are important in understanding the dynamics of neural networks [4]. In one review paper [14], the authors noted that the analysis of computational models at different space-time scales shows the fundamental mechanisms that relate neural processes to neuroscience data. Importantly, modeling at the microscopic level is necessary because neural information is transmitted between the neurons of the brain. Delay factors are inherent in many biological systems for two reasons: (1) the finite propagation speed of the signals and (2) the finite processing time in the synapses. The time delay is caused by the finite speed of signal transmission over a distance in processing information from various sensory systems to form a coherent and unified perception of the external world. Currently, relatively few studies have been conducted on biological systems with time delay. In mathematical biology, the effect of time delay on dynamical behaviors, such as stability, firing mode, and mixed-mode oscillations (MMOs), in a single neuron model remains unclear. Hence, it is important to study the behavior of a single neuron with a differential-difference model. chemical systems [36]. Several researchers have studied the bifurcation scenarios of MMOs. Koper [27] found that the MMOs occurred near two Hopf bifurcations, which were located near the saddle-node (SN) bifurcation in a cross-shaped phase diagram. These results differed from those of previous studies that showed MMOs occur in the vicinity of a Shil'nikov homoclinic orbit. Erchova and McGonigle [19] reviewed studies on MMOs in the brain. These studies included examples at the single-neuron level and possible instances of MMOs across local and global brain networks using a number of noninvasive techniques, such as electroencephalography and magnetoencephalography. Rubin and Wechselberger [39] investigated the generation of MMOs in a three-dimensional Hodgkin-Huxley neuron model. Krupa et al. [28] studied the generation of MMOs in a coupled system using a two-dimensional oscillator to model the dynamics of the membrane potential and the calcium concentration in dopaminergic neurons in the mammalian brain stem. Desroches et al. [15] investigated the organization of the MMOs in the self-coupled FitzHugh-Nagumo system and in a simplified model to study synchronization in a network of Hodgkin-Huxley neurons. However, relatively few studies have investigated MMOs in a two-dimensional, single-neuron model in a differential-difference form. Therefore, we chose to investigate the occurrence of MMOs in the 2DHR-type model with recurrent feedback.
The remainder of the paper is organized as follows. In Sec. 2, we introduce a general class of nonlinear neuron models in a differential-difference form. The model can be considered a generalization of the 2DHR-type model [42,10]. From the characteristic equation of this model, we analyze the local stability of an equilibrium as a function of the time delay. In Sec. 3, the direction and stability of the Hopf bifurcations are analyzed by reducing the system on a center manifold. In Sec. 4, we provide results of numerical simulations to assess the accuracy of the theoretical analysis. In addition, we suggest a possible reason for the emergence of MMOs in the planar system. Finally, conclusions are presented in Sec. 5.
2. The delay Hindmarsh-Rose-type model. Let us consider that the simplest possible modification of the systems described in [42,10] to simulate synaptic feedback is the following system of differential-difference equations: where x and y denote the cell membrane potential and a recovery variable, respectively, and a, b, c, d, τ and k 1 are positive constants. The parameter c represents the time scale. The parameter k 2 is positive for excitatory feedback and negative for inhibitory feedback, with the strength of the feedback given by the magnitude of k 2 . Because the current is ionic, we expect that its magnitude will be approximately proportional to the difference between x and the "resting potential" v 0 , which is similar to the assumption in [37]. The parameter I denotes the membrane current or an external stimulus. If the value of I is increased and the other parameters are unchanged, the cubic function is shifted up. However, if the value of the parameter a is reduced and the other parameter values are unchanged, the quadratic function is shifted down. Hence, the effect of I is reflected though parameter a.

SHYAN-SHIOU CHEN AND CHANG-YUAN CHENG
For convenience, the term I − k 2 v 0 is replaced with I, and it is assumed that k 1 + k 2 = 1 and k 1 = k in Eqs. (1) and (2). System (1)-(2) can be recast in the following form.
For k = 1, the system is the same as in [42,10]. We can identify the number of equilibria and their stability and several codimension-one and codimension-two bifurcations, such as the SN, Hopf, Bautin and Bogdanov-Takens bifurcations. Let the point (x 0 , y 0 ) be an equilibrium, and let x(t) = x(t) − x 0 and y(t) = y(t) − y 0 . To simplify the notation, the variable x is replaced with x, and the equations can be transformed as follows: The associated characteristic equation is or Eq. (7) can be expressed as where F (λ) = λ 2 + a 1 λ + a 0 and G(λ . This characteristic equation determines the local stability of an equilibrium. That is, the equilibrium is stable if and only if all of the characteristic roots λ have negative real parts.
The following two conditions determine whether the characteristic roots are purely imaginary.

SHYAN-SHIOU CHEN AND CHANG-YUAN CHENG
Substituting Eq. (23) into Eq. (24), the left side of Eq. (24) can be expressed as , and the right side can be expressed as Therefore, So, Similarly, Therefore, W 20 (θ), W 11 (θ), E 1 and E 2 can be determined. We use the following definitions [22]:  (4) with k = 1. For k = 1 and τ = 0, we will obtain the range of parameter values for which a Hopf bifurcation exists and discuss possible reasons for the occurrence of MMOs. In the following, we define b = 2 and c = 3, as in [10].
To begin, we recall several results for k = 1 from [10]. This is the case without any delay. For k = 1 and τ = 0, the system is the same as for k = 1. In Fig. (2), the four lines AH − , AH + , SN − and SN + separate the parameter domain (a, d) into    Figure 2. The number and stability of equilibria in the a-d parameter space with b = 2, c = 3, and τ = 0. The number and stability of the equilibria for the six areas Ω 1 -Ω 6 are as follows: Ω 1 has one stable equilibrium, Ω 2 has two stable equilibria and one saddle equilibrium, Ω 3 has one unstable equilibrium, Ω 4 has one stable, one unstable and one saddle equilibria, Ω 5 has one stable equilibrium, and Ω 6 has two unstable equilibria and one saddle equilibrium. The line SN − indicates a SN bifurcation for a saddle and an unstable node, the line SN + indicates a SN bifurcation for a saddle and a stable node, and the lines AH ± indicate AH bifurcations. the six areas denoted by Ω 1 -Ω 6 . The number and stability of the equilibria for the six areas are indicated in Fig (2). A codimension-one bifurcation can occur when the values of the parameters (a, d) cross one of these four lines. For example, when the parameter values cross the line AH + from Ω 3 to Ω 5 , the stability of the equilibrium changes. If the transversality conditions hold, this is a Hopf bifurcation. Similarly, when the parameter values cross the line SN − from Ω 1 to Ω 2 , one stable and one saddle equilibria emerge. If the transversality conditions hold, this is a SN bifurcation with a stable node. By [10], the Hopf (resp. Hopf, SN and SN) bifurcation in fact occurs at the point belongs to the line AH − (resp. AH + , SN − and SN + ). The locus of the saddle-homoclinic bifurcation, which is based on the BT bifurcation [29], is represented by the blue dashed line in Fig. (3). The Bautin bifurcation occurs at d = 2b − b 2 /c 2 = 32/9 ≈ 3.56. Hence, for d < 32/9, the Hopf bifurcation is subcritical, and for d > 32/9, the Hopf bifurcation is supercritical. A nondegenerate fold bifurcation of the cycles originating from the Bautin point occurs to the right of AH + and lower than d < 32/9; this condition is represented by the blue line in Fig. (3). Therefore, the number and stability of equilibria, the conditions for codimension-one bifurcations for SN ± , AH ± , the homoclinic bifurcations based on the BT bifurcation and fold bifurcations of cycles are known for the case k = 1.
Next, we examine the parameter values that lead to the existence of Hopf bifurcations and switching stability in Eqs. (5) and (6) with τ as the bifurcation parameter. For k > 1, the model has the form of recurrent excitation, and for k < 1, the model possesses the structure of recurrent inhibition. For these two recurrent feedback cases, we choose two value, k = 1.1 and k = 0.9. The parameters corresponding  The parameter regions corresponding to one pair of purely imagery roots ω + and two pairs of purely imagery roots ω ± are shown in pale gray and dark gray, respectively, in Fig. (3(a)). As the value of k approaches 1, the two regions reduce to the lines SN ± and AH + , respectively, in the a-d parameter domain. The regions of switching stability are shown in Fig. (4), which shows the maximum number of stability switches for k = 1.1, where the region for one switch is shown in blue and the region for two switches is shown in pale blue. Regions with greater numbers of switches are shown in black near the boundary. Therefore, the parameter regions for Hopf bifurcations with τ as the bifurcation parameter in the parameter space (a, d) are known for the system (5) and (6) for both recurrent inhibition and recurrent excitation.
To demonstrate the stability switches of the equilibria and dynamical behavior of the Hopf bifurcations, we will show results for two pairs of parameters, labeled ROI 1 and ROI 2 in Fig. (3(a)), where we have listed the Hopf bifurcations in Table 1. For the first pair, two switches, from stability to instability and back to stability, occur. For the second pair, only one switch, from stability to instability, occurs. The dynamical behavior for the parameter values d = 2.55 and a = 0.2 (i.e., ROI 2 ) is presented in Fig. (5); the results for ROI 1 are similar and are thus not shown. It follows from Fig. (2) that the system with the given parameter values has only one stable equilibrium (x 0 , y 0 ), where x 0 ≈ −0.968 and y 0 ≈ −0.666, when τ = 0. One switch from stability to instability occurs when a = 0.2, where delay parameters are listed in Table 1. For these parameters, we obtain β 2 ≈ 6.81 > 0 and µ 2 ≈ −26.7 < 0 at τ + 0 ≈ 1.7377. It follows from Theorem 3.1 that the Hopf bifurcation at τ + 0 is  Figure 4. The maximum number of stability switches in the ad parameter space with b = 2, c = 3 and k = 1.1. The colored regions indicate the existence of a pair of purely imaginary roots, λ ± = iω ± with ω + > ω − > 0, for Eq. (10). In other words, there exists a stability switch. The dark blue, blue, and green regions indicate those values for which the maximum number of stability switches is 1, 2 and 3, respectively. The number of stability switches is greater than 4 in the red region.
subcritical. That is, we can determine the existence of an unstable limit cycle. The equilibrium (x 0 , y 0 ) is stable if τ < τ + 0 and unstable if τ > τ + 0 . Time histories of the membrane potentials for Eqs. (5) and (6) with various delays, i.e., τ = 1.6, 1.7, 1.75 and 1.79, are plotted in Fig. (5). The trajectory converges to an equilibrium in Fig. (5(a)). In Fig. (5(b)), there exist a stable MMO, an unstable limit cycle and a stable equilibrium. The existence of the unstable limit cycle is confirmed by Theorem 3.1. That is, there exist bistable behaviors. In Fig. (5(c)), the system has a stable MMO, and the equilibrium becomes unstable. In Fig. (5(d)), the equilibrium becomes an unstable attractor, and the MMO decays. Therefore, the numerical results are consistent with the theoretical results for the Hopf bifurcations given by Theorem 3.1, and MMOs are exhibited by the model. The trajectories of the system with τ = 1.7 and 1.75 are canard-like, as can be observed in Fig. (6).
The range of parameter values leading to the occurrence of MMOs is indicated with the dotted black line in Figs. (3(a)) and (3(b)). We consider only the parameter values −5 ≤ a ≤ 5 and −2 ≤ d ≤ 5, as shown in Fig. (3(a)). The range of delay resulting in MMOs is between 0 and 5. For example, with the parameter values a = 0.2 and d = 2.55, 7 different MMOs are shown, i.e., from 1 1 to 1 7 , in Fig. (7(a)). Fig. (7(a)) shows two regions, "Non-Periodic" and "Non-MMO". The former is the area where only one stable equilibrium exists, and the latter is the area where there is an unstable equilibrium and a stable periodic orbit. In Fig. (7(b)), the fold bifurcation of cycles occurs at τ c ≈ 1.6616. The first Hopf bifurcation occurs at τ ≈ 1.7377, as shown in Table 1. The parameter pairs (a, b) shown with black dots in the parameter space are those for which MMOs exist. It follows from [10] that for the system without delay, the bifurcation of the limit cycle occurs near AH + and to the right of AH + , and the bifurcation of the homoclinic orbit occurs between SN + and SN − in Fig. (3(a)). When the Hopf bifurcation with delay τ as the bifurcation parameter occurs for certain values of the parameters a, b, c, d, k and τ , the MMOs can occur in two cases with τ = 0: the parameter is near the fold bifurcation of the limit cycle, or the parameter is near the bifurcation of the homoclinic orbit.
Finally, we suggest a possible reason for the occurrence of MMOs. Eqs. (1) and (2) can be considered a fast-slow (i.e., singularly perturbed) system because the system involves two variables with dynamics on different timescales. The ratio between the two timescales is given by 1/c 2 . For example, if we take c = 3, the ratio is 1/9. Hence, we can call the variables x and y fast and slow variables, respectively. The set of points that satisfy g(x, y) = 0 is a slow manifold of the differential-difference equation in the singular case. One solution of the singularly  τ max x(t)   1 (b). In (a), 7 types of the MMOs are shown, i.e., from 1 1 to 1 7 , where 1 1 , 1 2 and 1 3 are labeled between the curves and the remaining four labels are omitted because of space limitations. In (a), the region labeled "Non-periodic" indicates where the system has a stable equilibrium; in the "Non-MMO" region, there exist a stable period orbit and an unstable equilibrium. The black line indicates the locus for which the first Hopf bifurcation occurs. In (b), the abscissa is τ and ordinate is the maximum of the action potential.
perturbed system is a canard that follows an attracting slow manifold and then follows a repelling slow manifold near a bifurcation point. MMOs may also occur through the canard phenomenon [3]. The canard-induced MMO is a limit cycle resulting from a Hopf bifurcation that transitions from a small, nearly harmonic cycle to a large relaxation oscillation in a narrow parameter interval. For a twodimensional system with additional variables or noise, MMOs may occur in larger regions as the dynamics switch between small-and large-amplitude oscillations. In our study, the canard-induced MMOs of the two-variable differential-difference equation may occur in a Hopf bifurcation with delay.

5.
Conclusions. In this study, we analyzed Hopf bifurcations in a 2DHR-type model with recurrent feedback, where delay was the bifurcation parameter. This model is a generalization of the model presented in [42,10]. Depending on the parameters of the model, the number of equilibria and their stability may be the same as those of the model without delay [10]. The conditions for the existence of purely imaginary roots of the characteristic equation for the system were studied, and the ranges of parameter values for these conditions were obtained for two cases, k = 0.9 and k = 1.1. When the parameter values were near the values corresponding