# Non-Trivial Dynamics in the FizHugh-Rinzel Model and Non-Homogeneous Oscillatory-Excitable Reaction-Diffusions Systems

Benjamin Ambrosio, M A Aziz-Alaoui, Argha Mondal, Arnab Mondal, Sanjeev K Sharma, Ranjit Kumar Upadhyay

<figure><img src="image_1.png" /></figure>

## To cite this version:

Benjamin Ambrosio, M A Aziz-Alaoui, Argha Mondal, Arnab Mondal, Sanjeev K Sharma, et al. "Non-Trivial Dynamics in the FizHugh-Rinzel Model and Non-Homogeneous Oscillatory-Excitable Reaction-Diffusions Systems." Biology, 2023, 12 (7), pp.918. 10.3390/biology12070918. hal-04292144

HAL Id: hal-04292144

https://normandie-univ.hal.science/hal-04292144v1

Submitted on 22 Nov 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

Article

# Non-Trivial Dynamics in the FizHugh-Rinzel Model and Non-Homogeneous Oscillatory-Excitable Reaction-Diffusions Systems

Benjamin Ambrosio<sup>1,2,*</sup> ID, M. A. Aziz-Alaoui<sup>1 ID</sup>, Argha Mondal<sup>3,4</sup>, Arnab Mondal<sup>5</sup>, Sanjeev K. Sharma<sup>5 ID</sup> and Ranjit Kumar Upadhyay<sup>5 ID</sup>

**Simple Summary:** This article provides original numerical and mathematical insights about the FHR model and non-homogeneous FitzHugh-Nagumo reaction-diffusion systems.

**Abstract:** This article focuses on the qualitative analysis of complex dynamics arising in a few mathematical models in neuroscience context. We first discuss the dynamics arising in the three-dimensional FitzHugh-Rinzel (FHR) model and then illustrate those arising in a class of non-homogeneous FitzHugh-Nagumo (Nh-FHN) reaction-diffusion systems. FHR and Nh-FHN models can be used to generate relevant complex dynamics and wave-propagation phenomena in neuroscience context. Such complex dynamics include canards, mixed-mode oscillations (MMOs), Hopf-bifurcations and their spatially extended counterpart. Our article highlights original methods to characterize these complex dynamics and how they emerge in ordinary differential equations and spatially extended models.

**Keywords:** FitzHugh-Nagumo; FitzHugh-Rinzel model; fast-slow dynamics; bifurcation; canard; mixed mode oscillations; bursting oscillations; neuroscience; waves

## 1. Introduction

The goal of this article is to shed light on mechanisms related to the emergence of complex oscillations arising in a few dynamical systems derivated from the FitzHugh-Nagumo equations (FHN) [1,2]. A specific model of interest here is a non-homogeneous FitzHugh-Nagumo reaction-diffusion (RD)-type system (Nh-FHN) where a parameter is space-dependent. This model derivates from the classical FHN model, which is itself a simplified version of the Hodgkin-Huxley (HH) model [3]. The HH reaction-diffusion model is fundamental in the field of mathematical neuroscience because it stands for the first model generating action potentials through ionic fluxes dynamics and because it allows for their propagation along the axon. Before digging into a more detailed description of the model under consideration here, it is worth providing some elements of context about recent contributions that illustrate how dynamical systems have become a crucial theoretical tool to describe the time evolution of neuronal activity and dynamical brain states. In [4], the authors show how the Drosophila central complex, a region implicated in goal-directed navigation, performs vector arithmetic. The authors describe a neural signal in the fan-shaped body that tracks the allocentric traveling angle of a fly. Previous work had identified neurons in Drosophila that track the heading angle. On this basis,

<figure><img src="image_4.png" /></figure>

check for
updates

**Citation:** Ambrosio, B.; Aziz-Alaoui, M.A.; Mondal, A.; Mondal, A.; Sharma, S.K.; Upadhyay, R.K. Non-Trivial Dynamics in the FizHugh-Rinzel Model and Non-Homogeneous Oscillatory-Excitable Reaction-Diffusions Systems. *Biology* **2023**, *12*, 918. https://doi.org/10.3390/biology12070918

Academic Editor: Claudio R. Mirasso

Received: 4 May 2023

Revised: 16 June 2023

Accepted: 19 June 2023

Published: 27 June 2023

<figure><img src="image_5.png" /></figure>

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

<footnote n="1">1 UNIHAVRE, LMAH, FR-CNRS-3335, ISCN, Normandie University, 76600 Le Havre, France</footnote>

<footnote n="2">2 The Hudson School of Mathematics, New York, NY 10001, USA</footnote>

<footnote n="3">3 Department of Mathematics, Sidho-Kanho-Birsha University, Purulia 723104, India; arghamondalb1@gmail.com</footnote>

<footnote n="4">4 Department of Mathematical Sciences, University of Essex, Wivenhoe Park, Colchester CO4 3SQ, UK</footnote>

<footnote n="5">5 Department of Mathematics and Computing, Indian Institute of Technology (Indian School of Mines), Dhanbad 826004, India; arnab.17dr000410@am.ism.ac.in (A.M.)</footnote>

<footnote>* Correspondence: benjamin.ambrosio@univ-lehavre.fr</footnote>
they characterize a neuronal circuit that performs a coordinate transformation and vector
addition to compute the allocentric traveling direction. For modeling purposes, the authors
use a network of dynamical systems where traveling directions show up as bumps of
activity representing calcium concentrations. In [5,6], the authors propose a network model
of excitatory and inhibitory neurons to describe the activity of the visual cortex. Their
work highlights emergent properties such as Gamma rhythms. It also illustrates how the
network, doing its own mathematical computations, mimics the visual cortex response to
orientated visual stimuli. Again, they claim that the network dynamical system, and the
brain by extrapolation, computes mathematical convolutions to deliver cortical responses.
Finally, in [7], Volpert et. al. investigate the spatiotemporal dynamics of neural oscillations
observed in real EEG data acquired during a picture-naming task. Their theoretical models
rely on one hand on analytical solutions of non-autonomous Poisson equations and on the
other hand on optical flow patterns such as sources, sinks, spirals and saddles. Their work
contribute to a better understanding of the neural dynamics at the macroscopic scale with
the objective to characterize neural activity corresponding to a specific behavior. Overall,
all these works consecrate the use of Dynamical Systems to characterize neural activity
in a specific context. They appeal for a better knowledge of the mechanisms inducing
oscillations, neural rhythms and spatiotemporal organization. This list of relevant works
could be extended widely. Here, our focus is on how complex oscillations emerge and
propagate in FitzHugh-Nagumo inspired systems. This question arises naturally in the
context of propagation of an action potential along an axon or in excitable media. In recent
decades, a significant progress relying on the properties of FHN type models has been
made in the theoretical description of cardiac tachycardia and fibrillation, see [8,9] and
references cited therein. For axonal propagation, the article [10] provides a typical example
of a phenomenological study relying on excitable models. In this article, the authors indeed
rely on a non-homogeneous FHN RD-type system to give numerical insights about how to
distinguish between geometrical enlargements that lead to minor changes in propagation
from those that result in critical phenomena such as blockages of the original traveling
spike. They focus specifically on the inhomogeneity in the diffusion's coefficient to take
into account the discrepancy of the axonal diameter. Our approach is by several aspects
related to the contribution of [10], but in a more theoretical framework. We seek to highlight
the detailed dynamical mechanisms that can shape complex oscillations and affect their
propagation across spatial domains. Coming back to the Nh-FHN model, the space domain,
typically stands for an excitable tissue such as a neural tissue or a neuron's axon. One
interesting property of the FHN ODE is that upon variation of a parameter, a globally
attracting limit-cycle emerges from an Hopf-bifurcation. Depending on the value of this
parameter, FHN can be either excitable or oscillatory. For the specific case of Nh-FHN, at
a designated restricted location of the spacial domain, the cells are assumed to be in an
oscillatory state, and elsewhere in the domain, the cells are excitable, see Figure 1. The
excitability of the cells depend on the value of the aforementioned parameter, which now
depends on the space variable, and which in application is related to the external current.
In this context, for Nh-FHN, we observe the following phenomenon; upon the variation of
this parameter value related to excitability, waves of depolarization may propagate from
the center of the domain toward the boundaries when the cells are enough excitable, or the
solution can evolve to a stationary state if the cells are not enough excitable. For a parameter
value range in between, a bifurcation occurs, typically a Hopf-bifurcation. Within this
range of parameters, for some cells, mixed mode oscillations or other complex dynamics
can be observed in time. This scenario was first described in [11] and further analyzed and
discussed in a few papers such as [12,13]. Other scenarios can be drawn with Nh-FHN, for
example one could include damaged tissue parts in the domain, see Figure 1. Numerical
simulations have been developed for this case, but the details of it are beyond the scope of
the present work. The wave propagation phenomenon in excitable media and in FHN RD
in particular has been intensively investigated for a few decades, see for example [14–21]
and references therein. Since the seminal study of the Fisher KPP equations, a general
theory of traveling waves has also been developed, see for example [22–24]. The numerical simulations to be presented in this article suggest, however, that the vast literature does not account for the details of the dynamics illustrated here. Our focus here is to present few oscillatory complex phenomena that appear to be tractable by analogy with lower dimensional systems. Our approach is to emphasize the analogy between some dynamics observed in the spatial Nh-FHN system for some cells and the dynamics of the FHR ODE system. This idea gives the structure of the present article. After the introductory part, we provide an original theoretical and numerical analysis of the FHR system in Section 2. We then illustrate numerical simulations and theoretical results for Nh-FHN in Section 3. Before digging into a more detailed analysis, we first introduce the equations to be considered, recalling also some necessary background.

<figure><img src="image_6.png" /><figcaption>**Figure 1.** Diagram illustrating the principle of Nh-FHN. FHN ODE models can be excitable or oscillatory depending on a parameter value. Adding a space variable and a diffusive term, one can consider domains where the parameter is space-dependent, opening the possibilities for very rich dynamics. In the left panel (a), we illustrate a 1D domain with an oscillatory part at the center and excitable tissue elsewhere. In the right panel (b), we illustrate a 2D domain with oscillatory parts, a damaged part and excitable tissue elsewhere. The damaged part can be drawn randomly with a parameter value ranging from values corresponding to very low excitability to high excitability or oscillatory dynamics.</figcaption></figure>

## 2. Methods

### 2.1. The FHR System

The FHR ODE model was introduced as a simple model to generate bursting oscillations, (i.e., alternated phases of oscillations and quiescent states) in [25,26], and previously studied by J. Rinzel and R. FitzHugh in an unpublished paper in 1979. In those fundamental papers, the FHR model reads as

$$ \left\{ \begin{aligned} \frac{dv}{dt} &= f(v) - w + I, \\ \frac{dw}{dt} &= \phi(a + v - bw), \\ \frac{dy}{dt} &= \epsilon(c - v - dy) \end{aligned} \right. \qquad (1) $$

with

$$ f(v) = -(1/3)v^3 + v $$

For the following values of parameters, $I = 0.3125$, $a = 0.7$, $b = 0.8$, $c = -0.775$, $d = 1.0$, $\phi = 0.08$, $\epsilon = 0.0001$ the model exhibits nice bursting behavior. The mechanism is strikingly simple and intuitive. The model consists of a FitzHugh–Nagumo (FHN) system, represented by the two first equations and a super slow variable, whose dynamics are given by the third equation. A relevant approach here is to formally consider $y + I$ as a parameter for the two first equations. Then, the dynamics of the two first equations are known from the FHN analysis. But the variable $y$ moves slowly with the variation of the first variable, allowing the emergence of the bursting. The description made in the
original papers [25,26] reflects a precise numerical qualitative analysis of the FHR dynamics supported by numerical illustrations. The author suggested that the model would be suitable for a deeper rigorous mathematical analysis, and invited the community to study the model in this direction. This call would receive important feedback, and the model was indeed studied later in numerous papers, see for example [27–30] to cite only but a few. More generally, the idea of adding one or more variables with different scales to capitalize on the known dynamics of FHN was successful in many aspects. For example, in [31], see also [32], the dynamics of FHN at different scales were used to model the alternating pulse and surge pattern of gonadotropin-releasing hormone secretions. In [33], a three-dimensional FHN served as an explanatory example to illustrate the inflection sets in the context of the excitability threshold, a method previously developed for planar systems, see [34] and references therein. In particular, such 3D systems are reference models for the emergence of complex oscillations such as MMOs and bursting, which have been widely observed in biological contexts [31,35–38]. The analysis of systems with multiple scales is a very active topic of research with numerous applications in biology, physics and chemistry, see [39–41]. Here, we will consider a slightly different version of the original FHR, for which we will provide an original numerical and theoretical analysis. The system under consideration writes:

$$ \left\{ \begin{aligned} \epsilon \frac{du}{dt} &= f(u) - v + w + I, \\ \frac{dv}{dt} &= u - bv - c \\ \frac{dw}{dt} &= \epsilon(-u - w) \end{aligned} \right. \qquad (2) $$

with f a cubic function, b > 0, ε a small parameter and I a parameter to vary. For the ODE section, we will provide an analysis with f and b set as

$$ f(u) = -(1/3)u^3 + u, \quad b = 0.8, \quad c = 0. $$

Next, we move forward with a short presentation of the Nh-FHN model.

## 2.2. A Non-Homogeneous FHN Model (Nh-FHN)

The non-homogeneous FHN reaction-diffusion model (Nh-FHN) to be considered here writes

$$ \left\{ \begin{aligned} \epsilon \frac{du}{dt} &= f(u) - v + w + I(x) + d_u \Delta u, \\ \frac{dv}{dt} &= u - bv - c(x) + d_v \Delta v \end{aligned} \right. \qquad (3) $$

on a real bounded interval space domain (α, β) or in a two dimensional square with Neumann Boundary Conditions (NBC). The notation I(x) and c(x) are used on purpose to emphasize their dependence in the space variable. The notations Δu stands for $\sum_{i=1}^{n} \frac{\partial^2 u}{\partial x_i^2}$. Around the center of the interval (α, β), I(x) and c(x) are set to a value for which the diffusion-less ODE underlying system would generate relaxation oscillations. Out of this region, I and c are set to a value for which the diffusion-less ODE would be in a stationary stable but excitable state. We allow those functions to be varied as a function of a parameter; this leads to a bifurcation path from a stationary state to propagation of oscillations for Nh-FHN. Previously, various studies have been conducted by some of the authors of the present paper, see [11–13]. In Section 4, we will present numerical simulations of this PDE to be compared to the dynamics of Equation (2). We will consider the following set of parameters in Section 4.1

$$ \begin{aligned} n &= 2, d_u = d_v = 1, b = 0, f(u) = -u^3 + 3u \\ I(x) &= 0 \\ c(x) &= 0 \quad \text{in a small ball at the center} \\ c(x) &= c_0 = -1.21 \quad \text{otherwise} \end{aligned} $$
and also

$$ 
\begin{array}{l}
n = 1, \Omega = (\alpha, \beta), d_u = d = 1, d_v = 0, b = 0, f(u) = -u^3 + 3u \\
I(x) = 0 \\
c(x) = p(x/\beta)^4 - 2p(x/\beta)^2 \\
\beta > 0, \alpha = -\beta
\end{array}
$$

in Section 4.2. All the partial differential equations will be solved with a RK4 discretization in time and a finite difference scheme in space.

### 3. Analysis of the FHR System

#### 3.1. A Short Background on FHN

The mechanisms at play in the FHR and Nh-FHN models rely primarily on the dynamics of the FHN model. Consequently, it is worth to first briefly describe the dynamics of the following classical FHN system:

$$ 
\left\{ 
\begin{array}{l}
\epsilon \frac{du}{dt} = f(u) - v + I \\
\frac{dv}{dt} = u - bv,
\end{array}
\right.
\qquad (4) $$

where $f(u) = -\frac{u^3}{3} + u$, $\epsilon$ is a small parameter, $0 < b < 1$ and $I > 0$.

**Proposition 1.** Equation (4) admits a unique stationary solution $(u^*, v^*)$ given by

$$ v^* = u^* / b $$

and $u^*$ is the unique solution of

$$ f(u) - u/b + I = 0. $$

It follows that $u^*$ is an increasing function of $I$; and the map $I \to u^*$ is a bijection from $(0, +\infty)$ to $(0, +\infty)$. Furthermore, as $I$ increases from $0$ to $+\infty$, $(u^*, v^*)$ is successively a repulsive node, a repulsive focus, an attractive focus and an attractive node. At $f'(u^*) = b\epsilon$ (i.e., $u^* = \sqrt{1 - b\epsilon}$) a Hopf bifurcation occurs.

Furthermore, when $u^*$ is close to zero, the system is known to exhibit relaxation oscillations. With this result in mind, one principle at play to obtain MMOs becomes clear: if we add a third variable that moves slowly, and in such a way that the dynamics for the system of the two first equations follow the loop: attractive focus, repulsive focus, relaxation oscillation and then return mechanism, MMOs should appear. The next subsection illustrates in more detail this idea to obtain MMOs as well as other phenomena related to the so-called canard solutions.

#### 3.1.1. A System with MMOs

We consider here the Equation (2) introduced above

$$ 
\left\{ 
\begin{array}{l}
\epsilon \frac{du}{dt} = f(u) - v + w + I, \\
\frac{dv}{dt} = u - bv \\
\frac{dw}{dt} = \epsilon(-u - w)
\end{array}
\right.
$$

with

$$ f(u) = -(1/3)u^3 + u, \epsilon = 0.1, I \in (1.3, 1.5), b = 0.8. $$

We focus on this system because it exhibits MMOs. We describe here two distinct mechanisms leading to alternance of small and large oscillations. The first mechanism corresponds to small oscillations related to the focus nature of the fixed point in the two-dimensional FHN Equation (4) as described in Section 1. The second mechanism is different; in this case, the trajectories clearly exhibits canard-type solutions, and the trajectories jump from the repulsive manifold toward alternatively one side or the other side of the stable manifold. Of note, this distinctive evolution toward the left or right side of the stable manifold oc-
curs within a tiny region of the repulsive manifold. We end the section with a remark on
Shilnikov Chaos.

3.1.2. MMOs and Focus

The first mechanism to obtain MMOs is as follows. Considering *w* as a parameter, the two first equations represent a classical FHN system, which goes through a Hopf bifurcation as *w* is varied. If *w* varies very slowly, in an interval corresponding to a focus for the FHN system, then we will observe focus-like dynamics in a neighborhood of the fixed point of the FHN subsystem. The focus is first attractive, then repulsive, until the trajectory falls into a relaxation oscillation type. This corresponds to small oscillations followed by a large oscillation. The dynamics are such that during the relaxation oscillation *w* returns close to its initial value. This is a return mechanism, from which a new cycle follows. This behavior is illustrated in Figure 2.

<figure><img src="image_7.png" /><figcaption>**Figure 2.** Mixed Mode Oscillations. This figure illustrates the asymptotic behavior of solutions of Equation (2) for $I = 1.45$. (a) In the left panel, time series of $u, v$ and $w$ are presented. (b) In the center, we illustrate the projection of the same trajectory on the plane $(u,v)$. (c) In the right, we illustrate the trajectory in the three dimensional phase space along with the critical manifold $v = f(u) + w + I$. This figure illustrates classical MMOs related to the focus nature of the fixed point of the first two equations (considering $w$ as a fixed parameter). In this case, the observed dynamics suggest the following successive states: attractive focus, repulsive focus, relaxation oscillations and return mechanism.</figcaption></figure>

**Remark 1.** *Another approach to describe the above dynamics relies on canards. The small oscillations are also to be seen as trajectories successively following the attractive and repulsive parts of the manifold till they exit the vicinity of the fold through a relaxation oscillation. We wanted here to emphasize an approach by a dynamical Hopf bifurcation, which appears to be relevant.*

3.1.3. MMOs and Canards

The second mechanism we want to discuss gives rise to a quite different type of
dynamics. In this case, there are no multiple small focus-like oscillations. Instead, after a
relaxation oscillation, the trajectory may enter a canard-type trajectory after crossing the
apex of the right part of the critical manifold; it follows the unstable manifold and leaves it
either to the left side—where it reaches the left part of the attractive critical manifold—or
to the right side—and reaches the right part of the attractive manifold. The latter case
corresponds to a small (or middle) oscillation while the first case gives a large relaxation
oscillation. This is illustrated in Figure 3. It is worth to describe in more detail the trajectory
represented in this figure. Two large oscillations are followed by a single small oscillation
where the trajectory follows the repulsive manifold before being attracted by the right
side of the attractive manifold. This is repeated four times, the fifth time, the second large
oscillation is replaced by a small oscillation (but larger than the other small oscillations).
This corresponds to a brutal change of direction in the trajectory occurring in a tiny zone of
the phase space. After that, the cycle is repeated. Introducing the letter M for the medium
oscillation, L for large and S for small, we could denote this occurrence: LLS-LLS-LLS-LLS-
LM. Note that during successive cycles LLS, the canards occurring during the second L, are significantly different in the phase space, and move forward in a specific direction. At the fourth time the canard exits the unstable manifold to the right. Of note, this discrimination between left and right exits has been used to generate rich dynamics in various contexts, see the recent papers [27,42].

<figure><img src="image_8.png" /><figcaption>**Figure 3.** (Complex) Mixed Mode Oscillations. In this case, there are no multiple small focus-like oscillations. Instead, after a relaxation oscillation, the trajectory may enter a canard-type trajectory after crossing the apex of the right part of the critical manifold; it follows the unstable manifold and leaves it either to the left side—where it reaches the left part of the attractive critical manifold—or to the right side—and reaches the right part of the attractive manifold. The latter case corresponds to a small (in fact a middle) oscillation, while the first case gives a large relaxation oscillation. It is worth describing in more detail the trajectory represented in this figure. Two large oscillations are followed by a single small oscillation where the trajectory follows the repulsive manifold before being attracted by the right side of the attractive manifold. This is repeated four times, and the fifth time, the second large oscillation is replaced by a middle-size oscillation. This corresponds to a change in the trajectory occurring in a tiny value range (see bottom left). After that, the cycle is repeated. Introducing the letter M for the medium oscillation, L for large and S for small, we could denote this occurrence: LLSLLSLLSLLSLM. (a) time series of *u*, *v* and *w*. (b) projection on the *u* − *v* plane. (c) projection on the *u* − *w* plane of four successive segments of the trajectory along the repulsive branch and the subsequent fast trajectory. The first three segments exit towards the left part of the attractive manifold. The fourth one exits toward the other side. (d) trajectory in the three dimensional phase space along with the critical manifold *v* = *f*(*u*) + *w* + *I*.</figcaption></figure>

**Remark 2.** In the next section, computations will show that the regimes of MMOs exhibited below correspond to a range of parameter value *I* for which the fixed point has two complex eigenvalues with a positive real part and one negative real value. This is a signature of Shilnikov Chaos, which relates to dynamics alternating phases on the attractive manifold and phases on the repulsive manifold supported by the complex eigenvalues [43,44] and references therein. Figure 4 gives an interesting glimpse of it. For the initial condition considered here, at the beginning, the solution follows the stable manifold corresponding to the negative eigenvalue. But afterwards, the trajectory exits the neighborhood of the fixed point with a focus-like dynamics; this corresponds to the two complex eigenvalues with positive real parts. After that, the trajectory never comes back in a neighborhood of the fixed point— the trajectory following the stable manifold corresponds to a
transient behavior—and as such the asymptotic observed dynamics do not relate to Shilnikov chaos despite the eigenvalues.

<figure><img src="image_9.png" /><figcaption>**Figure 4.** This figure illustrates a remarkable transient behavior. For the parameter values considered here, the stationary point has two complex eigenvalue with positive real parts and one negative real eigenvalue. Accordingly, for the initial conditions chosen here, the trajectory follows the stable manifold associated with negative eigenvalue. After that, the trajectories are repelled from the fixed point. The asymptotic behavior, however, does not relate so much with those eigenvalues. As such, the Shilnikov chaos is not relevant for these parameter values. (a) I = 1.45. (b) I = 1.3.</figcaption></figure>

## 3.2. Basic Stability Analysis

The following proposition results from simple computations.

**Proposition 2.** For any $b > 0$, Equation (2) admits a unique fixed point $(u^*, v^*, w^*)$ given by

$$v^* = u^*/b, w^* = -u^*$$

where $u^*$ is the unique solution of

$$f(u) - (1 + \frac{1}{b})u + I = 0.$$

A local stability analysis provides the following proposition.

**Proposition 3.** There exists $I^* \in (1, 2)$ such that at $I = I^*$, an a Hopf bifurcation occurs.

**Proof.** The Jacobian matrix $J(u)$ writes

$$\begin{pmatrix} \frac{f'(u)}{\epsilon} & -\frac{1}{\epsilon} & \frac{1}{\epsilon} \\ 1 & -b & 0 \\ -\epsilon & 0 & -\epsilon \end{pmatrix} \qquad (5)$$

which gives

$$\begin{aligned} \mathrm{Det}(J(u) - \lambda I) &= -\lambda^3 + \lambda^2(-\epsilon - b + \frac{f'(u)}{\epsilon}) + \lambda(f'(u) + b\frac{f'(u)}{\epsilon} - b\epsilon - \frac{1}{\epsilon} - 1) + b f'(u) - 1 - b \\ &= -\lambda^3 + \lambda^2(-\epsilon - b + \frac{f'(u)}{\epsilon}) + \lambda(f'(u) + b\frac{f'(u)}{\epsilon} - b\epsilon - \frac{1}{\epsilon} - 1) - bu^2 - 1 \end{aligned}$$

Thanks to simple algebraic computations (Routh-Hurwitz criterion and Cardan formula), one can prove that in the interval $(1, 2)$ an eigenvalue is real negative and the two other are complex conjugate; moreover, for $I < I^*$, the stationary point is unstable, and for $I > I^*$, it is stable. □

See Figure 5 for numerical illustration.
<figure><img src="image_10.png" /><figcaption>**Figure 5.** This figure gives an illustration of Proposition 3. In both panels the x-axis represents $f'(u^*)$ for $I \in [1, 2]$. (a) illustration of $-\Delta$ as defined in the Cardan formula. Since $\Delta$ is negative, the determinant admits two complex eigenvalues. (b) illustration of $a_2a_1 - a_0$ as defined in the Routh-Hurwitz criterion. It is positive for $I > I^*$. We also have $a_0 > 0$ and $a_2 = \epsilon + b - \frac{f'(u)}{\epsilon}$, which is positive for $I > I^*$, ($I^*$ is the first value for which $a_2a_1 - a_0$ or $a_2$ equals zero).</figcaption></figure>

### 3.3. Absorbing Set, Existence of Periodic Solution, Numerical Approximation

We assume $\epsilon > 0$ is small enough.

**Proposition 4.** System (2) admits an absorbing bounded set.

**Proof.** We define

$$ \psi(t) = \epsilon u^2 + v^2 + w^2 $$

We have

$$ \psi' = -\frac{u^4}{3} + u^2 + (1 - \epsilon)uw + Iu - bv^2 - \epsilon w^2 $$

By using young inequalities, we can prove that there exist two positive constants $K_1$ and $K_2$ such that

$$ \psi' \leq -K_1g + K_2. $$

Multiplying both sides by $e^{K_1t}$ and integrating leads to

$$ \psi(t) \leq e^{-K_1t}g_0 + \frac{K_2}{K_1}(1 - e^{K_1t}). $$

This completes the proof. □

**Theorem 1.** There exists $I_0$, such that for $I \in [0, I_0)$, the system admits a non-constant periodic solution.

**Proof.** We assume $I = 0$. The result extends to $[0, I_0)$ by continuity. We consider the Poincaré map $F = (F_1, F_2, F_3)$ from the manifold $M = \{(u, v, w); u = 0, v < 0, w \in \mathbb{R}\}$ to itself defined thanks to the flow of the ODE. For any $w$ such that $|w|$ is not too large, and for $\epsilon$ is small enough, the map is well-defined thanks to the slow-fast theory since the trajectories are close to relaxation oscillations and $w'$ is of $O(\epsilon)$. We know also that the trajectories will return at $M$ with a $v$-coordinate at $f(-1) + O(\epsilon)$ (with $f$ defined in Equation (2), and since $f'(-1) = 0$, see classical works on slow-fast systems such as for example [45]). From the Brouwer theorem, we deduce that for each fixed $w$, there exists $v$ such that

$$ F_2(0, v, w) = v $$

i.e., the $v$-coordinate is a fixed point. By this way, we define a continuous function $w \to \varphi(w) = v$ such that

$$ F_2(0, \varphi(w), w) = \varphi(w). $$
Next, we look for $w^*$ such that,

$$F_3(0, \varphi(w^*), w^*) = w^*.$$

Let us start with an initial condition $w_0 < 0$ close to 0. We have

$$F_3(0, \varphi(w_0), w_0) = w_0 e^{-\epsilon T} - \epsilon \int_0^T e^{-\epsilon(T-s)} u(s) ds$$

where $T$ is the time at which the flow ensued from $(0, \varphi(w_0), w_0)$ returns to $\mathcal{M}$. Replacing the exponential by its order 1 Taylor expansion, we found that adopting the notation $F_3(w_0)$ for the sake of simplicity

$$F_3(w_0) = w_0 - \epsilon w_0 T - \epsilon \int_0^T u(s) ds + o(\epsilon)$$

which gives $F_3(w_0) > w_0$, for $w_0 < 0$ in a small neighborhood of 0 and $\epsilon$ small enough. Analog arguments allow us to prove that $F_3(w_0) < w_0$ for $w_0 > 0$ and $\epsilon$ are small enough. We omit here the details of the computations that can be made explicit by changing the variable of integration from $s$ to $u$ and rely on the fact that the time spent on the right or left part of the stable manifold depends on the sign of $w$; this corresponds geometrically to the relative position of the nullclines. It follows that there exists $w^*$ such that $F_3(w^*) = w^*$ and therefore $(0, \varphi(w^*), w^*)$ is a fixed point of the Poincaré map. □

### 3.4. A Numerical Approximation for Small Oscillations

The aim of this section is to illustrate how the small oscillations observed during typical MMOs (as in Figure 2) can be locally captured by the dynamics of a "moving" focus. Observation of the small oscillations of Figure 2 show that oscillations in the $(u, v)$ plane are first decreasing in amplitude, then increasing till the trajectory leaves the vicinity of the fold line. Since $w$ moves very slowly, this corresponds to the fact that for fixed $w$ in the corresponding range, the stationary point of the two first equations is a focus, first stable and then unstable, see Proposition 1. In this section, we will approximate the dynamics of the full system by the dynamics of a simpler system for which computations can be made quite explicitly. We will first operate a change variables to translate our attention on the dynamics around the focus. Then, we will look at the linearization and compare the dynamics of the simple system with the original one. The next proposition gives the dynamics in a system of coordinates around $(u^*(w), v^*(w))$ where $(u^*(w), v^*(w))$ is the "stationary" solution of the subsystem in $(u, v)$ considering $w$ constant. Since $w$ is super slow, this system of coordinates is relevant.

**Proposition 5.** Let $(u^*(w), v^*(w))$ be the unique solution of

$$\begin{cases} 0 = f(u) - v + w + I, \\ 0 = u - bv \end{cases} \qquad (6)$$

Then, after a change of variables around $(u^*(w), v^*(w))$, and with appropriate notations, the trajectory of the solutions of (2) is given by:

$$\left\{ \begin{aligned} \frac{du}{dt} &= g(u^*, u) - v + \epsilon^2 \frac{1}{f'(u^*) - \frac{1}{b}} (u^* + u + w), \\ \frac{dv}{dt} &= u - bv + \epsilon \frac{1}{f'(u^*) - \frac{1}{b}} (u^* + u + w) \\ \frac{dw}{dt} &= -\epsilon(w + u^* + u) \end{aligned} \right. \qquad (7)$$

with

$$g(u^*, u) = (1 - u^{*2})u - u^*u^2 - \frac{u^3}{3} \qquad (8)$$
Next, our goal is to replace Equation (7) with a simpler system that mimics the small oscillations. Now, since $w$ is very slow, when $(u, v)$ is small in Equation (7), it is geometrically relevant to approximate solutions of Equation (7), during the specific period of small oscillations, by the system:

$$ \left\{ \begin{aligned} \frac{d\epsilon}{dt} &= (1 - u^{*2})u - v, \\ \frac{dv}{dt} &= u - bv, \\ \frac{dw}{dt} &= -\epsilon(w + u^{*} + u) \end{aligned} \right. \qquad (9) $$

What is interesting with system (9) is that it gives you a practical way to compute the number of oscillations involved. Let

$$ A = \begin{pmatrix} \frac{1}{\epsilon}(1 - u^{*2}) & -\frac{1}{\epsilon} \\ 1 & -b \end{pmatrix} $$

And let

$$ \beta = 0.5 \sqrt{4 \det A - (TrA)^2} $$

The following proposition holds.

**Proposition 6.** The number of small oscillations $n_o$ occurring during a time interval $(0, T)$ in Equation (9) is given by

$$ n_o = \frac{1}{2\pi} \int_0^T \beta(t) dt $$

**Proof.** The idea is to work with an approximate solution of Equation (9). We divide the time interval $(0, T)$ into the subdivision

$$ (kdt)_{k \in \{0,1,2,...,N\}}, \quad \text{with} \quad N = \frac{T}{dt} $$

Then, at each time step, for fixed $w$, it is possible to compute the solution of the linear subsystem made of the first two equations of (9) on the interval $(kdt, (k+1)dt)$. Eigenvalues are given by

$$ 0.5(TrA^{+} - \sqrt{(TrA)^2 - 4 \det A}) $$

and are complex conjugate when

$$ 1 - u^{*2} \approx b\epsilon $$

Classically, denoting $\alpha + i\beta$ one of this eigenvalues, we have after a new change of variables,

$$ \begin{pmatrix} \tilde{u} \\ \tilde{v} \end{pmatrix} = \begin{pmatrix} \alpha + b & \beta \\ 1 & 0 \end{pmatrix} $$

The solution of this last system writes after dropping the tilde

$$ u(t) + iv(t) = \exp(\alpha t) \exp(i\beta t) $$

Now, the number of small oscillations is given by

$$ \sum_k \frac{dt}{2\pi} \beta_k $$

Now, since the approximation converges towards the solution, the results follow from

$$ n_0 = \lim_{N \to +\infty} \sum_{k=0}^{N} \frac{dt}{2\pi} \beta_k = \frac{1}{2\pi} \int_{0}^{T} \beta(t) dt $$

□
**Remark 3.** The interest of this approximation is that it provides a geometrical approach to interpret the dynamics of small oscillations occurring in MMOs: the amplitude of the solutions first decrease and then increase. For system (7), this corresponds to $\alpha(u^*(t)) < 0$, in which case the amplitude decreases and then to $\alpha(u^*(t)) > 0$ in which case the oscillations increase.

**Remark 4.** Additionally, note that system (9) provides a simple way to generate small oscillations and control the number $n_0$; this approach captures the essential phenomenon at play in the generation of small oscillations in system (2). Figure 6 illustrates both oscillations for system (9) and (2). System (9) could be extended to generate this small oscillation recurrently with a reset as it is carried out in the classical Leaky Integrate and Fire models often used in applications [5,6,46]. The difference being that in LIF models, the equation is linear and one-dimensional with respect to a variable V (the voltage) with a reset that occurs when a V reaches a threshold value—while other inputs are typically included to drive this variable V. Here, the idea would be to reset the value w when it reaches the desired threshold value. In between the resets, the dynamics are oscillatory in the variables $(u, v)$.

<figure><img src="image_11.png" /><figcaption>**Figure 6.** This figure illustrates small oscillations for systems (9), in green, and (2) in purple. At the beginning the solutions are very close, but since we dropped the nonlinear terms the solutions are distinct after some time. Note, however, that the solutions of (9) capture two aspects of the small oscillations occurring for (2). First, it allows to generate and control in a simple way small oscillations by mimicking the behavior near a focus. Next, it generates oscillations that are first decreasing and then increasing, which corresponds to a dynamical Hopf-bifurcation behavior (TrA changes its sign and $(TrA)^2 - 4 \det A < 0$).</figcaption></figure>

### 3.5. Slow-Fast Analysis

In this section, we shall give some insights about the dynamics of Equation (2) thanks to a slow-fast approach. Setting $\epsilon = 0$ in Equation (2) after different time scalings provides the main dynamical picture. First, setting $\epsilon = 0$ gives

$$ \begin{cases} 0 = f(u) - v + w + I, \\ v_t = u - bv \\ w_t = 0 \end{cases} \qquad (10) $$

Then, after the time rescaling $t = \epsilon t'$ in Equation (2), dropping the ' and setting $\epsilon = 0$ gives

$$ \begin{cases} u_t = f(u) - v + w + I, \\ v_t = 0 \\ w_t = 0 \end{cases} \qquad (11) $$

Finally, rescaling with $t = \frac{1}{\epsilon}t'$, dropping the ' and setting $\epsilon = 0$ gives

$$ \begin{cases} 0 = f(u) - v + w + I, \\ 0 = u - bv \\ w_t = -(w + u) \end{cases} \qquad (12) $$

Equations (10) and (11) can be seen, for a fixed constant $w$, respectively, as the reduced and the layer systems of the classical 2D FHN system. We first consider the fast dynamics.
Outside of the critical manifold, the trajectories are given by the layer system (11), which is a one dimensional ODE in $u$. For any initial condition (unless we start at the repulsive point or in the stable manifold of a saddle), the trajectory will reach one of the attractive parts of the critical manifold $v = f(u) + w + I$, where $f'(u) < 0$.

Then, we look at the slow dynamics. Equation (10) is an ODE on the critical manifold. If the stationary point of Equation (10) is on the attractive part of the critical manifold, the solutions are well defined, and they will evolve to it. To further capture the evolution of the system, we need to consider the very slow motion given by Equation (12). If this stationary point is on the repulsive part with $f'(u) > 0$, then when it reaches the fold line $f(u) = 0$, system (10) is not defined, and the derivative explodes in finite time; this corresponds to a jump from one side of the attractive part of the critical manifold to the other one. If the stationary point is on the fold line, a more complex behavior (MMOs) can be expected, see [47,48] and references therein. An interesting and less classical insight from slow-fast analysis relates to the emergence of the sequence LLSLLSLLSLLSLM described in Section 3.1.3. In this case, the transition between canards exiting to the left (large oscillation) and canards exiting to the right (middle oscillation) plays a crucial role. Following the ideas in [49,50], we highlight hereafter some relevant computations. Consider Equation (2)

$$ \begin{cases} \epsilon u_t = f(u) - v + w + I \\ v_t = u - bv \\ w_t = -\epsilon(u + w) \end{cases} $$

We use the change of variables:

$$ u = 1 + \bar{u}, v = \bar{v} + \bar{v}, w = \bar{w} + \bar{w} $$

with

$$ \begin{aligned} f(1) - \bar{v} + \bar{w} + I &= 0 \\ 1 - b\bar{v} &= 0 \end{aligned} $$

which leads to the equation (dropping the bars)

$$ \begin{cases} \epsilon u_t = -3u^2 - u^3 - v + w, \\ v_t = u - bv \\ w_t = -\epsilon(\mu + u + w) \end{cases} $$

with

$$ \mu = 1 + \bar{w} = -1 + \frac{1}{b} - I $$

Next, we proceed to the change of variables

$$ u = \sqrt{\epsilon \bar{u}}, v = \epsilon \bar{v}, w = \epsilon \bar{w} $$

Further applying the change of time

$$ t = \frac{t}{\sqrt{\epsilon}} $$

and dropping the bars, we obtain

$$ \begin{cases} u_t = -3u^2 - \sqrt{\epsilon u^3} - v + w, \\ v_t = u - b\sqrt{\epsilon v} \\ w_t = -\sqrt{\epsilon}(\mu + \sqrt{\epsilon u} + \epsilon w) \end{cases} $$

Setting $\epsilon = 0$ gives

$$ \begin{cases} u_t = -3u^2 - v + w, \\ \frac{dv}{dt} = u \\ \frac{dw}{dt} = 0 \end{cases} \qquad (13) $$
**Proposition 7.** The point $(u, v) = (0, w)$ is a stationary point of center type for the projection of Equation (13) into the $(u, v)$-plane. Every trajectory passing through $(0, c)$ with $w < c < w + \frac{1}{6}$ is a periodic solution. Every solution passing through $(0, c)$ with $c \ge w + \frac{1}{6}$ satisfies $u(t) \to -\infty$ and $v(t) \to -\infty$ as $t \to t^*$, where $t^*$ is the maximal time for which the solution is defined. Furthermore, as $c$ goes to $w + \frac{1}{6}$, the diameter of the periodic solution goes to $+\infty$.

**Proof.** Without loss of generality, we work with $w = 0$, since the general case follows from the change of variable

$$-\tilde{v} + w = -v$$

Next, we note that the function

$$G(u, v) = (-3u^2 - v + 1/6) \exp(6v).$$

is a first integral of Equation (13) (with $w = 0$). It follows that the solutions of (13) are contained in the level sets of the function $G$. The proposition then follows from quite long yet elementary computations to explicit the details of the level sets of $G$. We avoid the details here. An illustrative picture is provided in Figure 7.

<figure><img src="image_12.png" /><figcaption>**Figure 7.** Solutions of the integrable Equation (13). The fixed point lies in the middle of limit cycles (nonlinear center). Above some threshold, solutions are no longer periodic but rather escape the vicinity of the stationary point. This illustrates the behavior of Figure 3: some canards exit to the left, while others exit to the right.</figcaption></figure>

## 4. Dynamics in the Nh-FHN Model

In this section, we shall consider the following Nh-FHN system:

$$\left\{ \begin{array}{l} \epsilon u_t = f(u) - v + I(x) + d_u \Delta u, \\ v_t = u - bv - c(x) + d_v \Delta v \end{array} \right. \qquad (14)$$

on a regular bounded domain $\Omega \subset \mathbb{R}^n$, $n \le 3$ with Neumann Boundary Conditions (NBC).

### 4.1. Existence and Uniqueness of the Stationary Solution

As far as we know, the existence result for Equation (14) in the case $d_u > 0$, $d_v > 0$ has not been proved. We consider here the case $d_u = d > 0$, $d_v = 0$. In this case, the equation writes

$$\left\{ \begin{array}{l} 0 = f(u) - v + I(x) + d \Delta u, \\ 0 = u - bv - c(x) \end{array} \right. \qquad (15)$$

When $b = 0$ (and $d_v = 0$), the existence of a stationary solution is straightforward since Equation (15) becomes $u = c$, $v = f(c) + I + d \Delta u$. We consider hereafter the case $b > 0$.

**Theorem 2.** Equation (15) admits a unique solution.
**Proof.** Equation (15) is equivalent to

$$ 
\left\{ 
\begin{array}{l}
0 = f(u) - (u/b) + (1/b)c(x) + I(x) + d\Delta u, \\
v = (1/b)(u - c(x))
\end{array}
\right. 
\qquad (16) 
$$

Then, we apply the result in [51]. □

## 4.2. Numerical Experiments

The numerical simulations presented in this paragraph have been performed using our own C++ program.

### 4.2.1. Filtering of Frequencies and Local Mixed Mode Oscillations (MMOs)

In this paragraph, we shall discuss some qualitative dynamics arising near a Hopf-bifurcation. We illustrate the numerical simulation of system (14) with the following set of parameters:

$$ 
\begin{array}{l}
n = 2, d_u = d_v = 1, b = 0, f(u) = -u^3 + 3u \\
I(x) = 0 \\
c(x) = 0 \quad \text{in a small ball at the center} \\
c(x) = c_0 = -1.21 \quad \text{otherwise}
\end{array} 
\qquad (17) 
$$

It is known from numerical experiments (see [11,12]) that if $c_0$ is decreased, the stationary solution of (14) is stable, and if it is very close to $-1$, the propagation of relaxation oscillations occurs. The value of $c_0$ considered here is close to a bifurcation point, which leads to more complex phenomena. Specifically, in this case, we observe the so-called mixed-mode oscillations for a center cell. Relaxation oscillations will propagate at a frequency smaller than the natural frequency of FHN in the oscillatory regime. This is illustrated in Figure 8. In this figure, panels (a), (b), (c) correspond to a fixed value of $x$ near the center. At this space location, $c(x) = 0$. Panels (d), (e), and (f) correspond to a value of $x$ near the right border. For this latter value of $x$, we have $c(x) = c_0 = -1.21$. We observe propagation of oscillations from the center toward the border. Note, however, that only large oscillations propagate. Small oscillations occurring in the center cells are filtered. Panels (d), (e), and (f) suggest the following interpretation of the local dynamics related to wave propagation. In panel (d), we represent $u$ for fixed $x$ as a function of time. We observe at this space location relaxation oscillations at a lower frequency than the oscillatory FHN (i.e., the diffusion-less system for $c = 0$). In panel (e), we represent $u, v$ and $\Delta u$ as functions of time, respectively, with red, dashed purple and blue colors. Through this panel, one can see how the wave propagation is seen locally. It corresponds to a wave of depolarization. Note that the diffusion remains close to zero for almost all times, which corresponds to a solution locally constant in space. When the neighbor cell jumps, the term $\Delta u$ becomes positive, because at this time we are at a minimum in space. This induces a jump for the considered cell toward the right part of the stable manifold (see panel (e)). At this point, the diffusion becomes negative because we are at a local maximum. Then, the diffusion comes back to almost zero again as the considered cell and its neighbors have approximately the same value. The same dynamic occurs symmetrically (the diffusion in this case is first negative and then positive), when the cell goes from the right part of the stable manifold toward the left part of the stable manifold. The dynamics of $v$ follows from the fact that $v_t = u - c + d_v \Delta v$. Panel (f) illustrates the dynamics of $u, v$ and $\Delta u$ in the three-dimensional phase space. The “critical” manifold $v = f(u) + d_u \Delta u$ is also represented where $\Delta u$ is considered as a variable. It illustrates the oscillations of relaxation for the cells near the border. Illustrations in panels (a), (b), (c) are analog. For this cell the dynamics resemble MMOs often pictured in three dimensional systems.
<figure><img src="image_13.png" /><figcaption>**Figure 8.** Simulation of Equation (14) with the set of parameters given in (17). The first row corresponds to a fixed value of *x* near the center for which *c*(*x*) = 0. The second row corresponds to a value of *x* near the left border for which *c*(*x*) = *c*₀ = −1.21. In the first column, we represent *u* as a function of time. In the **second column**, we represent *u*, *v* and *Δu* as functions of time, respectively, in green, blue and violet. In the **third column**, we represent the trajectory *u*, *v*, *Δu* along with the manifold *v* = *f*(*u*) + *d*<sub>*u*</sub>*Δu* (where *Δu* is seen as the third variable). These figures illustrate the apparition of MMOs for center cells and the propagation of relaxation oscillations from the center cells toward the border with a filtration of small oscillations.</figcaption></figure>

4.2.2. Fade of Wave Propagation (Death Spot)

In this paragraph, we discuss another phenomenon arising as we vary a parameter (denoted by *p* in this paragraph). We present here the numerical simulation of system (14) with the following set of parameters

$$
\begin{align}
n &= 1, \Omega = (\alpha, \beta), d_u = d = 1, d_v = 0, b = 0, f(u) = -u^3 + 3u \nonumber \\
I(x) &= 0 \nonumber \\
c(x) &= p(x/\beta)^4 - 2p(x/\beta)^2 \tag{18} \\
\beta &> 0, \alpha = -\beta \nonumber
\end{align}
$$

In this case, again, as in the previous paragraph, if *p* is close enough to −1, waves propagate from the center at the frequency of the diffusion-less system (with *c* = 0). On the other hand, if *p* is decreased enough, the system generally converges to a stationary solution. In between, for some range of *p*, regular waves propagate from center but fail to propagate at some point in space. The existence of Hopf Bifurcation has been proved for this specific case, see [12]. This type of failure of propagation, also referred as death-spot has been described previously, see for example [52]. For relevance in biological context, we refer to [10]. Related phenomena have also recently been considered in chains of kicked FHN
neurons, see [42]. The dynamics under consideration are illustrated here in Figure 9. At the center, oscillations are as in the ODE diffusion-less system. But at some space location, the oscillation fails to propagate. For some intermediate cells, we observe alternations between medium and larger oscillations. In this case, every other oscillation, the amplitude will be shortened— note the difference in the *u*<sub>xx</sub> time series. In this figure, panels (a), (b), and (c) correspond to a fixed value of *x* near the center. At this space location *c(x)* = 0. Panels (d), (e), and (f) correspond to an intermediate value of *x* between the center and the left border. In panels (a) and (d), we represent *u* for fixed *x* as a function of time. In panels (b) and (e), we represent *u*, *v* and *u*<sub>xx</sub> as functions of time, respectively, with the colors red, dashed purple and blue. In panels (c) and (f), we illustrate the dynamics of *u*, *v* and *u*<sub>xx</sub> in the three-dimensional phase space. The “critical” manifold *v* = *f(u) + du*<sub>xx</sub> is also represented where *u*<sub>xx</sub> is considered as a variable. Figure 10 illustrates *u* as a function of time. From left to right, we represent four different space locations with the same range of amplitude: left panel corresponds to a center cell, the far right corresponds to a cell close to the left border, the two other panels correspond to intermediate cells. This illustrates the fade of the wave propagation. Note the difference of these complex oscillations as compared to the MMOs described in the previous paragraph.

<figure><img src="image_14.png" /><figcaption>**Figure 9.** This figure illustrates the simulation of Equation (14) with the set of parameters given in (18). Panels (a–c) correspond to a fixed value of *x* near the center for which *c(x)* = 0. Panels (d–f) correspond to a value of *x* at an intermediate location between the left border and the center. In (a) and (d), we represent *u* as a function of time. In (b) and (e), we represent *u*, *v* and *u*<sub>xx</sub> as functions of time, respectively, in green, blue and violet. In (c) and (f), we represent the trajectory *u*, *v*, *u*<sub>xx</sub> along with the manifold *v* = *f(u) + du*<sub>xx</sub> (where *u*<sub>xx</sub> is seen as the third variable).</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 10.** This figure illustrates the fade of the wave propagation. We represent here the value of *u* as a function of time. From left to right, we consider four different space locations with the same range of amplitude: left panel corresponds to a center cell, the far right corresponds to a cell close to the left border, the two other panels correspond to intermediate cells.</figcaption></figure>

**5. Discussion**

Dynamical systems have proven to be a crucial tool in the description of phenomena observed in neuroscience. This includes complex oscillations, rhythms, wave propagation and pattern formations. As discussed in the introduction, models such as the FHN equations have been used successfully to visualize the emergence of those phenomena and highlighted new ways of biological interpretation. For example, cardiac tachycardia has been associated with spiral waves, which are a classical phenomena studied in FHN systems. Yet, the mechanisms inducing the various dynamics are far from being well-understood. It is necessary to gain a better understanding of the mathematical principles that govern those dynamics to be able to reproduce the biological observed phenomena. In the present article, we focused on the analysis of the FHR and Nh-FHN models. We have provided original theoretical and numerical tools to understand the emergence of complex oscillations and their propagation across excitable media. This included an original theoretical and numerical description of small oscillations and canards leading to different kind of MMOs for the FHR system and an original description of local MMOS, filtering frequencies and fade of wave propagation for Nh-FHN. A better qualitative analysis of the various complex oscillations and their propagation arising in Nh-FHN is an entire topic of research. One possible fruitful way to investigate is to work on finite-dimensional subspaces in which solutions effectively lie. Such insights can be applied to better simulate wave propagation, patterns and better understand their failure and their complexity in the context of neuroscience and biology.

**6. Conclusions**

In this article we have considered the FHR and Nh-FHN models, discussed their relevance in neuroscience context and provided a theoretical and numerical analysis of complex oscillations and their propagation. Future work perspectives include a more detailed investigation of reduction of the Nh-FHN dynamics to ODE models and their various applications to neuroscience and biology. PDEs are infinite dimensional systems; however, asymptotically, the solutions may lie in an attractor, which can be of a much smaller size. Some preliminary unpublished results, see also [13], indicate that the complex oscillations observed in Nh-FHN might be obtained with finite-size ODE models. This is an interesting perspective that would allow a better understanding of the dynamics observed and their biological interpretation.
**Author Contributions:** Conceptualization, all authors.; methodology, B.A.; software, B.A.; validation, all authors; investigation, all authors; writing—original draft preparation, all authors.; writing—B.A. and M.A.A.-A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the LMAH, Région Normandie, CNRS project IEA00134, and the HSM.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

## References

1. FitzHugh, R. Impulses and Physiological States in Theoretical Models of Nerve Membrane. *Biophys. J.* **1961**, *1*, 445–466. [CrossRef]

2. Nagumo, J.; Arimoto, S.; Yoshizawa, S. An Active Pulse Transmission Line Simulating Nerve Axon. *Proc. IRE 1962*, *50*, 2061–2070. [CrossRef]

3. Hodgkin, A.L.; Huxley, A.F. A quantitative description of membrane current and its application to conduction and excitation in nerve. *J. Physiol.* **1952**, *117*, 500–544. [CrossRef] [PubMed]

4. Lyu, C.; Abbott, L.F.; Maimon, G. Building an allocentric travelling direction signal via vector computation. *Nature* **2021**, *601*, 92–97. [CrossRef]

5. Chariker, L.; Shapley, R.; Young, L.S. Orientation Selectivity from Very Sparse LGN Inputs in a Comprehensive Model of Macaque V1 Cortex. *J. Neurosci.* **2016**, *36*, 12368–12384. [CrossRef] [PubMed]

6. Chariker, L.; Shapley, R.; Young, L.S. Rhythm and Synchrony in a Cortical Network Model. *J. Neurosci.* **2018**, *38*, 8621–8634. [CrossRef] [PubMed]

7. Volpert, V.; Xu, B.; Tchechmedjiev, A.; Harispe, S.; Aksenov, A.; Mesnildrey, Q.; Beuter, A. Characterization of spatiotemporal dynamics in EEG data during picture naming with optical flow patterns. *Math. Biosci.* **Eng. 2023**, *20*, 11429–11463. [CrossRef]

8. Panfilov, L.L.A.; Dierckx, H. Impact of Electrode Orientation, Myocardial Wall Thickness and Myofiber Direction on Intracardiac Electrograms: Numerical Modeling and Analytical Solutions. *Front. Physiol. Sec. Card. Electrophysiol.* **2023**, *14*. [CrossRef]

9. Harlaar, N.; Dekker, S.O.; Zhang, J.; Snabel, R.R.; Veldkamp, M.W.; Verkerk, A.O.; Fabres, C.C.; Schwach, V.; Lerink, L.J.S.; Rivaud, M.R.; et al. Conditional immortalization of human atrial myocytes for the generation of in vitro models of atrial fibrillation. *Nat. Biomed. Eng.* **2022**, *6*, 389–402. [CrossRef]

10. Maia, P.D.; Kutz, J.N. Identifying critical regions for spike propagation in axon segments. *J. Comput. Neurosci.* **2013**, *36*, 141–155. [CrossRef]

11. Ambrosio, B.; Francoise, J.P. Propagation of bursting oscillations. *Philos. Trans. R. Soc. A Math. Phys. Eng. Sci.* **2009**, *367*, 4863–4875. [CrossRef]

12. Ambrosio, B. Hopf Bifurcation in an Oscillatory-Excitable Reaction–Diffusion Model with Spatial Heterogeneity. *Int. J. Bifurc. Chaos* **2017**, *27*, 1750065. [CrossRef]

13. Ambrosio, B. Qualitative analysis of certain reaction-diffusion systems of the FitzHugh-Nagumo type. *Evol. Equ. Control Theory* **2023**. [CrossRef]

14. Kopell, N.; Howard, L.N. Plane Wave Solutions to Reaction-Diffusion Equations. *Stud. Appl. Math.* **1973**, *52*, 291–328. [CrossRef]

15. Rinzel, J.; Keller, J.B. Traveling Wave Solutions of a Nerve Conduction Equation. *Biophys. J.* **1973**, *13*, 1313–1337. [CrossRef]

16. Rinzel, J. Spatial stability of traveling wave solutions of a nerve conduction equation. *Biophys. J.* **1975**, *15*, 975–988. [CrossRef] [PubMed]

17. Rinzel, J.; Terman, D. Propagation Phenomena in a Bistable Reaction-Diffusion System. *SIAM J. Appl. Math.* **1982**, *42*, 1111–1137. [CrossRef]

18. Jones, C.K.R.T. Stability of the travelling wave solution of the FitzHugh-Nagumo system. *Trans. Amer. Math. Soc.* **1984**, *286*, 431–469. [CrossRef]

19. Ermentrout, G.B.; Rinzel, J. Reflected Waves in an Inhomogeneous Excitable Medium. *SIAM J. Appl. Math.* **1996**, *56*, 1107–1128. [CrossRef]

20. Carter, P.; Sandstede, B. Fast Pulses with Oscillatory Tails in the FitzHugh-Nagumo System. *SIAM J. Math. Anal.* **2015**, *47*, 3393–3441. [CrossRef]

21. Cornwell, P.; Jones, C.K.R.T. On the Existence and Stability of Fast Traveling Waves in a Doubly Diffusive FitzHugh-Nagumo System. *SIAM J. Appl. Dyn. Syst.* **2018**, *17*, 754–787. [CrossRef]

22. Volpert, A.; Volpert, V.; Volpert, V. *Traveling Wave Solutions of Parabolic Systems*; American Mathematical Society: Providence, RI, USA, 1994. [CrossRef]
23. Sandstede, B. Chapter 18—Stability of Travelling Waves. In *Handbook of Dynamical Systems*; Fiedler, B., Ed.; Elsevier Science: Amsterdam, The Netherlands, 2002; Volume 2, pp. 983–1055. [CrossRef]

24. Volpert, V.; Petrovskii, S. Reaction-diffusion waves in biology. *Phys. Life Rev.* **2009**, *6*, 267-310. [CrossRef]

25. Rinzel, J. A Formal Classification of Bursting Mechanisms in Excitable Systems. In Proceedings of the International Congress of Mathematicians, Berkeley, CA, USA, 3–11 August 1986.

26. Rinzel, J. A Formal Classification of Bursting Mechanisms in Excitable Systems. In *Mathematical Topics in Population Biology, Morphogenesis and Neurosciences: Proceedings of an International Symposium Held in Kyoto, Japan, 10–15 November 1985*; Teramoto, E., Yumaguti, M., Eds.; Springer: Berlin/Heidelberg, Germany, 1987; pp. 267–281.

27. Desroches, M.; Rinzel, J.; Rodrigues, S. Classification of bursting patterns: A tale of two ducks. *PLoS Comput. Biol.* **2022**, *18*, e1009752. [CrossRef] [PubMed]

28. Izhikevich, E.M. Synchronization of elliptic bursters. *SIAM Rev.* **2001**, *43*, 315–344. [CrossRef]

29. Mondal, A.; Sharma, S.K.; Upadhyay, R.K.; Mondal, A. Firing activities of a fractional-order FitzHugh-Rinzel bursting neuron model and its coupled dynamics. *Sci. Rep.* **2019**, *9*, 15721. [CrossRef]

30. Wojcik, J.; Shilnikov, A. Voltage Interval Mappings for an Elliptic Bursting Model. In *Nonlinear Systems and Complexity*; Springer International Publishing: New York, NY, USA, 2015; pp. 195–213.

31. Clément, F.; Françoise, J.P. Mathematical Modeling of the GnRH Pulse and Surge Generator. *SIAM J. Appl. Dyn. Syst.* **2007**, *6*, 441–456. [CrossRef]

32. Vidal, A.; Clément, F. A Dynamical Model for the Control of the Gonadotrophin-Releasing Hormone Neurosecretory System. *J. Neuroendocrinol.* **2010**, *22*, 1251–1266. [CrossRef]

33. Albizuri, J.U.; Desroches, M.; Krupa, M.; Rodrigues, S. Inflection, Canards and Folded Singularities in Excitable Systems: Application to a 3D FitzHugh-Nagumo Model. *J. Nonlinear Sci.* **2020**, *30*, 3265–3291. [CrossRef]

34. Desroches, M.; Krupa, M.; Rodrigues, S. Inflection, canards and excitability threshold in neuronal models. *J. Math. Biol.* **2012**, *67*, 989–1017. [CrossRef]

35. Amir, R.; Michaelis, M.; Devor, M. Burst Discharge in Primary Sensory Neurons: Triggered by Subthreshold Oscillations, Maintained by Depolarizing Afterpotentials. *J. Neurosci.* **2002**, *22*, 1187–1198. [CrossRef]

36. Bertram, R. Mathematical Modeling in Neuroendocrinology. *Compr. Physiol.* **2015**. [CrossRef]

37. Fletcher, P.A.; Thompson, B.; Liu, C.; Bertram, R.; Satin, L.S.; Sherman, A.S. $Ca^{2+}$ release or $Ca^{2+}$ entry, that is the question: What governs $Ca^{2+}$ oscillations in pancreatic $\beta$ cells? *Am. J. Physiol.-Endocrinol. Metab.* **2023**, *324*, E477–E487. [CrossRef]

38. Desroches, M.; Guckenheimer, J.; Krauskopf, B.; Kuehn, C.; Osinga, H.M.; Wechselberger, M. Mixed-Mode Oscillations with Multiple Time Scales. *SIAM Rev.* **2012**, *54*, 211–288. [CrossRef]

39. Letson, B.; Rubin, J.E.; Vo, T. Analysis of Interacting Local Oscillation Mechanisms in Three-Timescale Systems. *SIAM J. Appl. Math.* **2017**, *77*, 1020–1046. [CrossRef]

40. Kaklamanos, P.; Popović, N.; Kristiansen, K.U. Bifurcations of mixed-mode oscillations in three-timescale systems: An extended prototypical example. *Chaos Interdiscip. J. Nonlinear Sci.* **2022**, *32*, 013108. [CrossRef] [PubMed]

41. Kuehn, C. *Multiple Time Scale Dynamics*; Springer International Publishing: Berlin/Heidelberg, Germany, 2015. [CrossRef]

42. Ambrosio, B.; Mintchev, S.M. Periodically kicked feedforward chains of simple excitable FitzHugh-Nagumo neurons. *Nonlinear Dyn.* **2022**, *110*, 2805–2829. [CrossRef]

43. Glendinning, P.; Sparrow, C. Local and global behavior near homoclinic orbits. *J. Stat. Phys.* **1984**, *35*, 645–696. [CrossRef]

44. Shilnikov, L.; Shilnikov, A. Shilnikov bifurcation. *Scholarpedia* **2007**, *2*, 1891. [CrossRef]

45. Szmolyan, P.; Wechselberger, M. Canards in R3. *J. Differ. Equ.* **2001**, *177*, 419–453. [CrossRef]

46. Abbott, L. Lapicque’s introduction of the integrate-and-fire model neuron (1907). *Brain Res. Bull.* **1999**, *50*, 303–304. [CrossRef]

47. Krupa, M.; Szmolyan, P. Extending Geometric Singular Perturbation Theory to Nonhyperbolic Points—Fold and Canard Points in Two Dimensions. *SIAM J. Math. Anal.* **2001**, *33*, 286–314. [CrossRef]

48. Krupa, M.; Ambrosio, B.; Aziz-Alaoui, M.A. Weakly coupled two-slow-two-fast systems, folded singularities and mixed mode oscillations. *Nonlinearity* **2014**, *27*, 1555–1574. [CrossRef]

49. Krupa, M.; Popović, N.; Kopell, N. Mixed-mode oscillations in three time-scale systems: A prototypical example. *SIAM J. Appl. Dyn. Syst.* **2008**, *7*, 361–420. [CrossRef]

50. De Maesschalck, P.; Kutafina, E.; Popović, N. Three time-scales in an extended Bonhoeffer-van der Pol oscillator. *J. Dyn. Differ. Equ.* **2014**, *26*, 955–987. [CrossRef]

51. Boccardo, L.; Croce, G. *Elliptic Partial Differential Equations*; De Gruyter: Gottingen, Germany, 2013.

52. Yanagita, T.; Nishiura, Y.; Kobayashi, R. Signal propagation and failure in one-dimensional FitzHugh-Nagumo equations with periodic stimuli. *Phys. Rev. E* **2005**, *71*, 036226. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.