Delay-induced self-oscillation excitation in the FitzHugh-Nagumo model: regular and chaotic dynamics

TL;DR

This paper explores how time-delayed feedback in the FitzHugh-Nagumo model induces both regular and chaotic self-oscillations, expanding understanding beyond noise control to deterministic chaos via numerical and analytical methods.

Contribution

It demonstrates that delayed feedback can generate chaotic dynamics and chaos routes in the FitzHugh-Nagumo model, not just control coherence resonance.

Findings

Delayed feedback induces regular and chaotic oscillations.

Chaos can be achieved via a Feigenbaum route.

Large delays produce pattern-like dynamics.

Abstract

The stochastic FitzHugh-Nagumo model with time delayed-feedback is often studied in excitable regime to demonstrate the time-delayed control of coherence resonance. Here, we show that the impact of time-delayed feedback in the FitzHugh-Nagumo neuron is not limited by control of noise-induced oscillation regularity (coherence), but also results in excitation of the regular and chaotic self-oscillatory dynamics in the deterministic model. We demonstrate this numerically by means of simulations, linear stability analysis, the study of Lyapunov exponents and basins of attraction for both positive and negative delayed-feedback strengths. It has been established that one can implement a route to chaos in the explored model, where the intrinsic peculiarities of the Feigenbaum scenario are exhibited. For large time delay, we complement the study of temporal evolution by the interpretation of the dynamics as patterns in virtual space.

Tables

Table 1. Table 1: Feigenbaum relation values ( 4 ) characterizing transition to chaotic regime in model ( 1 ) with negative feedback strength.

period, $2^n$	${\gamma }_n$	relations ${\delta }_n$	comparing with $\delta$
$2^1=2$	-0.1878	4.335	7.16%
$2^2=4$	-0.3846	2.64	43.47%
$2^3=8$	-0.43	4.3	7.91%
$2^4=16$	-0.4472	–	–
$2^5=32$	-0.4512	–	–

Equations

$$\begin{array}[]{l}\varepsilon\dot{x}=x-x^{3}/3-y+\gamma(x_{\tau}-x),\\ \ \dot{y}=x+a,\end{array}$$

$$\Lambda =\overline{\underset{T\to \infty }{\lim }}\frac{1}{T}\ln \frac{||\mathbf{\xi }(t=T)||}{||\mathbf{\xi }(t=0)||},$$

$$\left\{\begin{array}[]{l}\varepsilon\dot{\xi_{1}}=\xi_{1}-x^{2}\xi_{1}-\xi_{2}+\gamma({\xi_{1}}_{\tau}-\xi_{1}),\\ \ \dot{\xi_{2}}=\xi_{1},\end{array}\right.$$

$${\delta }_n=\frac{{\gamma }_n-{\gamma }_{n+1}}{{\gamma }_{n+1}-{\gamma }_{n+2}},$$

$$\det \left|\mathbf{J}+{\mathbf{J}}_{\tau }\cdot e^{-\lambda \tau }-\lambda \cdot \mathbf{I}\right|=0,$$

$$\mathbf{J}=\left[\begin{array}{cc}\frac{\partial f_x}{\partial x}& \frac{\partial f_x}{\partial y}\\ \frac{\partial f_y}{\partial x}& \frac{\partial f_y}{\partial y}\end{array}\right]=\left[\begin{array}{cc}(1-x_s^2-\gamma )/\epsilon & -1/\epsilon \\ 1& 0\end{array}\right],$$

$${\mathbf{J}}_{\tau }=\left[\begin{array}{cc}\frac{\partial f_x}{\partial x_{\tau }}& \frac{\partial f_x}{\partial y_{\tau }}\\ \frac{\partial f_y}{\partial x_{\tau }}& \frac{\partial f_y}{\partial y_{\tau }}\end{array}\right]=\left[\begin{array}{cc}\gamma /\epsilon & 0\\ 0& 0\end{array}\right].$$

$$\det \left[\begin{array}{cc}-\frac{1}{\epsilon }(b+\gamma )-\lambda +\frac{\gamma }{\epsilon }{e}^{-\lambda \tau }& -\frac{1}{\epsilon }\\ 1& -\lambda \end{array}\right]=0$$

$$\epsilon {\lambda }^2+\lambda (b+\gamma )-\gamma \lambda e^{-\lambda \tau }+1=0.$$

$$i\omega (b+\gamma )-\epsilon {\omega }^2-i\omega \gamma \cos \omega \tau -\omega \gamma \sin \omega \tau +1=0$$

$$\{\begin{array}{l}\cos \omega \tau =\frac{b+\gamma }{\gamma }\\ \sin \omega \tau =\frac{1-\epsilon {\omega }^2}{\omega \gamma }\end{array}$$

$$\begin{array}[]{c}\frac{(b+\gamma)^{2}}{\gamma^{2}}+\frac{(1-\varepsilon\omega^{2})^{2}}{\omega^{2}\gamma^{2}}=1,\\ \varepsilon^{2}\omega^{4}+\omega^{2}\big{(}b(b+2\gamma)-2\varepsilon\big{)}+1=0.\end{array}$$

$${\epsilon }^2{\omega }^4+\xi {\omega }^2+1=0$$

$$\xi \le -2\epsilon .$$

$$\left[\begin{array}[]{l}\omega^{2}_{1,2}=\frac{-\xi+\sqrt{D}}{2\varepsilon^{2}}\\ \omega^{2}_{3,4}=\frac{-\xi-\sqrt{D}}{2\varepsilon^{2}}\end{array}\right.$$

$$\left[\begin{array}[]{l}\begin{cases}2\varepsilon+\xi\geq 0\\ D\geq(2\varepsilon+\xi)^{2}\end{cases}\\ \begin{cases}2\varepsilon+\xi<0\\ D\geq 0\end{cases}\end{array}\right.\Leftrightarrow\left[\begin{array}[]{l}\xi=-2\varepsilon\\ \xi<-2\varepsilon\end{array}\right.\Leftrightarrow\xi\leq-2\varepsilon.$$

$$\{\begin{array}{l}D\le (\xi +2\epsilon {)}^2\\ D\ge 0\\ -\xi -2\epsilon \ge 0\end{array}\iff \{\begin{array}{l}\xi \ge -2\epsilon \\ \xi \le -2\epsilon \end{array}\iff \xi =-2\epsilon .$$

$$\begin{array}[]{c}\omega^{2}\geq 1/\varepsilon\Leftrightarrow\begin{cases}\gamma\leq-b/2\\ \left[\begin{array}[]{l}\omega=\omega_{1},\ \omega>0\\ \omega=\omega_{2},\ \omega<0\end{array}\right.\end{cases}\\ \omega^{2}\leq 1/\varepsilon\Leftrightarrow\begin{cases}\gamma\leq-b/2\\ \left[\begin{array}[]{l}\omega=\omega_{3},\ \omega>0\\ \omega=\omega_{4},\ \omega<0\end{array}\right.\end{cases}\end{array}$$

$$\begin{cases}\cos\omega\tau\geq 0\\ \sin\omega\tau\geq 0\end{cases}\Leftrightarrow\begin{cases}\gamma\leq-b\\ \left[\begin{array}[]{l}\begin{cases}\omega<0\\ \omega^{2}\leq 1/\varepsilon\end{cases}\\ \begin{cases}\omega>0\\ \omega^{2}\geq 1/\varepsilon\end{cases}\end{array}\right.\end{cases}$$

$$\begin{cases}\gamma\leq-b\\ \left[\begin{array}[]{l}\omega=\omega_{1}\\ \omega=\omega_{4}\end{array}\right.\end{cases}$$

$$\begin{cases}\cos\omega\tau\leq 0\\ \sin\omega\tau\geq 0\end{cases}\Leftrightarrow\begin{cases}-b\leq\gamma\leq-b/2\\ \left[\begin{array}[]{l}\omega=\omega_{1}\\ \omega=\omega_{4}\end{array}\right.\end{cases}$$

$$\begin{cases}\cos\omega\tau\leq 0\\ \sin\omega\tau\leq 0\end{cases}\Leftrightarrow\begin{cases}-b\leq\gamma\leq-b/2\\ \left[\begin{array}[]{l}\omega=\omega_{2}\\ \omega=\omega_{3}\end{array}\right.\end{cases}$$

$$\begin{cases}\cos\omega\tau\geq 0\\ \sin\omega\tau\leq 0\end{cases}\Leftrightarrow\begin{cases}\gamma\leq-b\\ \left[\begin{array}[]{l}\omega=\omega_{2}\\ \omega=\omega_{3}\end{array}\right.\end{cases}$$

$$\gamma =-\frac{b}{1-\cos \omega \tau }.$$

$$\begin{cases}\gamma\leq-b\\ \left[\begin{array}[]{l}\gamma=-\frac{b}{1-\cos\omega_{1}\tau},\ \ \ \mathrm{if}\ \ \omega_{1}\tau\in\mathrm{\ I^{st}\ quarter}\\ \gamma=-\frac{b}{1-\cos\omega_{4}\tau},\ \ \ \mathrm{if}\ \ \omega_{4}\tau\in\mathrm{\ I^{st}\ quarter}\end{array}\right.\end{cases}$$

$$\begin{cases}-b\leq\gamma\leq-b/2\\ \left[\begin{array}[]{l}\gamma=-\frac{b}{1-\cos\omega_{1}\tau},\ \ \ \mathrm{if}\ \ \omega_{1}\tau\in\mathrm{\ II^{nd}\ quarter}\\ \gamma=-\frac{b}{1-\cos\omega_{4}\tau},\ \ \ \mathrm{if}\ \ \omega_{4}\tau\in\mathrm{\ II^{nd}\ quarter}\end{array}\right.\end{cases}$$

Full text

††journal: Chaos, Solitons & Fractals

\affiliation

[inst1]organization=Institute of Physics, Saratov State University,addressline=83 Astrakhanskaya str., city=Saratov, postcode=410012, country=Russia

Introduction

Time delay is inherently present in a wide range of systems in physics, biology, economics, and other fields [1]. In particular, a manifold of biological processes associated with time delay ranges from individual [2] and collective [3] neural cell activity to the population dynamics [4, 5]. A large variety of applications associated with natural delay and imposed time-delay processes includes medical issues such as therapy for Parkinson’s disease management [6, 7, 8] and control of other pathological brain rhythms (for instance, epileptic seizures) [9, 10], cardiac rhythm control [11], epidemic propagation control [12, 13]. In addition, time-delay processes play a significant role in chemical reaction rates [14, 15], laser dynamics [16, 17, 18, 19], analog and digital electronic circuits [20, 21], traffic flow simulation [22], climate change modelling [23].

Besides the presence of delays in biological neural networks and their influence on the neural network behaviour [24, 25, 26, 3, 27], we would like to emphasize the impact of time-delay loops in the context of artificial intelligence. Time delays are widely used for the development of delay-based machine learning algorithms involving spiking neural networks [28, 29]. Long time-delay loops provide for physical realization of reservoir computing architecture [30, 31]. The introduction of time-delay reservoir computing enabled simple optical, electronic and opto-electronic hardware implementations, which led to improvements of computation time scales for supervised learning [32, 33, 34, 35]. Moreover, deep-learning architectures have been adopted for the reservoir computing [36, 37, 38]. The delay-based reservoirs were successfully applied to a wide range of tasks, such as chaotic time series forecasting or speech recognition. In addition, using a single bistable delayed-feedback dynamical node as a network, one can synthesize a physical device (for instance, see the delayed-feedback-based Ising machine [39]) for solving problems with non-deterministic polynomial-time hardness (so-called NP-hard problems).

In the context of nonlinear dynamics, time delay represents a factor being responsible for a broad spectrum of fundamental phenomena such as delay-induced regular oscillatory behavior [40, 41, 42, 43], delay-induced bifurcations [44, 45] and multistability [46, 47], stabilization of steady states and periodic orbits [48, 49], vibrational resonance [50, 51], delay-induced chaos [52, 53], delayed-feedback control of chaos [48] (including spatiotemporal patterns [54, 55]), delay-induced synchronization [56, 57, 58] and desynchronization [59, 60]. In addition to the deterministic effects, time delay represents a useful approach for controlling stochastic phenomena such as stochastic [61, 62, 63] and coherence [64, 65, 66, 67, 68, 69, 70] resonances, noise-induced chimera states in networks of nonlocally-coupled excitable oscillators [71] and noise-induced patterns in excitable media [72]. Thus, delayed feedback is a powerful tool for achieving a wide range of operating regimes, enhancing amplitude-frequency characteristics, and controlling the dynamics of nonlinear systems [73, 74, 75].

As was first mentioned in Ref. [76], there exists an analogy between the behavior of time-delayed systems and the dynamics of ensembles of coupled oscillators or spatially extended systems [77]. The similarity takes place when the delay time, $\tau$, is much longer than system’s response time, which allows the system to reveal spatio-temporal phenomena (for example, coarsening [78], chimera states [79, 80, 81], soliton dynamics [18, 21, 82]) in the purely temporal dynamics of a time-delay system. This space-time analogy can be obtained by implementing a space-time transformation of the delay-feedback system, where the temporal dynamics is mapped onto space-time ($\sigma$, $n$) by introducing the space-time map $t=n\tau+\sigma$ with an integer (slow) time variable $n$, and a pseudo-space variable $\sigma\in[0,\eta]$, where $\eta=\tau+\delta$ with a small quantity $\delta$. For each set of parameters a unique value $\eta$ can be chosen such that the oscillatory dynamics is periodic with the period $\eta$. The mentioned technique is described in details in review [83] as well as a manifold of spatio-temporal phenomena tracked down in the dynamics of delayed-feedback oscillators. In addition, the large delay time in one delay-feedback system can be interpreted as one-directional impact of the identical system on the considered oscillator without any feedback [71, 84]. In this case one may expect appearance of self-sustained oscillations.

It is known that delayed feedback provides for the realization of self-oscillatory dynamics being unachievable in the absence of delay. In particular, it has been established that delayed feedback loops allow to shift critical parameter values corresponding to the subcritical Andronov-Hopf bifurcation. This fact became a base for the delayed-feedback-based coherence resonance control scheme in non-excitable systems [69, 70]. In contrast, the mechanism of noise-induced oscillation control by using delayed-feedback in the FitzHugh-Nagumo system (which is an excitable oscillator) in the regime of coherence resonance is not associated with changing the critical parameter value and the occurrence of the supercritical Andronov-Hopf bifurcation [66, 64]. The delay-induced self-oscillation excitation in the FitzHugh-Nagumo neuron has been studied in the context of deterministic dynamics for the short-time-delay limit [85] as well as for noticeable time delay [86, 87]. In all the cases the time delay feedback induces bistability as the coexistence of the regular self-oscillatory dynamics (spiking) and a quiescent steady state regime, whereas the stable equilibrium is a single attractor at zero delay.

In the current study, we summarize the materials published in Refs. [85, 86, 87] and complement them by our results which describe the FitzHugh-Nagumo oscillator with time-delayed feedback for both negative and positive delayed-feedback strength. Based on the linear stability analysis of steady states and numerical modelling of the equations under study, we demonstrate that the action of delayed-feedback loop can result in oscillations arising in the FitzHugh-Nagumo model accompanied with loss of stability of equilibrium point. In addition, we analyse basins of attractions, the probabilities to observe the self-oscillations, and the evolution of Lyapunov exponents. Besides the delay being smaller or comparable with the system response time, we also consider the case of long delay. By this way, we combine the study in the context of temporal evolution with the research of structures in quasi space.

1 System under study

The FitzHugh-Nagumo system originally proposed for the description of processes in nerve membranes [88, 89] is a paradigmatic model for the type-II excitability. In the present paper, we consider the FitzHugh-Nagumo model with time-delayed feedback:

[TABLE]

where $x=x(t)$ is the activator variable, $y=y(t)$ is the inhibitor variable. A parameter $\varepsilon\ll 1$ is responsible for the time scale separation of fast activator and slow inhibitor variable, $a$ is the threshold parameter which determines the system dynamics: the system exhibits the excitable regime at $|a|>1$ and the oscillatory one for $|a|<1$. A parameter $\tau$ is the delay time, $x_{\tau}$ is a value of the variable $x$ at the time moment ($t-\tau$). A parameter $\gamma$ is the feedback strength. The feedback term $\gamma(x_{\tau}-x)$ was first introduced by K. Pyragas to stabilize unstable periodic orbits resulting in the stabilization of the chaotic dynamics [48], but is also used for controlling various kinds of the deterministic and stochastic dynamics (see the references in the introduction). In this paper, we consider the FitzHugh-Nagumo system in the excitable regime ($a>1$ and $\varepsilon>0$) for varying delay time and time-delayed feedback strength.

Besides the analytical approaches described in appendices, our results are based on numerical simulations carried out by the integration of the model equation by using the Heun method [90] with time step $h=0.001$ starting from random or specially prepared initial conditions ($x_{0}$,$y_{0}$), where $x_{0}=x(t\in[-\tau:0])=\mathrm{const}$.

Using the term ’initial conditions’ ($x_{0}$,$y_{0}$) we mean the following. When integrating Eqs. (1) during the numerical modelling, a value of variable $x$ at the time moment $(t+h)$ is determined by the previous values of the dynamical variables, $x(t)$ and $y(t)$, as well as by a value $x(t-\tau)$ (here, $h$ is the integration time step). This means that the integration of system (1) starting from the initial time moment $t=0$ requires the initialization of $\tau/h$ values of $x$ in the time range $[-\tau:0)$. We set all of them to be identical and equal to $x_{0}$. In addition, we initialize a single initial value of the $y$-variable, $y_{0}=y(t=0)$ for starting the numerical simulations.

2 Positive time-delayed feedback strength ($\gamma>0$)

2.1 Temporal dynamics

As noted above, we consider system (1) in the excitable regime, occurring when the threshold parameter $a$ is greater than unity. In this case the system (1) without time-delayed feedback is characterised by the presence of a single attractor in the phase space: a stable focus with coordinates $x_{s}=-a$, $y_{s}=\frac{a^{3}}{3}-a$. However, in the presence of time delay, one can observe trajectories tracing self-oscillations which are not associated with the initial single attractor. The possibility to realize the self-oscillatory dynamics depends both on the delayed-feedback parameters ($\gamma$ and $\tau$) and the initial conditions. In order to quantitatively characterize this fact, we introduce the probability $P$ to achieve the regimes where the phase point oscillates outside the vicinity of the equilibrium point $(x_{s},y_{s})$. The corresponding probability is given in Fig. 1(a) for varying feedback strength $\gamma$ and delay time $\tau$ and fixed parameters $a=1.001$ and $\varepsilon=0.05$. Here, the probability for each $\gamma$ and $\tau$ is calculated by modelling the system 1500 times starting from the grid initial conditions $x_{0}\in[-2.5;2.5]$ and $y_{0}\in[-1.5;1.5]$, where $x(t\in[-\tau:0])=x_{0}=\mathrm{const}$. To eliminate the action of transients, the integration time is sufficiently long, $t_{\text{total}}=10000$. After each integration cycle, the position of the phase point was analyzed for proximity to the equilibrium $(x_{s},y_{s})$. Finally, a number of trajectories outside the vicinity of equilibrium was divided by the total number of initial conditions. This method for calculating probability is called ’phase space probability’ in some literature, where a set of virtual copies of considered system starting from different initial conditions is called a statistical ensemble, and averages computed with this definition of probability are called ensemble averages [91].

Figure 1(a) shows a large area corresponding to $P=0$ and coloured in black. In this area, we observe the only quiescent steady state regime associated with the equilibrium $(x_{s},y_{s})$. Outside this area the probability takes values in the range $P\in(0:1)$, which indicates the coexistence of the steady state regime with undamped oscillations.

The multistability of a system without time-delayed feedback can be visualized by attractors in a phase space and their basins of attraction, showing a set of initial conditions leading to corresponding attractor. The system under consideration is of infinite dimension, and therefore we can show only a section of the basins of attraction in the functional space and projection of phase trajectories on the phase plane ($x,y$). Figure 2 shows these results for different $\tau$ values. The section of the basins of attraction was obtained in the next way. The selected ranges of ($x,y$) values were divided into a lattice of initial conditions ($x_{0},y_{0}$). It is described in the last paragraph of Sect. 1 how the initial conditions are specified. The initial conditions leading to stable steady state shown by orange cross are colored in yellow, while the initial conditions leading to the blue attracting set are colored in light blue.

In particular, Fig. 2 illustrates the evolution of attracting sets by projections on phase plane ($x,y$) and corresponing section of basins of attraction for fixed feedback strength, $\gamma=0.025$ (this level is indicated by the dashed line in Fig. 1 (a)) and increasing time delay starting from the value $\tau=3.08$ (marked by the symbol ’$\ast$’ in Fig. 1 (a)). The coexistence of two attractors is illustrated in Fig. 2. The first attractor, the stable steady state (the orange cross) is characterised by yellow basin of attraction. The second attractor, a stable limit cycle (the blue curve), is characterised by the basin coloured in light-blue. As demonstrated in Fig. 2, the section of basins of attraction evolve with increasing time delay.

The maximal size of the limit cycle’s basin can be found in the center of area I in Fig. 1(a). The corresponding section of basin of attraction is shown in Fig. 2(c). The further increase of the parameter $\tau$ inside part I leads to decrease of the basin (Fig. 2(d)) and its total disappearance when $\tau$ comes back again to the black area of Fig. 1(a). The same birth of oscillations and their disappearance can be observed in the next four parts: II, III, IV, V.

In addition to the visual distinction between areas I–V in Fig. 1(a), they differ in the type of oscillations. The spiking behaviour with period-1 limit cycle can be observed in the part I of blue area in Fig. 1(a). The corresponding time realization is given in Fig. 3(a). With further increasing the delay time $\tau$ and moving further into the parts II, III, IV, the number of periods is growing as period-2 oscillations (Fig. 3(b)) in Part II, period-3 oscillations (c) in Part III, and period-4 oscillations (Fig. 3(d)) in Part IV. A similar pattern persists with a further increase in the parameter $\gamma$.

The existence of observed oscillatory regimes is confirmed by calculation of maximal Lyapunov exponents (Fig. 1(b)) on the same parameter plane ($\tau$,$\gamma$). The map of Lyapunov exponents is prepared in the next way. The main Lyapunov exponent is determined as follows

[TABLE]

where the initial length $||\boldsymbol{\xi}(t=0)||$ is set to be equal to unity, while the final length $||\boldsymbol{\xi}(t=T)||$ is obtained by simulating the system (1) together with the following system of equations:

[TABLE]

where $\xi_{1},\xi_{2}$ are the components of the vector $\boldsymbol{\xi}$, and parameters $\varepsilon,\gamma,\tau,x,y$ have the same meaning as in Eq. (1).

The black areas in Fig. 1(b) with negative maximal Lyapunov exponents $\Lambda$ correspond to the existence of only stable equilibrium point without any oscillations. The blue region with $\Lambda=0$ confirms the existence of periodic oscillations for chosen parameters $\gamma$ and $\tau$.

We have calculated the Lyapunov exponents in a wider range of the feedback strength $\gamma\in[0;1]$ and have found that the previously obtained trend remains. Only two regimes were observed in the considered parameter plane including stable equilibrium point with $\Lambda<0$ for small $\gamma$ and $\tau$ and its coexistence with periodic oscillations with $\Lambda=0$ for the rest of the parameter values. Regimes with positive $\Lambda$ have not been found for $\gamma>0$, which confirms that the delay-induced self-oscillations are regular. These calculations are given in A.

2.2 Soliton structures at large time delay

If the delay time is sufficiently long, one can complement studying model (1) in the frameworks of the temporal evolution by the consideration of the system in the context of spatio-temporal structures. In particular, the mentioned approach applied to Eqs. (1) with positive delay feedback strength allows to reveal soliton structures [87]. It is important to note that such structures are observed in the nonlinear dissipative system, can persist in the presence of noise, save their shape and do not decay over time. For these reasons, they can be referred to the dissipative solitons. Choosing the initial conditions or using a temporal external influence, one can induce single [Fig. 4(a)] and multiple [Fig. 4(b)-(c)] dissipative solitons in system (1). The oscillatory regimes illustrated in Fig. 4(a)-(c) represent a particular kind of delay-induced self-oscillatory dynamics: the observed oscillations are obtained at $a=1.01$ where the self-oscillations cannot be exhibited in the absence of time delay. On projections on the phase plane, the soliton structures in Fig. 4(a)-(c) consist of self-sustained alternating spike loops and the rest states at the equilibrium [Fig. 4(d)].

The multiplicity of the dissipative solitons can be potentially used in the context of information storage. Indeed, adjusting the initial conditions, any number of solitons can be disposed in the virtual space and vice versa, each assigned soliton can be canceled by a slight external impact. However the minimal distance between the units is restricted. This is a limitation which bounds the maximal possible number of solitons in space. Nevertheless, the mentioned restriction can be overcame by increasing time delay. Then the information capacity of the system can grow.

3 Negative time-delayed feedback strength ($\gamma<0$)

Let us consider the regimes induced by the time-delayed feedback with negative strength $\gamma$. As the first step, we use the same technique as in Sec. 2 to find out any exit of trajectory outside the vicinity of the equilibrium point ($x_{s}$,$y_{s}$) using the probability analysis. The corresponding probability on the plane ($\gamma$,$\tau$) is shown in Fig. 5 for two parameter values $a=1.001$ panel (a)) and $a=1.01$ (panel (b)). There is the same color scheme as in Sec. 2, but comparing Figs. 1(a) and 5, one can see a clearly pronounced difference in maximum probability values. In contrast to probabilities calculated for $\gamma>0$ in Sect. 2, now for $\gamma<0$, there are clear parameter areas with $P\to 1$ and $P\to 0$. It means, that there is no bistability for $\gamma<0$. In the black area of parameter values, the system demonstrates a stable equilibrium point while the parameters from the white area lead to a completely different dynamics. The probability was prepared again with global approach with listing possible initial conditions $x_{0}\in[-2.5;2.5]$ and $y_{0}\in[-1.5;1.5]$, and, therefore, the probability $P=1$ means that no trajectories come to the vicinity of the equilibrium point.

In order to uncover the mechanism of birth of oscillations caused by the feedback with negative strength we fix the parameter $\gamma=-0.03$ (dashed line in Fig. 5(b)) and move from left to right increasing the delay time $\tau$. First, in the black area for $\tau\approx 0$ only the stable equilibrium point is observed. Next, in the beginning of the white area marked by the symbol $\ast$ in Fig. 5(b), a limit cycle of period-1 with a small amplitude is born around the previously stable equilibrium point. The corresponding temporal realization with corresponding attracting set in projection on phase plane ($x,y$) is given in Fig. 6(a). The amplitude of the obtained limit cycle grows with distancing from purple area. Therefore, the appearance of small yellow areas in Fig. 5 is caused not by the multistability, but by the smallness of the limit cycle which is comparable with the size of the vicinity used for probability calculations. These oscillations have large decay, and finally end with equilibrium point.

At some point, the period doubling bifurcation takes place, and a new loop of the limit cycle appears (see Fig. 6(b) for $\tau=0.45$). In the center of the white area situated between symbols $\ast$ and $\star$ the period-3 cycle can be obtained (see Fig. 6(c) for $\tau=0.55$). This regime combines oscillations near the previous fixed point and the spiking behaviour. With increasing $\tau$ and getting close to $\star$ from the left, the number of loops decreases and period-2 can be obtained (Fig. 6(d), $\tau=0.7$) with final disappearance of any oscillations when $\tau$ gets to the black area again.

With further increase of $\tau>1.5$ the system demonstrates a wide range of regimes: cycles with different period numbers and spiking behavior, like in Fig. 6(a–d). The spiking behavior is the most typical regime in the system with negative time-delayed feedback.

Such loss of stability of the equilibrium point is proven by the linear stability analysis (see B).

The obtained loss of stability has also been proven by the Lyapunov exponents (Fig. 7) described first in Sect. 2. The existence of the stable fixed point in the black area is characterised by $\Lambda<0$, while the rest part (blue) of parameter plane $\tau\in[0;10]$, $\gamma\in[-0.1;0]$ has $\Lambda=0$, corresponding to periodic oscillations. At the same time, if we enlarge the feedback strength in absolute value, there are additional interesting regimes accompanied by $\Lambda>0$, which indicated the chaotic behaviour (yellow regions in Fig. 7).

3.1 Large feedback strength. Chaotic oscillations

The largest part of Fig. 7(a) is colored in blue corresponding to periodic oscillations with $\Lambda\approx 0$. The main difference of Lyapunov exponents calculated for positive and negative $\gamma$ (Fig. 1(a) and Fig. 7) is that the parameter plane contains yellow regions with $\Lambda>0$ for negative $\gamma$. This indicates that chaotic oscillations can be observed for corresponding parameters $\gamma$ and $\tau$ from that regions. At the selected scale of parameter $\gamma$ (Fig. 7(a)), this covers a very small parameter area. However, when the absolute value of parameter $\gamma$ is increased (Fig. 7(b)), the number of such areas and their size grows. In this subsection, we take a closer look at these chaotic regimes.

The exact parameter values for which chaotic regimes can be observed, do not change significantly for $a=1.01$ and $a=1.001$. The similar calculations for $a=1.001$ are given in C, Fig. 14. In both $\gamma$ scales of Fig. 7, there are very tiny areas with $\Lambda>0$ for $\tau<4$. The corresponding representative regimes are illustrated by projection of attracting sets on phase plane ($x,y$) and temporal dynamics in Fig. 8 (a,b) for two values of delay time $\tau=0.05$ and $\tau=0.6$. In both cases, the system demonstrates chaotic oscillations of a small amplitude near the previously stable focus equilibrium point. Increasing the delay time $\tau$ one can observe a growth in the amplitude of chaotic oscillations (Fig. 8(c,d)). This is accompanied with enlarging the attractor in $(x,y)$ phase plane.

It is shown above that forward increase in $\tau$ value may transform the dynamics to the spiking regime (Fig. 6). The phase portrait and temporal evolution in Fig. 8(c,d) show that system (1) may demonstrate chaotic spiking behaviour. Positive value of the maximal Lyapunov exponent $\Lambda\approx 0.04$ confirms this fact.

In order to examine the mechanism of the chaos appearance we now consider the behaviour of the system (1) during transition from regular periodic regime at $\tau=7$, $\gamma=-0.15$, $a=1.01$ to the largest region of existing chaotic regimes for $\tau\in(7;9)$ and $\gamma\in(-1;-0.4)$ shown in Fig. 7(b). The projections of attracting sets on phase plane ($x,y$) and temporal dynamics obtained during this transition are shown in Fig. 9 and demonstrate that the transition to the chaotic behaviour appears through the cascade of period-doubling bifurcations. The dependencies of the first and the second Lyapunov exponents on the parameter $\gamma$ are given in Fig. 10. The second Lyapunov exponent reaches zero value every time when the period of cycle is doubled. The Feigenbaum relations

[TABLE]

calculated in Tab. 1 are relatively close (with error about 7% for $\delta_{1}$ and $\delta_{3}$) to value $\delta=4.669$ which correspond to $\lim_{n\to\infty}\delta_{n}$ in one-parameter maps. Finally, at $\gamma_{tr}=-0.45$ the system (1) comes to the chaotic regime, has been previously shown in Fig. 8(c,d).

3.2 Regular and chaotic waves at long delay

For negative $\gamma$ and long time delay the represented in quasi-space dynamics of system (1) is characterised by the excitation of wave-like patterns where the spatial impulses homogeneously fill the quasi space. Such structures do not involve the quiescent steady state regimes. This indicates that the delayed-feedback represents a factor for the steady state’s instability at negative values $\gamma$. In contrast to the soliton structures obtained for positive $\gamma$, here one cannot induce a certain number of impulses by varying the initial conditions. The coexistence of chaotic [Fig. 11 (a,b)] and regular [Fig. 11 (c,d)] self-oscillatory dynamics is achieved for big enough absolute values of the delay feedback strength.

4 Conclusion and discussion

Delayed-feedback provides for the self-sustained oscillation excitation in the FitzHugh-Nagumo model. In our paper, we demonstrate that the delay-induced self-oscillatory dynamics can be regular and chaotic. In particular, at positive delayed-feedback strength one observes the bistability: the quiescent steady state regime coexists with regular self-oscillatory motions along the limit cycle. Despite the presented analysis of basins of attraction is reduced and involves the special kind of initial conditions $x_{0}=x(t\in[-\tau:0])=$const, it allows to establish that both the steady state and the limit cycles are characterised by finite basins of attraction. Thus, they can be called ’attractors’ in a full sense.

For long time delay and positive delay feedback strength one can interpret the induced self-oscillations as dissipative solitons in the quasi-space. In such a case, the spatio-temporal dynamics involves both attractors and the delay-induced bistability is reflected in the ability to realize a certain number of solitons at certain times.

For positive delayed-feedback strength the behaviour is more complicated: after the steady state loses its stability, both regular and chaotic self-oscillations can be achieved. For long time delay, regular and chaotic self-oscillations are manifested in the quasi-space as wave patterns. Surprisingly, increasing the absolute value of the negative delayed-feedback strength, one can realize the route to chaos through a cascade of period-doubling bifurcations such that the Feigenbaum constants registered in numerical experiments are close to the theoretical value $\delta=4.669$. The exhibition of the Feigenbaum scenario by the delayed-feedback FitzHugh-Nagumo model seems to be a non-trivial results and will be theoretically analysed in further studies.

The obtained results are in a full correspondence with materials presented in Refs. [85, 86, 87] and other publications and complement a manifold of observed regimes by the chaotic dynamics.

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Acknowledgements

The main idea of the paper and the parts with probabilities, basins of attraction and linear stability analysis were prepared by N. Semenova, supported by The Council for grants of President of Russian Federation (project no. SP-749.2022.5). The calculations of Lyapunov exponents with finding the chaotic regimes were carried out by A. Bukh, supported provided by The Council for grants of President of Russian Federation (project no. SP-774.2022.5). Sects. 2.2, 3.2 and the main work with the text was done by V. Semenov, supported by the Russian Science Foundation (project no. 22-72-00038).

Appendix A Lyapunov exponents for positive time-delayed feedback strength $\gamma>0$

Appendix B Linear stability analysis for negative time-delayed feedback strength $\gamma<0$

The FHN system in excitable regime ($a>1$) without time delay demonstrates the stable equilibrium point of type “focus”. With decreasing parameter $a$ passing through $a=1$ this equilibrium point looses stability and becomes unstable. This can be illustrated by the roots of characteristic equation. The stable focus is accompanied by the complex roots $\lambda_{12}$ with negative real parts $\real(\lambda_{12})<0$ in excitable regime ($a>1$), and with positive $\real(\lambda_{12})>0$ in oscillatory regime ($0<a<1$). The bifurcation between these two regimes is accompanied with the case when the roots are purely imaginary with $\real(\lambda_{12})=0$. This equilibrium type is called “center” [92].

In the case of delayed system (1), the similar analysis can be applied. As it has been mentioned in the book [93], in the general case the characteristic equation for the system with time-delayed feedback can be obtained as follows:

[TABLE]

where $\mathbf{J}$ is the Jacobian matrix for considered system without time-delayed terms, $\mathbf{J_{\tau}}$ is the Jacobian matrix for time-delayed terms, $\mathbf{I}$ is the unity matrix, $\lambda$ are the roots of characteristic equation, which type determines the type of equilibria. In the case of our system (1) these matrices are:

[TABLE]

where $f_{x}$ and $f_{y}$ are the right and left parts of Eq. (1), $x_{s}$ is the $x$-value of the equilibrium state equal to $x_{s}=-a$. The Jacobian matrix for delayed part is:

[TABLE]

Let us denote $b=a^{2}-1$, which is positive when $a>1$. Then the characteristic equation is:

[TABLE]

[TABLE]

The purely imaginary $\lambda$ corresponds to the equilibrium of type “center”. Therefore, we set $\lambda=i\omega$, $\omega\in\mathbb{R}$, $\omega\neq 0$. Then $e^{-\lambda\tau}$ can be rewritten as $(\cos\omega\tau-i\cdot\sin\omega\tau)$ and (9) transforms to

[TABLE]

Dividing this equation into a system of two equations for the real and imaginary parts, we obtain

[TABLE]

The first equation can be rewritten in the form $\gamma=-b/(1-\cos\omega\tau)$, indicating that $\gamma<0$ when $a>1$. Therefore, next we will consider only the case of negative $\gamma$.

To find a solution $\omega$ to the system (10), we use the rule $\sin^{2}\omega\tau+\cos^{2}\omega\tau=1$:

[TABLE]

Denote $\xi=b(b+2\gamma)-2\varepsilon$. The parameters, which are used in this paper, can lead only to $\xi<0$. Then the equation trainsforms to

[TABLE]

with the discriminant $D=\xi^{2}-4\varepsilon^{2}$. The solution exists only if $D\geq 0$. And this leads to additional conditions for $\xi$: $\xi\geq 2\varepsilon$ or $\xi\leq-2\varepsilon$. Combining this with $\xi<0$, we get the final condition:

[TABLE]

Equation (12) produces four solutions $\omega_{1},\ \omega_{3}>0$ and $\omega_{2},\ \omega_{4}>0$:

[TABLE]

All four solutions exist when condition (13) holds. However, in what follows we will need to compare $\omega^{2}$ with the value $1/\varepsilon$.

1) $\omega^{2}\geq 1/\varepsilon$:

Let us find the condition for existence $\omega_{1,2}$. The inequality $\frac{-\xi+\sqrt{D}}{2\varepsilon^{2}}\geq\frac{1}{\varepsilon}$ simply transforms to $\sqrt{D}\geq 2\varepsilon+\xi$. This leads to the system of inequalities:

[TABLE]

Similarly, one can find a condition for the existence of $\omega_{3,4}$ from $\frac{-\xi-\sqrt{D}}{2\varepsilon^{2}}\geq\frac{1}{\varepsilon}$:

[TABLE]

The case $\xi=-2\varepsilon$ corresponds to $D=0$ and $\omega_{1,3}=\sqrt{-\xi/(2\varepsilon^{2})}$, $\omega_{2,4}=-\sqrt{-\xi/(2\varepsilon^{2})}$. Therefore, the condition (16) has already been taken into account in condition (15).

Substituting all the renaming into (15), we get $\gamma\leq-b/2$. Thus the condition $\omega^{2}\geq 1/\varepsilon$ means that $\gamma\leq-b/2$ and then only $\omega_{1,2}$ exists.  

 

2) $\omega^{2}\leq 1/\varepsilon$:

Similarly to the previous calculations, we get that condition $\omega^{2}\leq 1/\varepsilon$ means that $\gamma\leq-b/2$ and then only $\omega_{3,4}$ exists. This is illustrated by the next system of inequalities:

[TABLE]

Let us return to the system of equations (10) and consider the solutions depending on which quarter of the unity circle the angle $\omega\tau$ belongs to.

I) The first quarter: $0+2\pi n\leq\omega\tau\leq\pi/2+2\pi n$, $n\in\mathbb{Z}$.

[TABLE]

Including (17), we get

[TABLE]

II) The second quarter: $\pi/2+2\pi n\leq\omega\tau\leq\pi+2\pi n$, $n\in\mathbb{Z}$.

[TABLE]

III) The third quarter: $\pi+2\pi n\leq\omega\tau\leq 3\pi/2+2\pi n$, $n\in\mathbb{Z}$.

[TABLE]

IV) The fourth quarter: $3\pi/2+2\pi n\leq\omega\tau\leq 2\pi(n+1)$, $n\in\mathbb{Z}$.

[TABLE]

One can simply get the equation for $\gamma$ from (10):

[TABLE]

The solution $\omega$ depends on the parameter $\gamma$, and the bifurcation line can be therefore plotted in the parameter plane ($\gamma,\tau$), if Eq. (23) is solved graphically. The function $\cos$ is even and it is enough to consider only $\omega\tau$-quarters with existing $\omega_{1}$ and $\omega_{4}$. Figure 13 contains the graphical solution of the next systems of equations and inequalities:

[TABLE]

and

[TABLE]

where $b=a^{2}-1$, $\xi=b(b+2\gamma)-2\varepsilon$,

$\omega_{1}=\sqrt{\frac{-\xi+\sqrt{\xi^{2}-4\varepsilon^{2}}}{2\varepsilon^{2}}}$,

$\omega_{4}=-\sqrt{\frac{-\xi-\sqrt{\xi^{2}-4\varepsilon^{2}}}{2\varepsilon^{2}}}$

The dashed line in Fig. 13 separates the area upper the bifurcation line with stable focus and the lower area with unstable focus. This line purely coinsides with the results of numerical simulation (Fig. 5).

Appendix C Lyapunov exponents for positive time-delayed feedback strength $\gamma<0$