Qualitative Analysis of certain Reaction-Diffusion Systems of the FitzHugh-Nagumo type

TL;DR

This paper conducts a qualitative analysis of nonlinear Reaction-Diffusion systems, specifically the FitzHugh-Nagumo type, exploring solution behaviors, bifurcations, and their implications in Neuroscience.

Contribution

It introduces a non-homogeneous FitzHugh-Nagumo model and analyzes a related toy model to understand solution convergence and bifurcation phenomena.

Findings

Solutions converge to fixed points or periodic orbits

Existence of a cascade of Hopf bifurcations

Insights into neural excitability and oscillations

Abstract

This article aims to provide insights into the qualitative analysis of some nonlinear Reaction-Diffusion (RD) systems arising in Neuroscience. We first introduce a non-homogeneous FitzHugh-Nagumo (nhFHN) featuring excitability and oscillatory properties. Then, we discuss the qualitative analysis of a toy model related to nhFHN. In particular, we focus on the convergence of solutions of the toy model toward different solutions (fixed point, periodic) and show the existence of a cascade of Hopf bifurcations. Finally, we connect this analysis to the nhFHN system.

Figures

Equations

$$\left\{\begin{array}[]{rcl}\epsilon u_{t}&=&f(u)-v+du_{xx},\,(x,t)\in(a,b)\times(0,+\infty)\\ v_{t}&=&u-c(x)\\ \end{array}\right.$$

$$\mathcal{H}=L^2(a,b)\times L^2(a,b)$$

$$\mathcal{V}=H^1(a,b)\times H^1(a,b)\text{mbox}where{H}^1(a,b)\text{mbox}denotestheusualSobolevspace.$$

$$||\cdot ||\text{mbox}denotestheusualnormon\mathcal{H}.$$

$$(u,v)\in C({\mathbb{R}}^{+},\mathcal{H})$$

$$\forall \mu >0,\forall T>0,\forall t\in (\mu ,t)u(t)\in L^{\infty }(a,b)$$

$$u\in L^2((0,T),H^1)\cap L^4((a,b)\times (0,T))$$

$$\left\{\begin{array}[]{rcl}u_{t}&=&\alpha u-u^{3}-v+u_{xx}\\ v_{t}&=&u\end{array}\right.$$

$$\left\{\begin{array}[]{rcl}u_{t}&=&\alpha u-v+u_{xx}\\ v_{t}&=&u\end{array}\right.$$

$${\phi }_0(x)=1,\text{mbox}and\forall k\in {\mathbb{N}}^{\ast }{\phi }_k(x)=\sqrt{2}\cos (k\pi x).$$

$$-({\phi }_k{)}_{xx}={\lambda }_k{\phi }_k$$

$$({\phi }_k{)}_x(0)=({\phi }_k{)}_x(1)=0,$$

$${\lambda }_k=k^2{\pi }^2.$$

$$u(t)=\sum _{k=0}^{\infty }{u}_k(t){\phi }_k,v(t)=\sum _{k=0}^{\infty }{v}_k(t){\phi }_k$$

$$(E_{k})\left\{\begin{array}[]{rcl}u_{kt}&=&(\alpha_{k}-\lambda_{k})u_{k}-v_{k}\\ v_{kt}&=&u_{k}\end{array}\right.$$

$$A_k=\left(\begin{array}{cc}\alpha -{\lambda }_k& -1\\ 1& 0\end{array}\right)$$

$$\sigma_{k}^{1}=\frac{1}{2}\bigg{(}\alpha-\lambda_{k}-\sqrt{(\alpha-\lambda_{k})^{2}-4}\bigg{)},\,\sigma_{k}^{2}=\frac{1}{2}\bigg{(}\alpha-\lambda_{k}+\sqrt{(\alpha-\lambda_{k})^{2}-4}\bigg{)}.$$

$$\underset{k\to +\infty }{\lim }{\sigma }_k^1=-\infty \text{mbox}and\underset{k\to +\infty }{\lim }{\sigma }_k^2=0^{-},$$

$${\sigma }_k^2=\frac{1}{\alpha -{\lambda }_k}+o(\frac{1}{\alpha -{\lambda }_k}).$$

$$\underset{t\to +\infty }{\lim }||(u,v)(t)||=0.$$

$$\sum _{k=0}^N(|u_k(t){|}^2+|v_k(t){|}^2)\le e^{-\delta t}\sum _{k=0}^N(|u_k(0){|}^2+|v_k(0){|}^2),$$

$$u_k(0)={\int }_0^{1}u(x,0){\phi }_k(x)dx,v_k(0)={\int }_0^{1}v(x,0){\phi }_k(x)dx.$$

$$\sum _{k=N+1}^{+\infty }(|u_k(0){|}^2+|v_k(0){|}^2)<\frac{ϵ}{2}.$$

$$\frac{d}{dt}(|u_k(t){|}^2+|v_k(t){|}^2)\le 2(\alpha -{\lambda }_k)|u_k(t){|}^2\le 0,$$

$$\forall t>0,\sum _{k=N+1}^{+\infty }(|u_k(t){|}^2+|v_k(t){|}^2)<\frac{ϵ}{2}.$$

$$||(u,v)(t)||<ϵ.$$

$$\underset{t\to +\infty }{\lim }||(u,v)(t)-{\phi }_k(u_k(t),v_k(t))||=0.$$

$$\underset{t\to +\infty }{\lim }||(u,v)(t)||=+\infty .$$

$$\underset{t\to +\infty }{\lim }||(u,v)(t)||=0.$$

$$\underset{t\to +\infty }{\lim }||(u,v)(t)||=+\infty .$$

Full text

1. Introduction

This article is motivated by the qualitative analysis of the following RD system

[TABLE]

with $f(u)=-u^{3}+3u$, $\epsilon$ small, and Neumann Boundary conditions (NBC). In the Neuroscience context, $u$ represents a potential and $v$ a recovery variable while $\Omega=(a,b)$ relates to the spatial domain, typically the axon of one single neuron. This system provides a simple model generating excitable and oscillatory behavior. Thanks to the non-homogeneity in $c(x)$, it can feature both properties in one single PDE, see [Amb-2009, Amb-2009-thesis, Amb-2017]. This is well suited for Neuroscience applications where it is relevant to characterize the response of an excitable tissue as the stimulus intensity or frequency varies. It also raises deep mathematical questions. Before delving into more details, let us recall some contextual facts which motivate the introduction of this non-homogeneous model. In 1952, see [HH], Hodgkin and Huxley (HH) provided their seminal model, describing the electrical pulse propagation along the giant squid neuron axon. They were awarded with the 1963 prize in Physiology and Medicine (together with J. Eccles) for their work which included both a modeling part and physiological experiments. Of note, their original article contained a four ODE model as well as a spatially extended RD system (HH model) based on the cable analogy. We refer, for example to [BookDay2001, BookErm2010, BookIzh2006] for further details on the model. Later, in 1961, see [Fit-1961], inspired by the HH model, R. FitzHugh introduced a two-dimensional ODE to provide a simpler model that can reproduce excitable and oscillatory behavior. He started from the Van der pol model (VDP), see [VDP26], known for its oscillatory behavior. On the other hand, in [Bon-1948], K.F. Bonhoeffer described the dynamics of a two-dimensional ODE system featuring oscillatory or excitable properties. Following the ideas of Bonhoeffer, Fitzhugh modified the VDP model to obtain a 2d equation able to exhibit the oscillatory and excitable dynamics observed in the HH model. He named it the BonHoeffer Van der Pol model. One year later, Nagumo et al. provided an analog electrical circuit in [Nagumo]. The model is now widely known as the FitzHugh-Nagumo model. With these ideas in mind, a simple model which can provide excitable and oscillatory properties is the diffusion-less underlying ODE system of (1) with a constant $c$. Indeed, it is well known that if $|c|<1$, the unique stationary point is unstable and that there exists a unique limit cycle that attracts all the non-constant trajectories. At $|c|=1$, a supercritical Hopf bifurcation occurs, and for $|c|\geq 1$ the stationary point is globally asymptotically stable. Furthermore, in the latter case, the system is excitable because choosing an initial condition below a certain threshold will induce a large excursion in the phase plane, which is akin to a pulse for the first variable. Our motivation in studying the diffusive system (1) is to investigate the qualitative behavior when $c(x)$ is set within a range of values giving raise to excitable and oscillatory dynamics. For example, will the signal propagate if some cells are oscillatory at the center of the domain and anywhere else the cells are excitable? This question, which relates to an interplay between excitability and oscillatory properties, has been partially addressed in a series of papers, see [Amb-2016, Amb-2009, Amb-2014]. This typically leads to a bifurcation path from a stationary state to the propagation of oscillations when peripheral cell’s excitability increases. A complementary question is the characterization of the basin of attraction of patterns, see [Amb-2016]. Studying the propagation of oscillations initiated at the center throughout the domain relates to the theory of traveling waves. While our approach here is quite different, it is worth mentioning some contributions of this theory that applies to the FHN equations. Since the first studies on the Fisher-KPP equation [Fis1937, Kol1937, McK1975], the topic has indeed aroused considerable interest, see for example [BookVolpert1994]. The particular case of the wave propagation phenomenon in the diffusive FHN, with $\Omega={\mathbb{R}}$ has also been intensively studied for a few decades. One of the first contributions can be found in [Rin1973, Rin1975]. In these articles, J. Rinzel and coauthors, following the ideas in McKean [McK1970] studied a FHN RD system with piecewise linear nonlinearity (instead of the cubic nonlinearity). They provided explicit computations of periodic solutions, pulses, propagation speeds, and stability results. The case where the underlying ODE is bistable was also considered later in [Rin1982a]. Around the same period, in a series of four articles, J.W. Evans provided a more theoretical analysis of a general model of nerve conduction, see [Eva1971, Eva1972a, Eva1972b, Eva1975]. In [Eva1972a], he defines the so-called Evans function, which would become an essential tool for traveling wave stability analysis see [San2002, Bar2018]. In 1978, Rauch and Smoller, see [Rau1978], provided an analysis of FHN with $\Omega={\mathbb{R}}$ which includes the well-posedness of the problem in various Banach spaces, as well as the stability of the stationary solution $(0,0)$ with various techniques (linearization, contracting rectangles…). A few years later, see [Jon1984], C. Jones, relying on the papers of Evans, provided a detailed analysis focusing on the stability of the traveling waves of the FHN diffusive equation. Since then, many studies have been devoted to the characterization of different properties of traveling waves. This includes the proof of the existence of pulses of monotonic and periodic tales and the construction of solutions in a slow-fast context thanks to singular perturbation theory or via asymptotic expansions. Relevant articles to these matters are for example: [Car1977, Has1976, Bez2017, Car2015, Car2016, Car2018, Cor2018, Erm1996, Jon1991, Kal1993, Guc2009, Guc2010, Den1991, Cha2007, Kru1997]. Those ideas and techniques are relevant to the present work, but our approach here is somewhat different. First, we work on a bounded interval subjected to NBC and not in ${\mathbb{R}}$. Moreover our first motivation is to understand the dynamics that can emerge from the inhomogeneity in $c$. This relates by many aspects to the simpler case treated in [Amb-2014] where two coupled ODE FHN systems with different values of $c$ were studied. There is in fact several approaches to track the complex oscillations arising in nhFHN. The major contribution of the present article is to consider a toy model from which interesting oscillatory can be proved rigorously by projection in the eigenfunctions and to emphasize the analogy with nhFHN. Before to proceed further, it is also necessary to mention classical articles which consider RD systems related to (1) in a bounded domain and NBC. Those include for example [Marion, BookTemam, BookRobinson] and references therein cited. One must remark however, that in comparison to the systems considered there, system (1) has two characteristics that, combined, bring technical difficulties:

1. (1)

the system is partially dissipative in the sense that only the first equation diffuses 2. (2)

there is no dissipative term of the form $-\delta v$ in the second equation

As mentioned before, we are particularly interested in bifurcation phenomena from which periodic solutions arise as well as in the stability of these solutions. To this end, we will first introduce a toy model which allows proving interesting results with explicit computations. Then, we will emphasize how this relates to system (1). Technically, for (1), due to the nonhomogeneity arising with $c(x)$, the spectrum analysis, which relies on the eigenvalues of $u_{xx}$ for the toy model, is replaced by a Sturm-Liouville analysis. The remaining of the article is divided as follows: we analyze the toy model in section two; section three is devoted to the analysis of (1). The concluding remarks follow in section four.

**Notations and General Framework

**

A classical approach is to work in the $L^{2}$ setting, see [Amb-2009-thesis, BookRobinson, BookTemam, BookJLL, BookRothe]. Accordingly, we set the following notations:

[TABLE]

[TABLE]

[TABLE]

Subscripts $t$ and $x$ denote time and space derivatives, respectively. We will not write $(a,b)$ for simplicity if there is no ambiguity. For the reader’s convenience, we recall a result for the wellposedness of (1). We omit the proof which relies on the Galerkin method; see, for example, [Amb-2009-thesis, BookRobinson, BookTemam, BookJLL].

Theorem 1**.**

For initial conditions (IC) $(u,v)(x,0)$ given in ${\mathcal{H}}$, there exists a unique solution $(u,v)(x,t)$ satisfying

[TABLE]

[TABLE]

[TABLE]

*For every fixed $t$, the mapping between IC and the solution at time $t$ is continuous on ${\mathcal{H}}$.

There results are valid for all the systems considered in this paper.

2. Analysis of a toy model

In this section, we focus on the qualitative analysis of the following system:

[TABLE]

on the domain $(0,1)$ with NBC. System (2) provides a toy example allowing a cascade of Hopf-bifurcations. It is a simple RD system for which interesting analytical results can be obtained. Note that for $\epsilon=d=1$, and $c$ constant, system (2) results from (1) after a change of variables and taking out the square term in $u$. Here, we focus on asymptotic behavior. We are particularly interested in the co-existence of solutions, toward which convergence depends on initial conditions (IC). To this end, we will first give a detailed analysis of the linearized system of (2) around $(0,0)$. While the analysis might look as straightforward computations, it emphasizes the qualitative dynamics within the ODEs associated with the projection on eigensubspaces and details valid proofs of stability even though zero belongs to the closure of the set of eigenvalues. After that, we provide an analysis of (2). We prove that a global stability result persists for $\alpha<0$. For $\alpha>0$, small, while the $(0,0)$ becomes unstable, we explicit linear subspaces of ${\mathcal{H}}$, for which the solution still evolves towards [math]. We also provide a more general local stability result. Finally, we illustrate some of the dymanics with numerical simulations.

2.1. The linear case

Note first that $(0,0)$ is a constant solution of (2). The linearized system around this point is given by:

[TABLE]

on the domain $(0,1)$ with NBC. The spectral decomposition allows to provide a comprehensive analysis of (3). Classically, we set:

[TABLE]

We recall that the family $(\varphi_{k})_{k\in{\mathbb{N}}}$ is an orthonormal basis of $L^{2}$, and that the functions $\varphi_{k}$ satisfy:

[TABLE]

and

[TABLE]

with

[TABLE]

Looking for solutions expressed as,

[TABLE]

leads by projection on the eigenspace generated by $(\varphi_{k},\varphi_{k})$ to the resolution of the two-dimensional ODE $E_{k}$:

[TABLE]

The eigenvalues of the matrix

[TABLE]

are

[TABLE]

The following proposition summarizes the properties of $\sigma_{k}^{1}$ and $\sigma_{k}^{2}$.

Proposition 1**.**

When $\alpha$ crosses $\lambda_{k}$ from left to right, $\sigma_{k}^{1}$ and $\sigma_{k}^{2}$ cross the imaginary axis from left to right. Furthermore,

[TABLE]

with

[TABLE]

In the two next theorems, we state the main results describing the behavior of (3).

Theorem 2**.**

For $\alpha<0$, for any initial condition $(u(\cdot,0),v(\cdot,0))$ in ${\mathcal{H}}$, we have

[TABLE]

Proof.

For any fixed $N>0$, one can prove that there exists $\delta>0$ such that

[TABLE]

where $(u_{k},v_{k})(t)$ is the solution of $E_{k}$ with

[TABLE]

Note that since $\sigma_{k}^{2}\rightarrow 0$, our computations do not allow us to take $N=\infty$ and conserve the exponential decay. However, for any $\epsilon>0$ there exists $N$ large enough such that

[TABLE]

Since for $N$ large enough and $k>N$ we have,

[TABLE]

the following inequality holds:

[TABLE]

Combining the above results, we can deduce that for any $\epsilon>0$, there exists $T$ such that for $t>T$,

[TABLE]

∎

Theorem 3**.**

*Let $k\in{\mathbb{N}}^{*}$.

For $\alpha=\lambda_{k}$, $(0,0)$ is a center for system $E_{k}$, a source for $E_{l}$ if $l<k$ and a sink for $E_{l}$ if $l>k$. Furthermore, if $u_{l}(0)=v_{l}(0)=0$ for $l\in\{0,...,k-1\}$ then*

[TABLE]

Otherwise,

[TABLE]

For $\lambda_{k}<\alpha<\lambda_{k+1}$, $(0,0)$ is a source for $E_{l}$ si $l\leq k$ and a sink for $E_{l}$ if $l>k$. Furthermore, if $u_{l}(0)=v_{l}(0)=0$ for $l\in\{1,...,k\}$ then

[TABLE]

Otherwise

[TABLE]

Remark 1**.**

It is worth noting that for the linear approximation, the above simple computations allow to characterize eigensubspaces of IC leading to convergence toward the fixed point, periodic solutions, or infinity.

The spectral decomposition of $L^{2}$ has allowed an exhaustive study of the asymptotic behavior. Next, we consider the nonlinear case and show that spectral decomposition remains useful for asymptotic qualitative analysis.

2.2. The Nonlinear case

We now consider the system (2)

[TABLE]

on the domain $(0,1)$ with NBC.

For $\alpha<0$, the fixed point $(0,0)$ is still attracting all the IC in ${\mathcal{H}}$. Indeed, we have:

Theorem 4**.**

For $\alpha<0$, for all IC in ${\mathcal{H}}$

[TABLE]

Proof.

The proof relies on LaSalle’s principle, see [BookTemam, BookCaz1998]. For the reader’s convenience, we give some details of the rigorous proof. We consider arbitrary IC $(u(x,0),v(x,0))\in{\mathcal{H}}$. And consider the solution $(u^{n}=\sum_{k=0}^{n}u^{n}_{k}\varphi_{k},\sum_{k=0}^{n}v^{n}_{k}\varphi_{k}))$ of the approximated problem:

[TABLE]

with IC $u^{n}(x,0)=\sum_{k=0}^{n}\varphi_{k}\int_{0}^{1}u(x,0)\varphi_{k}dx,v^{n}(x,0)=\sum_{k=0}^{n}\varphi_{k}\int_{0}^{1}v(x,0)\varphi_{k}dx$.

The solution of the problem is the limit as $n$ goes to infinity of $(u^{n},v^{n})$. Next, first observe that, for fixed $n$,

[TABLE]

Furthermore,

[TABLE]

which proves that the trajectories are bounded in ${\mathcal{V}}$. Thanks to the compact injection from $H^{1}$ into $L^{2}$, we have compacity allowing us to apply LaSalle’s principle. But here, the $\omega-limit$ set of the trajectory ensued from $(u^{n}(x,0),v^{n}(x,0)$ is reduced to $(0,0)$. It follows that

[TABLE]

in ${\mathcal{H}}$.

Now, fix $\epsilon>0$. And fix $n$ such that

[TABLE]

A simple computation shows that for any $p>n$,

[TABLE]

which leads to, letting $p$ go to infinity to

[TABLE]

Which yields

[TABLE]

∎

trajectories distinct from $(0,0)$. It follows that, for $\alpha>0$, $(0,0)$ becomes unstable for any solution constant in space different from $(0,0)$ will evolve towards the limit cycle of the $ODE$. However, this system allows for the construction of specific solutions of interest.

Theorem 5**.**

For $0<\alpha<\lambda_{1}$, if $u(x)=-u(1-x)$ and $v(x)=-v(1-x)$ then for all IC in ${\mathcal{H}}$

[TABLE]

Proof.

By symmetry, $\int_{0}^{1}u(x,t)dx=\int_{0}^{1}v(x,t)dx=0$. Then, we apply LaSalle’s Principle as in the proof of 4, within the positive invariant subspace defined by $\int_{0}^{1}u(x,t)dx=\int_{0}^{1}v(x,t)dx=0$. Indeed, we have:

[TABLE]

and,

[TABLE]

∎

Remark 2**.**

Note that although $(0,0)$ is unstable for $0<\alpha<\lambda_{1}$, this theorem characterizes a set included in the basin of attraction of $(0,0)$. The arguments provided here apply to larger dimensions. For example, when the space dimension is 2, this set contains IC leading to pattern formation such as spirals (for $\alpha>\alpha_{1}$); see [Amb-2016] where equation (1) was studied theoretically and numerically for $c=0$.

The following lemma proves that a large number of nonlinear terms have a zero integral.

Lemma 1**.**

Let $k,m,n,l\in{\mathbb{N}},k,m,n,l>0$, then

[TABLE]

if and only if

[TABLE]

And

[TABLE]

if and only if one of the subscripts is the sum of the three others or the sum of two of them equals the sum of the other two.

The next theorem provides a local stability result for $0<\alpha<\lambda_{1}.$

Theorem 6**.**

For $0<\alpha<\lambda_{1}$, there exists a sequence $(\mu_{k})_{k\in{\mathbb{N}}}$ such that if

[TABLE]

then

[TABLE]

where $B(0,\mu_{k})\subset{\mathbb{R}}^{2}$ is the ball of center $(0,0)$ and radius $\mu_{k}$.

Proof.

We can write:

[TABLE]

where

[TABLE]

and for $k\geq 1$,

[TABLE]

where

[TABLE]

It follows that the nonlinear terms are bounded by:

[TABLE]

where $C$ is a constant. Thanks to the dynamics of each 2d ODE system $E_{k}$, one can ensure that for any arbitrary small $\epsilon$, for all $k\in{\mathbb{N}}$, there exists $\mu_{k}$ such that if $(u_{k}(0),v_{k}(0))$ belongs to $B(0,\mu_{k})$ the following estimates are valid:

[TABLE]

Thanks to the values of $\lambda_{k}=k^{2}\pi^{2}$, this actually implies that the series $\sum|u_{k}|$ is convergent. Note that, we do not have anymore bounds in $H^{1}$, in this case, and we cannot apply the LaSalle’s principle as previously. We need to construct arguments by hand. Note further that,

[TABLE]

which yields

[TABLE]

Thanks to the bound $|u_{k}|\leq 2\frac{\epsilon}{|\alpha-\lambda_{k}|}$ and the convergence of $\sum\frac{1}{|\alpha-\lambda_{k}|}$ , we obtain for a certain constant $C$:

[TABLE]

Therefore one can choose $\epsilon$ such that:

[TABLE]

as well as

[TABLE]

as long as,

[TABLE]

This implies that $\sum_{k=1}^{+\infty}(u_{k}^{2}+v_{k}^{2})$ converges to a real value as $t$ goes to infinity. Further arguments omitted here allow to prove that the limit is [math].

∎

2.3. Numerical simulations

In this section, we provide a few numerical simulations for system (2). Our simulations have been performed using our own $C^{++}$ program with an $RK4$ numerical scheme. We use a time step of $10^{-5}$ and a space step of $0.02$. We illustrate simulations for values of the parameter $\alpha\in\{1,15\}$. For each value of $\alpha$, we show two pictures : one corresponding to a solution with symmetry which implies that $\int udx=\int vdx=0$ and another one corresponding to a solution for which those integrals are non zero.

In figure $1$, we simulate system (2) for $0<\alpha=1<\lambda_{1}$. IC satisfy the symmetric condition $u(x)=u(1-x)$ and $v(x)=v(1-x)$. ore precisely, or this figure $u(x,0)=v(x,0)=0.5$ on $(0,0.5)$, and $u(x,0)=v(x,0)=-0.5$ on $(0.5,1)$. This implies a symmetry for the solution for all positive times. According to theorem 4, the solution converges toward $(0,0)$ in ${\mathcal{H}}$. The observation of $u$ in panel $A$ illustrates this theoretical result. For $v$, as illustrated in panel B, the evolution is slower. In panel C, we represent two trajectories corresponding to the time evolution of $u$ for two fixed values of $x$. One at $x=0.04$ in red and the other at $x=0.96$ in blue.

In figure $2$, we simulate again system (2) for $0<\alpha=1<\lambda_{1}$. But here, IC do not satisfy the symmetric condition $u(x)=u(1-x)$ and $v(x)=v(1-x)$. In this figure, we set $u(x,0)=v(x,0)=0.5$ on $(0,0.5)$, $u(x,0)=v(x,0)=0$ on $(0.5,1)$. The simulations show that the solution $u$ reaches asymptotically a constant function in space with periodicity in time. There is also periodicity in time for $v$, but the evolution of the shape in space is slower. Asymptotically, since solutions are constant in space, the qualitative behavior is given by the ODE diffusion-less system.

In figure $3$, we simulate system (2) for $\lambda_{1}<\alpha=15<\lambda_{2}$. IC satisfy $u(x)=u(1-x)$ and $v(x)=v(1-x)$ which implies that $\int udx=\int vdx=0i.e.$ $u_{0}(t)=v_{0}(t)=0$. We observe that the solution $u$ evolves non constantly in space with periodicity in time.

In figure $4$, IC do not satisfy the symmetric condition $u(x)=u(1-x)$ and $v(x)=v(1-x)$. We observe that the solution $u$ reaches a constant function in space with periodicity in time. We can note the difference of the amplitude of the limit cycle with the previous simulation; for these IC, $u_{0}(t)$ and $v_{0}(t)$ are no longer zero.

Lastly, in figure 5, we illustrate $u$ as a function of $x$ ant $t$.

3. Analysis and discussion around the Nh-FHN model

3.1. Equation and Assumptions

We now consider the system (1):

[TABLE]

with $f(u)=-u^{3}+3u$, $\epsilon$ small, and Neumann Boundary conditions (NBC). We are interested in the case where the function $c$ take values lesser than $-1$ (excitable property for the diffusion-less system) on some parts of the domain $(a,b)$ and between [math] and $1$ (oscillatory property for the diffusion-less system) on other parts of the domain. More specifically, we consider in this article a function $c$ which is oscillatory in a neighborhood of the center $\frac{a+b}{2}$ and excitable anywhere else; we further assume that excitability depends upon a given parameter $p>0$, see [Amb-2009, Amb-2017]. For sake of simplicity, we assume that the function $c(x)$ (which depends on $p$) is regular and satisfies the following conditions:

[TABLE]

The first assumption 7 implies that the diffusion-less system at $x=0$ has a limit cycle. The second assumption 8 implies that $c$ admits a maximum at $x=0$, furthermore, it implies that the excitability strength decreases from the border toward the center. The third assumption 9 is a compatibility condition with NBC. Finally, assumptions 10-12 relate to the strength of excitability with respect to $p$ for fixed $x$; the higher the value of $|c|$, the weaker the excitation’s strength.

The spectrum analysis of the linearized operator has been carried out in [Amb-2017]. We will rely here on some of these results while following a different presentation. The stationary solution is given by:

[TABLE]

We rewrite system (13) around $(\bar{u},\bar{v})$. We obtain

[TABLE]

We would like to proceed with projection on appropriate subspaces as in previous sections. To that end, we are interested in solutions of the following equation

[TABLE]

Note that equation (15) is a regular Sturm-Liouville problem. We have the classical following spectral theorem, see for example [BookJost2013, BookTrudinger, BookTeschl].

Theorem 7**.**

There exists an increasing sequence of real numbers $(\lambda_{k})$ and an orthogonal basis $(\varphi)_{k\in{\mathbb{N}}}$ of $L^{2}(a,b)$ such that:

[TABLE]

Furthermore,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

(Weyl assymptotics)

The linearized system of (14) writes

[TABLE]

and projection on the $kth$ subspace writes

[TABLE]

The eigenvalues of the matrix associated to $(E_{k})$ are given by

[TABLE]

The following result holds:

Corollary 1**.**

The number of eigenvalues with positive real part is finite while the number of eigenvalues with negative real part is infinite. If we assume that

[TABLE]

then $\sigma_{1}^{0}$ and $\sigma_{2}^{0}$ have a positive real part. And,

[TABLE]

The more elaborated theorem holds, see

Theorem 8**.**

For $p$ small enough, $(0,0)$ is a source for $E_{0}$. For $p$ large enough, $(0,0)$ is a sink . There is an Hopf bifurcation: there exists a value $p^{*}$ for which as $p$ crosses $p*$ from right to left, the real part of $\sigma_{1}^{0}$ and $\sigma_{2}^{0}$ increases from negative to positive. The other eigenvalues remaining with negative real parts.

Remark 3**.**

This theorem relies on the analysis of the Sturm-Liouville problem (15). The first Hopf bifurcation occurs when $\lambda_{0}$ becomes negative. Successive Hopf bifurcations occur when the the successive $\lambda_{k}$ become negative. The number of possible Hopf bifurcation is bounded. The limit-case is given by $c(x)=0$.

3.2. Discussion and Comparison with the Toy model

For the toy model, we provided the following results:

1. (1)

existence of a cascade of Hopf bifurcations 2. (2)

existence of solutions converging toward the unstable stationary point $(0,0)$ for $0<\alpha<\lambda_{1}$ 3. (3)

a local stability result for $0<\alpha<\lambda_{1}$.

The point 1 resulted from the spectral analysis of the linearized operator. The points 2 and 3 were possible thanks to the property that $\int\varphi_{0}\varphi_{k}=0$. Thanks to that property, despite the nonlinearity, the first projection spanned by $\varphi_{0}$ can be controlled in the others subspaces. The symmetry of the reaction term was also necessary for the point 2. Dealing with Nh-FHN, one lose the points 2 and 3. The point 1, i.e. the spectral analysis of the linearized operator remains partially valid and gives good insights. Since an infinite number of eigenvalues is negative, looking at the behavior through the projection into a finite space spanned by the eigenfunctions of the Laplacian provides an interesting way to gain qualitative insights. Beyond that natural perspective, a question that arises is the existence of initial conditions leading to the convergence toward the stationary solution $(c,f(c))$ when the first Hopf bifurcation arises. Such a result was proven for the toy model. These questions are to be addressed in a forthcoming article.

3.3. Numerical simulations

For the numerical simulations, we set $a=-1$, $b=1$, $d=1$ and

[TABLE]

Note that $c(0)=0$ which corresponds to oscillatory behavior for the ODE system, while $c(a)=c(b)=-p<-1$ which corresponds to a stable steady state for the ODE. We simulate equation (1) on $(a,b)$ with an explicit scheme of Runge-Kutta 4 type, with a time step of $10^{-5}$. The value of $\epsilon$ is set to $0.1$. Note that the dynamics depend on the size of the domain and the diffusion coefficient $d$. Here, we illustrate numerical simulations with these parameters fixed. We vary the parameter $p$. We plan to proceed to a further qualitative numerical study in a forthcoming article. Note also that some complementary results have been obtained for a chain with discrete kicks in [Amb-2020-2]. We illustrate here solutions for three values of the parameter $p$. If $p=2.2$, the solution converges toward a stationary solution; see figures 6 and 9. As $p$ is decreased, we observe waves. The amplitude of the waves depend on the parameter $p$, see figures 7,8,9 ($p=2$ and $p=1.9$).

4. Concluding remarks and Perspectives

In this article, we studied a few RD systems with a cubic nonlinearity arising from a Neuroscience context. We provided significant theoretical and numerical qualitative analysis of such infinite-dimensional systems highlighting dynamics arising under the stable regime and near Hopf-Bifurcations. We also emphasized the dependence on initial conditions. As for the nature of the solutions observed, our work relates with previous studies on traveling waves. Our approach is however different by may aspects as this was discussed in introduction. As for the biological motivation, it goes back to the HH RD model and the oscillatory and excitability properties featured in the FitzHugh-Nagumo model. Our first contribution was to consider a non-homogeneous FHN model (Nh-FHN) as simple as possible. We introduced the inhomogeneity with a function $c$ depending upon the space variable $x$. This dependence allows to introduce both excitable and oscillatory properties within the single PDE. Note that in this context, the function $c(x)$ indicates the degree of excitability at the point $x$. Naturally, an immediate application is to provide a model for propagation of waves trough the axon of a single neuron, as this was the case in the original HH model. More generally, Nh-FHN is a model for electrical wave propagation along excitable tissue. It is worth noting that RD excitable systems derivated from FHN type systems have played an important role to model the electrical activity in heart. Spirals, scroll-waves and patterns have been used to characterize arrhythmias; see for example [Maj2011, Pra2015, Prav2021] and references therein cited. There is also a large literature about numerical models used to develop techniques of low voltage fields to induce defibrillation. Those are based on RD ionic current models, see for example [Are2007, Tra2014]. The study of waves, spirals and patterns of electrical cortical activity has also been a very active field of research this recent years, see for example [Mul2018, Mey2022]. In this context, a better knowledge of the underlying mathematical mechanisms from which waves arise appears to be relevant. Here, we illustrated for example how waves may or may not exist depending on the initial conditions. Other simulations not illustrated here show also that the size of the domain an the diffusivity has an impact on the qualitative dynamics of the observed waves. In forthcoming articles, we plan to investigate further these phenomena. This includes the comparison between the solutions observed for the full PDE and finite dimensional solutions, a systematic description of the qualitative dynamics of the observed waves and how the results provided here extend to higher space dimensions.

Supplementary material

The codes used for simulations are available at https://github.com/benjamin-ambrosio/EvEqCoTh2023ToyModel and https://github.com/benjamin-ambrosio/EvEqCoTh2023NhFHN

Acknowledgments

I would like to thank Region Normandie France, the ERDF (European Regional Development Fund) project XTERM, IEA CNRS project IEA00134 for funding.