Propagation of bursting oscillations in coupled Hodgkin-Huxley Reaction-Diffusion systems

TL;DR

This paper develops a mathematical framework for Hodgkin-Huxley reaction-diffusion networks, demonstrating solution properties and analyzing the propagation of bursting oscillations through coupled systems with bifurcation phenomena.

Contribution

It introduces a general mathematical framework for Hodgkin-Huxley reaction-diffusion systems, proving existence, uniqueness, and invariant regions, and explores oscillation propagation and bifurcations.

Findings

Existence and uniqueness of solutions established.

Propagation of bursting oscillations demonstrated.

Bifurcation phenomena observed in coupled systems.

Abstract

We consider networks of reaction-diffusion systems of Hodgkin-Huxley type. We give a general mathematical framework, in which we prove existence and unicity of solutions as well as existence of invariant regions and of the attractor. Then, we illustrate some relevant numerical examples and exhibit bifurcation phenomena and propagation of bursting oscillations through one and two coupled systems.

Equations

$$\left\{\begin{array}[]{rcl}CV_{t}&=&I+\overline{g}_{Na}m^{3}h(E_{Na}-V)+\overline{g}_{K}n^{4}(E_{K}-V)+\overline{g}_{L}(E_{L}-V)\\ n_{t}&=&\alpha_{n}(V)(1-n)-\beta_{n}(V)n\\ m_{t}&=&\alpha_{m}(V)(1-m)-\beta_{m}(V)m\\ h_{t}&=&\alpha_{h}(V)(1-h)-\beta_{h}(V)h,\\ \end{array}\right.$$

$$\begin{array}[]{rcl}\alpha_{n}(V)=0.01\frac{-V+10}{\exp{(1-0.1V)}-1},&&\beta_{n}(V)=0.125\exp(-V/80),\\ \alpha_{m}(V)=0.1\frac{-V+25}{\exp{(2.5-0.1V)}-1},&&\beta_{m}(V)=4\exp(-V/18),\\ \alpha_{h}(V)=0.07\exp(-V/20),&&\beta_{h}(V)=\frac{1}{1+\exp(-0.1V+3))}.\end{array}$$

$$E_K=-12mV,E_{Na}=120mV,E_L=10.6mV$$

$${\overline{g}}_K=36mS/c{m}^2,{\overline{g}}_{Na}=120mS/c{m}^2,{\overline{g}}_L=0.3mS/c{m}^2.$$

$$E_i=\frac{RT}{zF}\ln \frac{[i{]}_{out}}{[i{]}_{in}},$$

$$E_{Na}\simeq 55,$$

$$\left\{\begin{array}[]{rcl}V_{it}&=&dV_{ixx}+I+\overline{g}_{Na}m_{i}^{3}h(E_{Na}-V_{i})+\overline{g}_{K}n_{i}^{4}(E_{K}-V_{i})+\overline{g}_{L}(E_{L}-V_{i})\\ &&+H_{i}(V_{1},...,V_{N}),\qquad i\in\{1,...,N\},\\ n_{it}&=&\alpha_{n}(V_{i})(1-n_{i})-\beta_{n}(V_{i})n_{i}\\ m_{it}&=&\alpha_{m}(V_{i})(1-m_{i})-\beta_{m_{i}}(V)m_{i}\\ h_{it}&=&\alpha_{h}(V_{i})(1-h_{i})-\beta_{h}(V_{i})h_{i}.\\ \end{array}\right.$$

$$H_i(u)=\sum _{j\in \{1,...,N\}}{\alpha }_{ij}(S-u_i)\Gamma (u_j),$$

$$\Gamma (s)=\frac{1}{1+{\exp }^{-\lambda (s-\theta )}},$$

$$S=100,\lambda =20,\theta =60,$$

$$u_t=u_{xx}$$

$$u_x(a)=u_x(b)=0$$

$$\left\{\begin{array}[]{rcl}V_{i}(t)&=&\displaystyle S(t)V_{i}(0)+\int_{0}^{t}S(t-s)(F_{V}(V_{i}(s),n_{i}(s),m_{i}(s),h_{i}(s))\\ &&+H_{i}(V_{1},...,V_{N}))ds\\ n_{i}(t)&=&\displaystyle n_{i}(0)+\int_{0}^{t}F_{n}(V_{i}(s),n_{i}(s))ds\\ m_{i}(t)&=&\displaystyle m_{i}(0)+\int_{0}^{t}F_{m}(V_{i}(s),m_{i}(s))ds\\ h_{i}(t)&=&\displaystyle h_{i}(0)+\int_{0}^{t}F_{h}(V_{i}(s),h_{i}(s))ds,\\ \end{array}\right.$$

$$\begin{array}[]{rcl}F_{V}(V,n,m,h)&=&I+\overline{g}_{Na}m^{3}h(E_{Na}-V)+\overline{g}_{K}n^{4}(E_{K}-V)+\overline{g}_{L}(E_{L}-V),\\ F_{n}(V,n)&=&\alpha_{n}(V)(1-n)-\beta_{n}(V)n,\\ F_{m}(V,n)&=&\alpha_{m}(V)(1-m)-\beta_{m}(V)m,\\ F_{h}(V,n)&=&\alpha_{h}(V)(1-h)-\beta_{h}(V)h.\\ \end{array}$$

$${\Phi }_{U_0}:C([0,T],X^{4N})\to C([0,T],X^{4N})$$

$$||\Phi (U)-U_0||\le KT,$$

$$||\Phi (U_1)-\Phi (U_2)||\le K_{2}T||U_1-U_2||,$$

$$\frac{V(t+h)-V(t)}{h}$$

$$n_t=\alpha (V)(1-n)-\beta (V)n$$

$$A(t)={\int }_0^t\alpha (V)ds,B(t)={\int }_0^t\beta (V)ds.$$

$$n(t)=\exp(-A(t)-B(t))\big{[}n_{0}+\int_{0}^{t}\alpha(V(s))\exp(A(s)+B(s))ds\big{]}$$

$$n(t)\ge 0.$$

$$\begin{array}[]{rcl}n(t)&\leq&\displaystyle\exp(-A(t)-B(t))\big{[}n_{0}+\displaystyle\int_{0}^{t}(\alpha(V)+\beta(V))\exp{(A(s)+B(s)})ds\big{]}\\ &=&\exp({-A(t)-B(t)})\big{[}n_{0}+(\exp{(A(t)+B(t))}-1)ds\big{]}\\ &=&1+(n_{0}-1)\exp(-A(t)-B(t))\\ &\leq&1.\end{array}$$

$$||V_i(t)|{|}_X\le ||V_i(0)|{|}_X+{\int }_0^t(C_1||V_i(s)|{|}_X+C_2)ds$$

$$\Vert |V_i(t)|{|}_X\le ||V_i(0)|{|}_X\exp (C_{1}t)+\frac{C_2}{C_1}.$$

$$V_i(t,x)\in [E_K,E_{Na}]\text{mbox}and{n}_i(t,x),m_i(t,x),h_i(t,x)\in [0,1].$$

$$\omega (U_0)=\bigcap _{t>0}\overline{\bigcup _{s\ge t}T(s)U_0},$$

$$\omega (K)=\bigcap _{t>0}\overline{\bigcup _{s\ge t}T(s)K}.$$

$$\mathcal{A}=\omega (K).$$

$$\underset{t\to +\infty }{\lim }\underset{z\in \mathcal{A}}{\inf }||T(t)U_0-z|{|}_X=0.$$

Full text

1 Introduction

Action potential propagation is a crucial phenomenon for information process in the nervous system and particularly in the brain. One of the paradigmatic models, which has been proposed to describe action potential propagation, is the Hodgkin-Huxley model. Written initially in 1952, see HH , for the description of the electrical activity of the squid giant axon, its formalism has served as basis of number of models widely used in Mathematical Neuroscience, see for example Ambrosio0 ; Ambrosio1 ; Ambrosio2 ; C1 ; Fitz ; Hin ; Mor ; Nagumo and references therein cited. At that time, Hodgkin and Huxley used the new voltage clamp technique to maintain constant the membrane potential. This technique allowed them to elaborate and fit the functional parameters of their model of four ODEs with their experiments. Basically, the model is obtained by considering the cell as an electrical circuit. The membrane acts as a capacitor whereas ionic currents result from ionic channels acting as variable voltage dependent resistances. The model takes into account three ionic currents: potassium ($K^{+}$), sodium ($Na^{+}$) and leakage (mainly chlorure, $Cl^{-}$). The Hodgkin-Huxley (HH) system reads as:

[TABLE]

where subscript $t$ stands for derivative $\frac{d}{dt}$ and where $I$ is the external membrane current, $C$ is the membrane capacitance, $\overline{g}_{i},\ E_{i},i\in\{K,Na,L\}$ are respectively the maximal conductances and the (shifted) Nernst equilibrium potentials. The functions $\alpha(V)$ and $\beta(V)$ describe the transition rates between open and closed states of channels. They read as:

[TABLE]

The (shifted) Nernst equilibrium potentials are given by:

[TABLE]

[TABLE]

These values are taken from Izh , p 37-38, and correspond to those of the Hodgkin-Huxley original paper HH , after a change of variables $V=-V$. Recall that the Nernst equilibrium potentials are obtained by solving, for each ion $i\in\{K,Na,L\}$ the equation:

[TABLE]

where $[i]_{in}$ and $[i]_{out}$ are concentrations of the ions inside and outside the cell. $R=8.315$ is the universal gas constant, $T$ is temperature in Kelvin, $F=96,48$ is the Faraday’s constant, $z$ is the valence of the ion. For example, this computation gives for sodium, with $T=293$, $[i]_{out}=440$, $[i]_{in}=40$ (see Izh p 50):

[TABLE]

which with a shift of $+65$ gives the value of $120$ used here. Originally, Hodgkin and Huxley used the shift to obtain a potential at rest of approximately [math]. Before going into more theoretical aspects, we give some interpretation about the form of conductances. The proportion of open potassium channels is $n^{4}$. This comes from the fact that $4$ opening gates of potassium are required to open the potassium channel. Hence, the $n$ gives the probability of the gate to be in active state and results in the $n^{4}$ term for potassium. For the sodium, it is supposed that there are there are three gates which open the channels and one which close them. Hence, the proportion of sodium opened channels is given by $m^{3}h$, where it is supposed that $m$ stands for the probability of sodium opening gates to be active while $h$ stands for probability of sodium closing gates to be active. For more details on various aspects of the HH model, we refer to Cro ; Erm ; Izh and the original paper HH . In the present paper, we consider a general network of reaction-diffusion (RD) HH equations. Note that the study of networks of reaction-diffusion systems start to attract an increasing number of mathematical analysis, see for example, Ambrosio1 ; Ambrosio2 ; Yang . They appear naturally in the neuroscience context. The general network reads as:

[TABLE]

in a bounded domain $\Omega=(a,b)\subset\mathbb{R}$ and with Neumann boundary conditions $V_{ix}(a)=V_{ix}(b)=0$ and where,

[TABLE]

[TABLE]

with

[TABLE]

and $\alpha_{ij}\geq 0$.

The functions $H_{i}$ represent excitatory nonlinear coupling (i.e. chemical synaptic coupling) between neurons, see Ambrosio3 ; Ambrosio4 ; Bel ; C1 and references therein cited. Without loss of generality, we have set the constant $C$ equal to $1$. Hence, the value $S=100$ allows the neuron $i$ to receive an excitatory input each time one of his presynaptic neurons crosses the treshold value $\theta=60$. The value $\lambda=20$, implies that $\Gamma$ approximates the Heaviside function. The Neumann condition is chosen here to not induce boundary effects. Here, $I$ and $\alpha^{\prime}s$ are regular non negative bounded function of $x$. Our first aim is to provide a mathematical framework for solutions of system (3). This is done in section 2. We prove the global existence of solutions as well as the existence of invariant regions and the existence of an attractor. Then, in the third section, we deal with numerical simulations for a single non-homogeneous HH RD system and two coupled HH RD systems. In this last part, we exhibit some relevant bifurcation phenomena and propagation of bursting oscillations from first neuron to the second neuron. Note that system (3) is autonomous, which means that the bursting oscillations are generated here without external driving component as in Ambrosio0 . See also Cle .

2 Mathematical Framework

2.1 Existence of the solution

Let $X=C(\bar{\Omega})=C([a,b])$ the space of continuous functions defined on the real interval $[a,b]$. We start with the proof of the existence of the solution of (3). The proof is divided into four parts. First, we prove the local existence and uniqueness of a mild-solution, thanks to a fixed point theorem. Then, we prove the regularity of this mild solution. Finally, we prove the bounds for $n_{i},m_{i},h_{i}$ as well as the global existence.

Theorem 2.1

For all initial conditions $(V_{i}(0),n_{i}(0),m_{i}(0),h_{i}(0))\in X^{4N}$ such that $\forall x\in[a,b]$, $n_{i}(0,x),m_{i}(0,x),h_{i}(0,x)\in[0,1]$, there exists a unique solution $U\in C([0,+\infty[,X^{4N})$ of (3). Furthermore, $U\in C^{1}(]0,+\infty[,X^{N})$ and for all $t>0$, $n_{i}(t,x),m_{i}(t,x),h_{i}(t,x)\in[0,1]$.

Proof

It is known that the equation

[TABLE]

with

[TABLE]

generates an analytical semigroup on $X$, see Lunardi , Theorem 3.1.22 p 98. Let us denote by $S(t)_{t\geq 0}$ this semigroup. The first step is to look for mild solutions of (3). Mild-solutions are solutions of the variation of constants formula of (3), and written with the semigroup notation. The general network equation for mild-solutions reads as:

[TABLE]

where:

[TABLE]

The local existence follows by application of the fixed point theorem, to the function:

[TABLE]

defined by the right-hand side of (6), in which we fix the initial condition $U_{0}\in X^{4N}$ (below, for sake of simplicity, we drop the subscript $U_{0}$ on $\phi$). Let $U$ denote the function $(V_{i},n_{i},m_{i},h_{i})(t),i\in\{1,...,N\},t\in[0,T]$. For all $U\in B(U_{0},r)$, the ball of center $U_{0}$ and radius $r$ in $C([0,T],X^{4N})$, we have, for fixed $T>0$

[TABLE]

where $||\cdot||$ stands here for the sup norm on $C([0,T],X^{4N})$ and where $K$ is a constant depending on $r$ which comes from the boundedness theorem. We have also used the property of the linear semigroup $S(t)$: $||S(t)u||_{X}\leq||u||_{X}$ and $\lim_{t\rightarrow 0}S(t)u_{0}=u_{0}$, see Lunardi p 35. It follows that, for $T<\frac{r}{K}$, $\Phi(U)\in B(U_{0},r)$. By analog computations, we prove that

[TABLE]

which proves that for $T<\min(\frac{1}{K_{1}},\frac{1}{K_{2}})$, $\phi$ is a contraction mapping. Therefore, $\Phi$ has a unique fixed point which provides the existence and uniqueness of the local mild solution.

Now, we deal with the regularity. We already know, that, solutions belong to $C([0,T],X^{4N})$. Since the $H$ term, does not complicate the proof, without loss of generality, we drop the subscript $i$. A computation shows that the quantity

[TABLE]

is well defined, for $t>0$, as $h\rightarrow 0$ and that the resulting function is contiuous on $X$. For $t=0$, the derivative may not exist, since $\frac{S(h)V_{0}-V_{0}}{h}$ admits a limit as $h\rightarrow 0$ if and only if $v_{0}\in C^{\prime\prime}([a,b])$. For functions $m,n$ and $h$ the same computation shows that they belong to $C^{1}([0,T],X)$. Hence, the regularity statement follows.

Now, we deal with the global existence. First, we prove that $n_{i},m_{i}$ and $h_{i}$ remain in $[0,1]$. Without loss of generality, we drop the $i$ subscript and consider only the $n$ variable.

It satisfies,

[TABLE]

Let

[TABLE]

Then,

[TABLE]

with $0\leq\alpha(V)$ and $0\leq\beta(V)$.

It follows that:

[TABLE]

Also,

[TABLE]

Then, the global existence follows from the inequality,

[TABLE]

which comes from the first equation of (3) and which by Gronwall inequality yields:

[TABLE]

2.2 Existence of an invariant region

We assume that the values of parameters are as specified in (2) and (5). This implies in particular that are $E_{K}<E_{L}<S<E_{Na}$. We also assume that $\bar{g}_{L}(E_{Na}-E_{L})>\sup_{x\in[a,b]}I(x)$111We can drop this last asumption. In this case, the upper bound $E_{Na}$ of the invariant region for $V_{i}$ is replaced by $\frac{1}{\bar{g}_{L}}\sup_{x\in[a,b]}I(x)+E_{L}$..

Theorem 2.2

We assume that $\forall i\in\{1,...,N\}$, $\forall x\in[a,b]$ $V_{i}(0,x)\in[E_{K},E_{Na}]$ and $n_{i}(0,x),m_{i}(0,x),h_{i}(0,x)\in[0,1]$, then $\forall t>0$ and $\forall x\in(a,b)$,

[TABLE]

Proof

We have already proved that $n_{i}(t,x),m_{i}(t,x),h_{i}(t,x)\in[0,1]$. To prove $V_{i}(t,x)\in[E_{K},E_{Na}]$ we use the following argument. Let $t_{0}$ the first time that $V_{i}$ reaches the value $E_{K}$ at the point $x_{0}$, then $\frac{\partial^{2}V_{i}}{\partial x^{2}}\geq 0$ which implies $V_{it}(t_{0},x_{0})>0$ since all the other terms are non negative and $\bar{g}_{L}(E_{L}-E_{K})>0$. Analogously, let $t_{0}$ the first time that $V_{i}$ reaches the value $E_{Na}$ at the point $x_{0}$ then $\frac{\partial^{2}V_{i}}{\partial x^{2}}\leq 0$ which implies $V_{it}(t_{0},x_{0})<0$ since all the other terms are non positive and $\bar{g}_{L}(E_{L}-E_{Na})<0$.

2.3 Existence of an attractor

We denote by $K=C([a,b],([E_{K},E_{Na}]\times[0,1]^{3})^{N})$. In this section, we use the framework of Tem to establish the existence of an invariant, compact and connected set $\mathcal{A}\subset K$ which attracts all tre trajectories starting in $K$. We first recall the definition of the $\omega-$limit set. We denote by $(T(t))_{t\geq 0}$ the semi-group associated with equation (3).

Definition 1

For $U_{0}\in K$, we define the $\omega-limit$ set of $U_{0}$ by:

[TABLE]

and similarly the $\omega-limit$ set of $K$ by:

[TABLE]

The following theorem holds.

Theorem 2.3

Let

[TABLE]

Then $\mathcal{A}$ is a nonempty invariant, compact and connected set in $X^{4N}$. Furthermore, for all initial condition $U_{0}$ in $K$, the solution $T(t)U_{0}$ of (12) starting at $U_{0}$ verifies:

[TABLE]

Proof

We use the framework of Tem . Hence, we need to establish some compacity on trajectories. For convenience, we drop the subscript $i$. For variable $V$, without loss of generality, we write

[TABLE]

since the coupling term $H$ does not change the proof arguments. We prove that $\overline{\cup_{t\geq t_{0}}T_{1}(t)K}$, is compact in $X$, for some $t_{0}>0$. By Ascoli’s theorem, it is sufficient to prove the equicontinuity of

[TABLE]

for some $t_{0}>0$. Then, we deal with variable $n$. For, variables $m,h$ the same arguments apply. We prove that:

[TABLE]

where

[TABLE]

and

[TABLE]

According to Tem , this suffices to prove theorem 2.3. Let us start with $V$. Recall that for all $t>0$, $V(t)\in C^{\infty}([a,b])$. Now, we establish a bound on $V_{x}$. For all $u\in C^{2}([a,b])$ it holds,

[TABLE]

since we assume $u^{\prime}(a)=0$. It follows that

[TABLE]

and that:

[TABLE]

with $C=(1+(b-a)^{p+1})^{\frac{1}{p}}$, where $||\cdot||_{k,p}$ denotes the norm in the Sobolev space $W^{k,p}$ and $||\cdot||_{p}$ denotes the norm in the Lebesgue space $L^{p}$. Note that inequality (9), is actually a particular case of Agmon-Douglis-Niremberg estimates, see Lunardi p 72 and AgDoNi . However, since we consider here the 1-dimensional case the proof is much more simpler. Next, we adapt some arguments from Har . First, note that the following classical inequalities hold:

[TABLE]

and

[TABLE]

for some generic constants $C_{p}$, see Lunardi p 35. This implies that:

[TABLE]

Now, let $\theta\in]0,1[$, $p$ and $r$ such that:

[TABLE]

for example, we can choose $\theta=\frac{3}{4}$, $p>3$ and $r>3$. Then, we have ,see Har and also Henry ,

[TABLE]

Finally, we obtain that for all $t\in]0,1]$:

[TABLE]

Applying this result in (8) implies:

[TABLE]

where $C_{2}$ comes from the fact that $V,m,n,h$ remain in $[E_{K},E_{Na}]\times[0,1]^{3}$. This implies that

[TABLE]

where the constant $C$ depends on $K$ but not on $t_{0}$, neither in $U_{0}$. This proves the equicontinuity of

[TABLE]

Now, we deal with $n$. We compute,

[TABLE]

we obtain (for convenience, we only write the dependence on $x$ and $y$),

[TABLE]

which by integration leads to

[TABLE]

where $\xi=\inf_{y\in[a,b]}(\alpha(y)+\beta(y))>0$ and $C$ is a generic constant. This gives the announced statement.

Remark 1

The difficult point in the proof was to obtain the bound (11).

3 Bifurcations and propagation of oscillations in one and two coupled non-homogeneous neurons

3.1 Bifurcation in a single non-homogeneous neuron

In this section, we deal with equation (3) with $N=1$. It reads:

[TABLE]

We have proceeded to numerical simulations in the domain $\Omega=(a,b)=(0,100),d=1$, on time interval $[0,500]$, using our own C++ program with a uniform finite difference scheme in space and a Runge-Kutta method integration in time. Hence, the space-mesh reads, $(ih)_{i\in\{0,100\}},h=1$, and the time-mesh reads, $(n\Delta t)_{n\in\{0,50000\}},\Delta t=0.01.$ Note that since the equation for the potential is linear in $V$ and the equation for channels are linear in $n,m,h$. One can, from step to step, compute $n(t+\Delta t),m(t+\Delta t),h(t+\Delta t)$ for fixed $V(t)$. And use a similar argument for $V$ computation. A discussion on this aspects can be found in Bal . The particularity is that we choose a regular non homogeneous $I$, approximating the discontinuous function $\tilde{I}$ defined by:

[TABLE]

In the following descriptions, for simplicity we identify $\tilde{I}(x)$ with $I(x)$. We will use the same convention for $\alpha_{21}$ on section 3.2.

We emphasize here three different regimes corresponding to different region of parameter value $I_{0}$.

Stationary and periodic solutions for $I_{0}\simeq 5.2$

**Convergence toward stationary solution for $I_{0}=5.2$

**For several initial conditions and for small $I_{0}$, i.e. for $0\leq I_{0}\leq 5.2$, we obtain asymptotically a convergence toward a steady state. In figure 1, we have illustrated this phenomenon for $I_{0}=5.2$, for initial condition $(1,1,1,1)$. The first panel illustrates the stationary space dependent solution, which is attained asymptotically. The second panel illustrates the time evolution for $x=0$ and $x=100$.

**Convergence toward periodic solution for $I_{0}=5.3$

**For the same initial conditions as before and for $I_{0}=5.3$, we obtain asymptotically a convergence toward a periodic solution. This periodic solution consists of large solutions starting at $x=0$ and propagating form left to right. We have illustrated this phenomenon in figure 2, for initial condition $(1,1,1,1)$. The first panel illustrates the propagation of oscillations from left to right. In fact, it shows the solution along the space for a fixed time $t=500$. As in the previous figures, the second panel illustrates the time evolution for $x=0$ and $x=100$. It shows periodic solutions.

**Convergence toward stationary solution for $I_{0}=5.3$

**We note for this value the coexistence of attractive stationary and periodic solutions. Indeed, for $I_{0}=5.3$, and initial condition $(0,0,0,0)$, we observe a convergence toward a stationary solution. The phenomenon is depicted as before in figure 3.

Propagation of busting oscillations

**Convergence toward a solution with propagation of busting oscillations for $I_{0}=130$

**For $I_{0}=130$, we observe propagation of busting oscillations. The left cells oscillate continuously, but not all the oscillations propagate trough the right-cells. In fact, the right-cells alternate quiescent phases with oscillating phases, which is a bursting phenomenon. We illustrate this phenomena in figure 4. The left-top panel illustrates the time evolution at $x=0$. We observe oscillations. The right-top panel illustrates the time evolution at $x=100$. We observe bursting phenomena. The left-middle panel illustrates the time evolution at $x=8$. We observe the so-called mixed-mode oscillations: an alternation of small and large oscillations. The right-middle panel illustrates the time evolution at $x=8$, in the (V,n,m) variables. The left-down panel illustrates the propagation of oscillations in space at fixed time $t=200$. The right-down panel illustrate the quiescent phase in space at fixed time $t=250$.

Death-spot

**Convergence toward a solution with death-spot $I_{0}=145$

**For $I_{0}=145$, we observe this another interesting phenomenon. The left cells oscillate continuously, but the wave propagation fails to propagate along the whole space. This is known as the death-spot phenomena, see Ambrosio0 ; Yan1 ; Yan2 and references therein cited. The right cells seem to have reached a stationary state whereas the left cells are still oscillating. The left-top panel illustrates the time evolution at $x=0$. We observe oscillations. The right-top panel illustrates the time evolution at $x=100$. We observe a near-stationary solution. The left-middle panel illustrates the time evolution at $x=8$. We observe small oscillations. The right-middle panel illustrates the time evolution at $x=8$, in the (V,n,m) variables. The left-down panel illustrates the oscillations of left cells and the near stationary state for right cells in space at fixed time $t=200$ and time $t=250$.

3.2 Propagation of bursting oscillations in 2 coupled neurons

In this section, we deal with equation (3) with $N=2$, $I_{2}=0$, $\alpha_{ij}=0$ for $(i,j)\neq(2,1)$ and,

[TABLE]

It reads:

[TABLE]

We have proceeded to numerical simulations as before with

[TABLE]

As the neuron $1$ is not affected by the coupling, the dynamic of this neuron has already been described. For the second neuron, we observe the propagation of bursting oscillations from neuron $1$ to neuron $2$. As the coupling acts only on the right side of the neuron $2$, oscillations propagate from right to left. It is interesting to note that there is also a short propagation toward the right boundary. We can observe this on the spike with a double head at the right side on the right-down panel of figure 6.

4 Conclusion

In this paper, we have considered a general network of HH RD systems. We have proved existence and uniqueness of solutions as well as the existence of invariant region and of the attractor in the space of continuous functions. We have also exhibited bifurcation phenomena and propagation of bursting oscillations along one and two coupled neurons. These numerical results show that a rich behavior may be obtained thanks to the space-inhomogeneities. In future works, we aim to develop theoretical tools to describe and analyze these phenomena.

Acknowledgments

The authors would like to thank Région Normandie and FEDER XTERM for financial support.