Noise-tuned bursting in a Hedgehog burster

Noise can shape the firing behaviors of neurons. Here, we show that noise acting on the fast variable of the Hedgehog burster can tune the spike counts of bursts via the self-induced stochastic resonance (SISR) phenomenon. Using the distance matching condition, the critical transition positions on the slow manifolds can be predicted and the stochastic periodic orbits for various noise strengths are obtained. The critical transition positions on the slow manifold with non-monotonic potential differences exhibit a staircase-like dependence on the noise strength, which is also revealed by the stepwise change in the period of the stochastic periodic orbit. The noise-tuned bursting is more coherent within each step while displaying mixed-mode oscillations near the boundaries between the steps. When noise is large enough, noise-induced trapping of the slow variable can be observed, where the number of coexisting traps increases with the noise strength. It is argued that the robustness of SISR underlies the generality of the results discovered in this paper.


. Introduction
How can we change the spike counts of a bursting neuron? From the viewpoint of dynamical systems, a burst consists of a series of spikes followed by a period of quiescence, and the spiking activity is caused by the fast currents while modulated by the slow currents through a variety of bifurcations (Izhikevich, 2007). For a single unit, varying the bifurcation parameter or the timescale separation parameter may lead to changes in spike counts of a burst (Wang, 1993;Barrio and Shilnikov, 2011;Zhu and Liu, 2017). For coupled bursting neurons, the effects of time delay and coupling can also induce various bursting patterns with different numbers of spikes in each burst (Dhamala et al., 2004;Zhu and Liu, 2017;Jia et al., 2018) or even termination of bursting such as amplitude death (Thottil and Ignatius, 2017).
Noise is inevitable and may play functional roles in neuronal systems (Lindner et al., 2004;Faisal et al., 2008;McDonnell and Ward, 2011;Schmerl and McDonnell, 2013;Bauermann and Lindner, 2019). For excitable neurons close to bifurcation, noise of moderate strength can give rise to coherent oscillations. This noise-induced oscillatory phenomenon was named coherence resonance (CR) (Hu et al., 1993;Longtin, 1997;Pikovsky and Kurths, 1997). For noise acting on the fast variable, DeVille et al. (2005) and Muratov et al. (2005) found a stochastic resonance-like coherent behavior and referred .
to it as self-induced stochastic resonance (SISR). Both CR and SISR can be significantly relevant to the dynamical behaviors in neuronal systems (Yamakou and Jost, 2019;Touboul et al., 2020;Yamakou et al., 2020;Baspinar et al., 2021;Zhu et al., 2022). In particular, SISR is proven to be more robust than CR since SISR does not require the system to be close to bifurcation (DeVille et al., 2005). Further, SISR can occur in a wider range of noise strengths and different noise strengths can induce different coherent stochastic orbits (DeVille et al., 2005;Zhu and Nakao, 2021). This property is highly relevant to our question, that is, noise may also tune the spike counts in a bursting neuron. The mathematical mechanism of the SISR oscillator was first analyzed by Freidlin regarding the motion of a light particle in the force field perturbed by small noise (Freidlin, 2001). Later, Muratov et al. (2005) showed its ubiquity in randomly perturbed excitable systems and termed it as SISR. The large timescale separation is important for the realization of SISR in that larger timescale separation enables the coherent oscillations to occur in larger intervals of the noise strength Yamakou and Jost, 2018). In the limit of vanishing noise and large timescale separation of fast-slow systems, the stochastic dynamics exhibit almost deterministic periodic oscillations, which can be predicted by the timescale matching condition (DeVille et al., 2005;Muratov et al., 2005;Yamakou and Jost, 2018;Yu and Liu, 2021). However, this condition does not consider the differences in timescales on different branches, which can be important to the transition positions, e.g., in the FitzHugh-Nagumo (FHN) neuron model (Zhu and Nakao, 2021;Zhu et al., 2022). Recently, we proposed a distance matching condition, which solves this problem by defining the so-called mean first passage velocity (Zhu and Nakao, 2021). In this paper, we use this condition to analyze the bursting neuron whose potential difference between the slow manifold and the barrier is non-monotonic. Our analysis reveals the possibility to control the spike counts in the bursting neurons solely by the noise strength.
. Materials and methods

. . Hedgehog burster
We consider a modified version of the FHN neuron model, which can exhibit bursting behaviors in two dimensions. This model is named Hedgehog burster by Izhikevich (2000) due to the hedgehog-like limit-cycle attractor. The governing equations are as follows: where f (x, y) = x − x 3 3 − y + 4 L(x) cos(40y), g(x, y) = x + a, and x, y represent the fast membrane potential and the slow recovery variable, respectively. The timescale separation is characterized by the small parameter ε = 0.0001. The bifurcation parameter a = −0.2 is fixed so that the fixed point is unstable and there is a stable limit cycle as shown in Figure 1A. The Hedgehog burster described by Equation (1) is different from the original FHN model because of the cosine term cos(40y) with the logistic function L( The logistic function L(x) is chosen so that the left branch of the x nullcline remains nearly unchanged while the right branch fluctuates, which is key to the realization of the bursting behavior. The large timescale separation forces the Hedgehog limit cycle to stick to the slow manifolds (the left and right branches of the x nullcline). As such, the wavy right branch of the x nullcline determines the spike counts in each burst (six spikes in Figure 1B), which can be also observed in the time series of y(t) as in Figure 1C. It can also be observed that there is a closed loop of the x nullcline on top of the right branch (in fact, there are more such loops above). Initial states on this loop remain on it for a while and then jump to the left nullcline, giving a single spike. However, the dynamics related to this loop are transient, which are out of the scope of this paper.

. . Stochastic periodic orbits by SISR
Noise can not only modify the system's deterministic behavior, but can also induce new phenomena that are not observed in deterministic systems. DeVille et al. (2005) and Muratov et al. (2005) showed that noise acting on the fast subsystem can give rise to a stochastic resonance-type phenomenon, which was named SISR. In SISR, the system's state follows the slow manifold while making an almost deterministic transition at a critical position so that the resulting stochastic oscillation can be very coherent. SISR resembles the stochastic resonance (SR) phenomenon in that both SR and SISR exhibit barrier-crossing behaviors. However, the potential difference is modulated by the external periodic signal in SR while by the slow variable in SISR.
The noise-induced coherent orbits caused by SISR are completely different from the deterministic limit cycle, as we will see later for the Hedgehog burster. Here, we show for the Hedgehog burster with non-monotonic potential differences, the SISR can still be predicted by using our previously proposed distance matching condition (Zhu and Nakao, 2021).
The governing equation for the noise-perturbed Hedgehog burster is as follows: Frontiers in Computational Neuroscience frontiersin.org . /fncom. .  where f (x, y) and g(x, y) are the same as in Equation (1), ξ (t) is the Gaussian white noise satisfying ξ (t) = 0, and , where σ represents the noise strength. Figures 5, 6 illustrate typical noise-tuned bursting dynamics of the system (Equation 2) obtained by the Monte Carlo simulations with several noise strengths. Depending on the noise strength, the system exhibits different stochastic periodic orbits.
To theoretically predict the stochastic periodic orbit of the SISR phenomenon, we follow the similar procedure as in Zhu and Nakao (2021).

. . Noise-tuned bursting in hedgehog burster
In what follows, we will apply the distance matching condition (Zhu and Nakao, 2021) to predict the transition positions on the stable branches of the x nullcline. The distance matching condition compares the noise-induced displacement of the state away from the stable branch and the distance .
/fncom. . from the stable branch to the unstable one. When these two distances are equal to each other, transition happens with a large probability.
When the state is on the left branch of the x nullcline, for each fixed y, there is a corresponding mean first passage time (MFPT) T e obtained from the Kramers rate (Kramers, 1940;Gardiner, 1985): where the potential function U(x; y) = − f (x, y)dx + C (C is a constant and y is regarded as a fixed parameter). The double prime over the potential function represents the second-order derivative with respect to x. Denoting the potential function on the left, middle, and right branches as U l , U m , and U r , respectively, the potential difference between the middle and left branches, dU ml = U m − U l , and the middle and right branches, dU mr = U m − U r , can be easily calculated. The landscape of the potential function (the constant C can be safely set to zero as it is eliminated in Equation 3) and the potential differences are shown in Figure 2. The potential difference dU ml monotonically decreases for decreasing y. However, the potential difference dU mr largely fluctuates while decreases on average for increasing y. Thus, the potential differences have more than one intersection points, which is significantly different from the FHN neuron model with monotonic potential differences (Zhu and Nakao, 2021). The intersection points in Figure 2B give the values of y at which the two potential differences are equal to each other. At these points, the boundary crossing of x from the left branch to the middle branch and that from the right branch to the middle branch are equally difficult. From the MFPT, the mean first passage velocity (MFPV) function can thus be defined as in Zhu and Nakao (2021): where S(y) = x m (y) − x l (y) is the distance between the middle and the left branches of the x nullcline for fixed y. The distance matching condition for the transition from the left branch to the middle branch is as follows: for the system state starting from y 0 (t 0 ) and reaching the critical transition position y * (t * ), the total displacement in the x direction from the left branch can be obtained by integrating the MFPV over time, which should be equal to the distance S(y) between the branches. By using Equations (3) and (4), we can express the distance matching condition in a self-consistent manner (Zhu and Nakao, 2021) as where we have changed the variable of integration from t to y by using dy = ε(x + a)dt. By solving Equation (5), we can obtain the critical transition position y * . Figure 3A plots the integral on the left-hand side (lhs) and the distance on the right-hand side (rhs) of Equation (5) for the left branch as functions of y for different noise strengths. The intersection points between the lhs (blue) and rhs (black) correspond to the critical transition position y * . The starting .
/fncom. . position of the system state is chosen as y 0 = 0.221, which is the maximum value of y for the deterministic limit cycle on the left branch shown in Figure 1A. We will see later that the starting position is not important as long as it is far enough from y * . The critical transition position on the left branch is plotted against the noise strength in Figure 4A (blue triangles). As expected, the transition position gradually increases with the noise strength.
The critical transition position on the right branch can be similarly calculated. However, due to the non-monotonic potential difference as shown in Figure 2, the transition process along the right branch should be treated differently. To this end, we divide the right branch into six regions as shown in Figure 1A and apply the distance matching condition separately. In each region, we start the system state from the bottomrightmost position on the branch (where the potential is locally minimum) and seek the self-consistent solution to Equation (5), where the lower limit y 0 of integration is taken as the y coordinate of the starting position. This procedure is started over again in every region. The results on the right branch are illustrated in Figure 3B for different noise strengths, where the first intersection of the lhs (red) and rhs (black) gives the critical transition position. It should be noted that the intersection takes place only on the left half of each region before the minimum of the distance curve (rhs, black). This means that if the transition cannot occur before the position with the minimum potential difference, then the transition after that position is extremely difficult in each region. Figure 4A also displays the obtained critical transition positions for various noise strengths on the right branch (red circles), which shows a significantly different tendency from those on the left branch (blue triangles). The transition positions exhibit several stages and decrease stepwisely with the noise strength, corresponding to the different regions on the right branch in Figure 1A. In Figure 4A, the transition positions on the left and right branches intersect at (σ , y * ) ≈ (0.173, −0.253), which corresponds to the case that the transitions occur at the same position (y * l = y * r ≈ −0.253). Above this noise strength, i.e., σ > 0.173, the predicted transition positions cannot constitute a complete orbit and are meaningless.
In numerical simulations, further increasing the noise strength (σ > 0.173) will make the actual transition positions on the left and right branches asymptotically approach each other (but never cross). The crossing at (σ , y * ) ≈ (0.173, −0.253) of the predicted transition positions in Figure 4A is due to the fixed initial position y 0 . A more appropriate value for the starting position y 0 is to choose the critical transition position on the other branch. We modify the starting positions on the left branch to the transition positions on the right branch and obtain the corrected y * as in Figure 4A. The modified results almost coincide with the original ones except when the two transition positions on the left and right branches get too close. This is because only a small interval of y contributes to the transition process, as can be seen from the steep curves in Figure 3. Therefore, y 0 can be chosen arbitrarily above but not close to the transition position (Zhu and Nakao, 2021).
After computing the transition positions on each branch, the complete stochastic periodic orbit can be accordingly obtained by gluing the slow motions along the slow manifolds at these transition positions. Figure 5 plots the predicted stochastic periodic orbits for several representative noise strengths that correspond to the black dashed lines labeled as S1-S4 in Figure 4A. The theoretical predictions are in good agreement .
/fncom. . with the results of Monte Carlo simulations. Even for large enough noise (S4), the range of y can be well predicted despite the relatively poor prediction on x. Increasing the noise strength will decrease the size of the stochastic periodic orbit. Therefore, it can be seen that noise can tune the spike counts in every burst at different stages in Figure 4A. From S1 to S4, there are 6,5,3, and 1 spikes in each burst. Considering the large timescale separation, the period of the bursting orbit can be approximated by the motion along the slow manifolds on the left and right branches, which is given as where y cl and y cr denote the critical transition positions on the left and right branches, respectively. The predicted periods are also in good agreement with the simulation results in different transition stages as shown in Figure 4B. The small standard deviation within each stage implies highly coherent oscillations, which can be observed in Figure 5. We can see that for each noise strength, the transition occurs near the leftmost tip of the right branch of x nullcline in each region. There is only one transition position on the right branch, which we call single-escape transition. On the other hand, the period near the boundary between two different stages has a larger standard deviation. The overestimation and underestimation of the bursting period near the boundaries imply the emergence of the mixed-mode oscillations. Three typical noise-induced mixed-mode oscillations are illustrated in Figure 6. In contrast to Figure 5, the system state can transit at two tips of the right branch, which we call double-escape transition. Indeed, for the . /fncom. . noise strength near the boundary between the two stages of coherent oscillations, the critical transition position begins to switch from one region to the other (i.e., the first intersection point switches from one valley to the other of the black distance curve in Figure 3B). Therefore, for these noise strengths, the transition can occur at different tips on the right branch for different realizations.

. . Large noise-induced slow variable traps
As discussed in the previous section, for large noise, the critical transition positions on the left and right branches approach each other asymptotically. For the FHN model, we have shown in Zhu and Nakao (2021) that, as the noise strength increases, the size of the stochastic periodic orbit will shrink and the period will approach zero. This is also true in the present case as we can see from Figure 4B or Figure 5. Furthermore, for sufficiently large noise, there appears a stable region where the slow variable y is almost clamped at a fixed value in the steady state after the transient in the FHN system (Zhu and Nakao, 2021) (see also Touboul et al., 2020 for large noise-induced asynchrony in interacting neuronal ensembles).
For the Hedgehog burster investigated in this paper, the nonmonotonic potential difference enables more than one stable regions where the slow variable y is restricted in a small interval, which we call slow variable traps. Due to the large noise strength, the transition occurs in a small window of the slow variable y (below the small window, the transition from the left branch is instant and above it the transition from the right branch . /fncom. . is instant). In Figure 7, the slow variable trap phenomenon is shown for three large values of the noise strength. For σ = 0.5, 0.65, and 0.8, there are 3, 4, and 6 stable slow variable traps, respectively. The widths of the orbits, i.e., the fluctuations of the slow variable in the traps are slightly different from each other. The thinnest orbit in the left panel of Figure 7 is the yellow one near y = −0.25, which corresponds to the first intersection point of the potential difference curves dU ml and dU mr in Figure 2B (or the crossing of the transitions positions on different branches in Figure 4A). For upper traps, the potential difference dU ml is large and it requires longer interval for the transition from the left branch to occur; while for lower traps, the potential difference dU mr is large and it requires longer interval for the transition from the right branch to take place. At the middle trap at y ≈ −0.25, the transitions from the left and right branches are balanced, resulting in smaller fluctuations in y.
In Figure 7, it is interesting to note that before the state becomes trapped into the stable region, the state remains in another metastable region for certain time in the case with σ = 0.5 (cyan and green). This kind of transient trap phenomenon may also occur in coupled SISR systems and influence the collective dynamical behavior in large ensembles.

. Discussion
We have investigated the noise-tuned bursting in a Hedgehog burster, where the noise is applied to the fast variable. Through the proposed distance matching condition, the stochastic periodic orbits were well predicted and the results of Monte Carlo simulations were reproduced. It was found that increasing the noise strength leads to the shrinking of the orbits. This implies that the spike counts in each burst can be tuned by varying the applied noise strength. Moreover, mixed-mode oscillations with the double-escape transition from the right branch can be realized near the boundaries between different stages. Finally, large noise-induced slow variable traps were analyzed, where the number of stable traps depends on the noise strength. The coexistence of multiple traps for the Hedgehog burster is a notable phenomenon that cannot be observed in the FHN model that permits only one trap when large enough noise is applied on the fast variable.
. /fncom. . It should be noted that although the Hedgehog burster considered in this paper is in the oscillatory case, noisetuned bursting can also occur in the excitable situation as shown in Figure 8 (see Appendix A). In that case, a lower bound of the noise strength exists for the coherent oscillations to occur, which depends on the bifurcation parameter a. Larger values of a lead to higher transition positions on the left branch as shown in Figure 8, which illustrates that the bifurcation parameter in the excitable Hedgehog burster also has a significant influence on the noisetuned bursting.
Despite that the Hedgehog burster studied in this paper is quite artificial, it shares similar dynamics with other high-dimensional bursting neurons, such as the Hindmarsh-Rose neuron model. It should be noted that the SISR mechanism allows coherent oscillations to occur away from local bifurcations, which leads to robustness of the phenomena investigated in this paper to parameter variations. Therefore, it is promising to observe these phenomena in other fast-slow dynamical systems experimentally. Besides, neurons behave collectively in vivo. It has been shown that SISR can play an important role in the interacting excitable FitzHugh-Nagumo systems with synchronization and anticoherence (Touboul et al., 2020). As synchronization manifests its significance in both normal and pathological (e.g., in Parkinson's disease and epilepsy) cases, the control of bursting patterns on a single neuron and the ensembles related to the contents investigated in this paper is of great interest. These open problems will be investigated in our future works.