# Noise-tuned bursting in a Hedgehog burster

^{1}School of Mechanical Engineering, Nanjing University of Science and Technology, Nanjing, China^{2}Department of Systems and Control Engineering, Tokyo Institute of Technology, Tokyo, Japan

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.

## 1. 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 (Muratov et al., 2005; 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.

## 2. Materials and methods

### 2.1. 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-\frac{{x}^{3}}{3}-y+4L(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(40*y*) with the logistic function $L(x)=\frac{1}{1+exp\left[5(1-x)\right]}$. 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.

**Figure 1**. Phase portrait and time series of the Hedgehog burster (Equation 1). **(A)** Phase portrait. The blue curves and the red line represent the *x* and *y* nullclines, respectively. The black curve with arrows illustrates the stable limit cycle, where the dots indicate the states on the limit cycle with equal time intervals. **(B,C)** Time series of the membrane potential *x*(*t*) and the recovery variable *y*(*t*).

### 2.2. 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:

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 〈ξ(*t*)ξ(*t*′)〉 = σδ(*t*−*t*′), 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).

## 3. Results

### 3.1. 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 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*)d*x*+*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.

**Figure 2. (A)** Potential function *U*(*x*; *y*), where *y* is regarded as a fixed parameter. The red curve is the *x* nullcline, which follows the local extrema of the potential for each fixed *y*. **(B)** Potential differences between the middle and left branches, *dU*_{ml} = *U*_{m}−*U*_{l} and between the middle and right branches, *dU*_{mr} = *U*_{m}−*U*_{r}, where *U*_{l}, *U*_{m}, and *U*_{r} denote the potential values on the left, middle, and right branches for given *y*, respectively.

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 d*y* = ε(*x*+*a*)d*t*. 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 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.

**Figure 3**. Identification of the critical transition positions by the distance matching condition, Equation (5). **(A)** Left branch. The blue curves and the black curve represent the lhs and rhs of Equation (5), respectively. **(B)** Right branch. The red curves and the black curve represent the lhs and rhs of Equation (5), respectively. Six patches in gray and yellow correspond to those in Figure 1A. Noise strengths for the integration curves of the lhs of Equation (5) are from 10^{−3} to 10^{−0.5} with equal logarithmic intervals [**(A)**: from left to right; **(B)**: from right to left]. The red dashed curves in **(B)** display the result for σ = 10^{−0.5}.

**Figure 4**. Critical transition positions and bursting periods. **(A)** Critical transition positions on the left branch (blue triangles) and right branch (red circles). The black crosses represent the transition positions on the left branch obtained by modifying the lower limit of integration as explained in the text. **(B)** Bursting periods. Theoretical prediction, Equation (6) vs. the Monte Carlo simulations (averaged for *t* = 200 with error bars showing the standard deviation). Noise strengths are the same as in Figure 3.

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 bottom-rightmost 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 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.

**Figure 5**. Stochastic periodic orbits and time series with the single-escape transition for different noise strengths. The bold orange and the thin gray curves are the theoretical periodic orbits and the results of Monte Carlo simulations, respectively. From top to bottom, the parameters are: σ = 0.00455 (S1), σ = 0.0207 (S2), σ = 0.0695 (S3), and σ = 0.16 (S4). They correspond to the black dashed lines S1–S4 in Figure 4A.

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 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.

**Figure 6**. Stochastic periodic orbits and time series with the double-escape transition for different noise strengths. The thin gray curves are the results of Monte Carlo simulations. From top to bottom, the parameters are: σ = 0.013 (D1), σ = 0.04 (D2), and σ = 0.1 (D3). They correspond to the pink dot-dashed lines D1–D3 in Figure 4A. The red and blue arrows in the left panel display the two transition positions, which lead to different spike counts in the burst. The red and blue lines in the right panel illustrate mixed-mode oscillations of the slow variable *y*(*t*).

### 3.2. 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 non-monotonic 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 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.

**Figure 7**. Slow variable traps for large noise. From top to bottom, σ = 0.5, 0.65, and 0.8. Left panels show the stable phase portraits and orbits; right panels show the time series from five different initial conditions.

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.

## 4. 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.

It should be noted that although the Hedgehog burster considered in this paper is in the oscillatory case, noise-tuned 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 noise-tuned bursting.

**Figure 8**. Noise-induced and noise-tuned bursting in excitable Hedgehog bursters. **(A–E)** The bifurcation parameters are: *a* = 1.5, 1.4, 1.3, 1.2, and 1.1, respectively. Black and blue curves represent trajectories for noise strengths above and below the threshold values shown in the inset of **(F)**. Red dots denote the stable equilibrium points. **(F)** Variation of the *y* nullcline with the bifurcation parameter *a*. The inset displays the threshold values of noise strength obtained from the distance matching condition (Equation 5).

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.

## Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

## Author contributions

JZ conceived the study, conducted the analysis and simulations, and wrote the manuscript. HN conceived the study and wrote the manuscript. Both authors contributed to the article and approved the submitted version.

## Funding

JZ acknowledges support from JSPS KAKENHI JP20F40017, Natural Science Foundation of Jiangsu Province of China (BK20190435), and the Fundamental Research Funds for the Central Universities (No.30920021112). HN thanks JSPS KAKENHI JP17H03279, JP22K11919, JP22H00516, JPJSBP120202201, and JST CREST JP-MJCR1913 for financial support.

## Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

## References

Barrio, R., and Shilnikov, A. (2011). Parameter-sweeping techniques for temporal dynamics of neuronal systems: case study of Hindmarsh-Rose model. *J. Math. Neurosci*. 1, 6. doi: 10.1186/2190-8567-1-6

Baspinar, E., Schülen, L., Olmi, S., and Zakharova, A. (2021). Coherence resonance in neuronal populations: mean-field versus network model. *Phys. Rev. E* 103, 032308. doi: 10.1103/PhysRevE.103.032308

Bauermann, J., and Lindner, B. (2019). Multiplicative noise is beneficial for the transmission of sensory signals in simple neuron models. *Biosystems* 178, 25–31. doi: 10.1016/j.biosystems.2019.02.002

DeVille, R. E. L., Vanden-Eijnden, E., and Muratov, C. B. (2005). Two distinct mechanisms of coherence in randomly perturbed dynamical systems. *Phys. Rev. E* 72, 031105. doi: 10.1103/PhysRevE.72.031105

Dhamala, M., Jirsa, V. K., and Ding, M. (2004). Enhancement of neural synchrony by time delay. *Phys. Rev. Lett*. 92, 074104. doi: 10.1103/PhysRevLett.92.074104

Faisal, A. A., Selen, L. P., and Wolpert, D. M. (2008). Noise in the nervous system. *Nat. Rev. Neurosci*. 9, 292–303. doi: 10.1038/nrn2258

Freidlin, M. I. (2001). On stable oscillations and equilibriums induced by small noise. *J. Stat. Phys*. 103, 283–300. doi: 10.1023/A:1004827921214

Gardiner, C. W. (1985). *Handbook of Stochastic Methods: For Physics, Chemistry and the Natural Sciences*. New York, NY: Springer.

Hu, G., Ditzinger, T., Ning, C. Z., and Haken, H. (1993). Stochastic resonance without external periodic force. *Phys. Rev. Lett*. 71, 807–810. doi: 10.1103/PhysRevLett.71.807

Izhikevich, E. M. (2000). Neural excitability, spiking and bursting. *Int. J. Bifurcat. Chaos* 10, 1171–1266. doi: 10.1142/S0218127400000840

Izhikevich, E. M. (2007). *Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting*. London: MIT Press.

Jia, B., Wu, Y., He, D., Guo, B., and Xue, L. (2018). Dynamics of transitions from anti-phase to multiple in-phase synchronizations in inhibitory coupled bursting neurons. *Nonlinear Dyn*. 93, 1599–1618. doi: 10.1007/s11071-018-4279-x

Kramers, H. (1940). Brownian motion in a field of force and the diffusion model of chemical reactions. *Physica* 7, 284–304. doi: 10.1016/S0031-8914(40)90098-2

Lindner, B., García-Ojalvo, J., Neiman, A., and Schimansky-Geier, L. (2004). Effects of noise in excitable systems. *Phys. Rep*. 392, 321–424. doi: 10.1016/j.physrep.2003.10.015

Longtin, A. (1997). Autonomous stochastic resonance in bursting neurons. *Phys. Rev. E* 55, 868–876. doi: 10.1103/PhysRevE.55.868

McDonnell, M. D., and Ward, L. M. (2011). The benefits of noise in neural systems: bridging theory and experiment. *Nat. Rev. Neurosci*. 12, 415–425. doi: 10.1038/nrn3061

Muratov, C. B., Vanden-Eijnden, E., and Weinan, E. (2005). Self-induced stochastic resonance in excitable systems. *Physica D* 210, 227–240. doi: 10.1016/j.physd.2005.07.014

Pikovsky, A. S., and Kurths, J. (1997). Coherence resonance in a noise-driven excitable system. *Phys. Rev. Lett*. 78, 775–778. doi: 10.1103/PhysRevLett.78.775

Schmerl, B. A., and McDonnell, M. D. (2013). Channel-noise-induced stochastic facilitation in an auditory brainstem neuron model. *Phys. Rev. E* 88, 1–13. doi: 10.1103/PhysRevE.88.052722

Thottil, S. K., and Ignatius, R. P. (2017). Nonlinear feedback coupling in Hindmarsh-Rose neurons. *Nonlinear Dyn*. 87, 1879–1899. doi: 10.1007/s11071-016-3160-z

Touboul, J. D., Piette, C., Venance, L., and Ermentrout, G. B. (2020). Noise-induced synchronization and antiresonance in interacting excitable systems: applications to deep brain stimulation in Parkinson's disease Disease. *Phys. Rev. X* 10, 011073. doi: 10.1103/PhysRevX.10.011073

Wang, X. J. (1993). Genesis of bursting oscillations in the Hindmarsh-Rose model and homoclinicity to a chaotic saddle. *Physica D* 62, 263–274. doi: 10.1016/0167-2789(93)90286-A

Yamakou, M. E., Hjorth, P. G., and Martens, E. A. (2020). Optimal self-induced stochastic resonance in multiplex neural networks: electrical vs. chemical synapses. *Front. Comput. Neurosci*. 14, 62. doi: 10.3389/fncom.2020.00062

Yamakou, M. E., and Jost, J. (2018). Coherent neural oscillations induced by weak synaptic noise. *Nonlinear Dyn*. 93, 2121–2144. doi: 10.1007/s11071-018-4311-1

Yamakou, M. E., and Jost, J. (2019). Control of coherence resonance by self-induced stochastic resonance in a multiplex neural network. *Phys. Rev. E* 100, 022313. doi: 10.1103/PhysRevE.100.022313

Yu, Q., and Liu, X. (2021). Self-induced stochastic resonance in an excitable potential well. *Phys. Lett. A* 410, 127520. doi: 10.1016/j.physleta.2021.127520

Zhu, J., Kato, Y., and Nakao, H. (2022). Phase dynamics of noise-induced coherent oscillations in excitable systems. *Phys. Rev. Res*. 4, L022041. doi: 10.1103/PhysRevResearch.4.L022041

Zhu, J., and Liu, X. (2017). Delay-induced locking in bursting neuronal networks. *Chaos* 27, 083114. doi: 10.1063/1.4998927

Zhu, J., and Nakao, H. (2021). Stochastic periodic orbits in fast-slow systems with self-induced stochastic resonance. *Phys. Rev. Res*. 3, 033070. doi: 10.1103/PhysRevResearch.3.033070

## Appendix

## Noise-induced and noise-tuned bursting in excitable Hedgehog burster

Although the Hedgehog burster investigated in this paper has a stable limit cycle without noise, we have shown that the stochastic periodic orbits induced by SISR are different from the deterministic limit cycle. In fact, even when there is no deterministic limit cycle in the Hedgehog burster, i.e., when the bifurcation parameter *a* > 1 in system (Equation 1), noise can still induce coherent oscillations.

Different from the oscillatory case, the excitable Hedgehog burster has a lower bound of the noise strength (Yamakou and Jost, 2018; Zhu and Nakao, 2021), below which the transition to the right branch approximately obeys a Poisson process. The lower bound of the noise strength is given by the point at which the intersection of the lhs and rhs ceases to exist. The critical values are plotted in the inset of Figure 8F. As expected, the lower bound increases with *a* since an earlier transition requires a larger noise strength. Figures 8A–E illustrate stochastic trajectories above and below the lower bound of the noise strengths, where coherent oscillations can be induced in the former case while the noise is too weak to initiate transitions in the latter case.

The noise-tuned bursting behaviors can be similarly investigated as in the main text. The difference is that there is a lower bound of the noise strength in the excitable case in contrast to the oscillatory case and coherent oscillations can only be observed above this bound. As a result, the spike counts in each burst can be also influenced by the bifurcation parameter. For example, there are at most five spikes for *a* = 1.3 (Figure 8C) whatever the noise strength is. Thus, in the excitable Hedgehog burster, the bifurcation parameter *a* can also have significant influences on the noise-tuned bursting. On the one hand, the transition position on the left branch can be raised by lifting the equilibrium point. On the other hand, the lower bound of the noise strength for coherent oscillations will also modify the transitions on the right branch.

Keywords: self-induced stochastic resonance, noise-induced phenomenon, bursting, mixed-mode oscillation, critical transition, fast-slow system

Citation: Zhu J and Nakao H (2022) Noise-tuned bursting in a Hedgehog burster. *Front. Comput. Neurosci.* 16:970643. doi: 10.3389/fncom.2022.970643

Received: 16 June 2022; Accepted: 11 July 2022;

Published: 28 July 2022.

Edited by:

Germán Mato, Bariloche Atomic Centre (CNEA), ArgentinaReviewed by:

Marius E. Yamakou, University of Erlangen Nuremberg, GermanyEmre Baspinar, Université Paris Saclay, France

Copyright © 2022 Zhu and Nakao. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Jinjie Zhu, zhu.j.ag@m.titech.ac.jp