## Abstract

We investigate a system of partial differential equations of reaction–diffusion type which displays propagation of bursting oscillations. This system represents the time evolution of an assembly of cells constituted by a small nucleus of bursting cells near the origin immersed in the middle of excitable cells. We show that this system displays a global attractor in an appropriated functional space. Numerical simulations show the existence in this attractor of recurrent solutions which are waves propagating from the central source. The propagation seems possible if the excitability of the neighbouring cells is above some threshold.

## 1. Introduction

We introduce in this article a system that exhibits propagation patterns of bursting oscillations. Previous motivations for such mathematical developments appeared in the literature of physiological models. For instance, in Aslanidi *et al*. (2001), excitation wave propagation was proposed as a possible mechanism for signal transmission in pancreatic islets of Langerhans. There are certainly different possible biophysical scenarios to generate such propagations. We focus here on one possible mechanism that seems rather mathematically tractable. This article includes numerical simulations and summarizes the mathematical tools developed to analyse the system. When a system is spatially extended, it is known that a localized stimulus may generate stable propagating pulses. However, the effect of periodic stimuli incorporating the spatially extended system has been studied much less. For a review on travelling waves in excitable media, see, for instance, Tyson & Keener (1988). For a general introduction to emergence of coherent structure, see Scott (1999).

## 2. On an ordinary differential equation of FitzHugh-Nagumo type

We focus here on the differential equations 2.1

Parameter *ϵ* is small and the dynamics is of fast–slow type. The function *f* is a cubic whose precise shape is irrelevant. We fix conveniently *f*(*u*)=−*u*^{3}+3*u* so that, with this scaling, *f*(0)=0 and *f* displays *u*=−1 and *u*=1 as critical points. Parameter *δ* is also small and it will appear later to be quite relevant to the mathematical discussion. We assume that *δ* is small enough so that the cubic *v*=*f*(*u*) intersects the nullcline *u*=*c*+*δ**v* in a single point. The parameter *c* is assumed to vary and accordingly system equation (2.1) undergoes bifurcations. Assume *δ*=0, then, if −1<*c*<1, equation (2.1) displays an attractive limit cycle. If *c*<−1, equation (2.1) displays a single stable stationary point which is globally attractive. A remarkable fact is the existence of a threshold such that, if the initial data are above the threshold, the corresponding orbit goes straight to the stationary point; if instead the initial data are below the threshold, the orbit undergoes a long incursion in the phase plane before returning back to the stable stationary point. This system is a paradigmatic model for the so-called excitability property often observed with the action potential of electrophysiology. As the parameter *c* varies, the value *c*=−1 is a bifurcation value; the system displays a Hopf bifurcation. The limit cycle emerges very quickly with a large size in relation to the fast–slow nature of the differential system. Such type of Hopf bifurcation is usually called singular as it vindicates the singular perturbation context of the parameter *ϵ*. This dynamical system is often used to describe the electrical activity of cells (such as cardiomyocytes, neurons, etc.). In such a modelling context, the electrical activity of the cell is either a pacemaker (if −1<*c*<1) or excitable (if |*c*|>1).

We consider finally a situation of the so-called dynamical bifurcation where the bifurcating parameter *c* is replaced by a slowly varying function of the time. For instance, this yields
2.2
with
moving back and forth through the critical value −1. The generic solution of the system displays alternatively pulsatile phases separated by quiescent phases. This is characteristic of the so-called bursting oscillations. These complex oscillations are ubiquitous in physiology and can be currently observed in neurophysiology, cardiophysiology, hormone secretions, etc. System (2.2) is one of the simplest that generates this type of oscillations (e.g. Françoise 2005).

Assume *δ*>0. We observe the following important fact:
2.3
Note then that
Note furthermore that
This yields
2.4
If we set
by integration, we deduce
This shows the existence of a bounded absorbing set and that equation (2.2) generates a global dynamical system. Note the importance of the term *δ* in bounds.

## 3. On a semi-linear partial differential equation of reaction–diffusion type

We now consider the partial differential equation (PDE) system
3.1
with
which is of reaction–diffusion type with a spatial heterogeneity
The domain *Ω* is open and bounded with a Lipschitz boundary ∂*Ω*. This system can be thought of as representing a set of cells constituted by a small nucleus of pacemakers near the origin immersed among an assembly of excitable cells. The parameter *c*_{0} defines the excitability of the neighbouring cells. The excitability of the neighbouring cells is higher and higher as the parameter *c*_{0} is close to −1. We denote the initial data (*u*(0,*x*),*v*(0,*x*)) and assume boundary conditions of Neumann type on ∂*Ω*. There are several different functional spaces that are interesting to develop the analysis of the equation. We consider, for instance, the space *Y* =*L*^{2}(*Ω*)×*L*^{2}(*Ω*) or the Sobolev space *X*=*H*^{1}(*Ω*)×*H*^{1}(*Ω*).

The next observation is that essentially the same computation previously done for the ordinary differential equation extends to the PDE. Consider more precisely 3.2 Note that the result above is actually true both for Neumann and for Dirichlet boundary conditions. Denote equation (3.2) yields hence

This shows the existence of an absorbing set in *Y* .

The next step is that actually the orbit is bounded in *X*=*H*^{1}(*Ω*)×*H*^{1}(*Ω*). This can be proved by multiplying the first equation by −Δ*u*, the second equation by −Δ*v*, integrating over *Ω* and summing up the two equations. The calculus used to obtain the bounds are a bit too involved to be included here, although they can be adapted from Marion (1989) and Temam (1988); see also Ambrosio (2009). By the Rellich–Kondrachov theorem, the orbit is thus precompact in *Y* . This shows that the *ω*-limit set of the orbit is compact (and connected) and defines an attractor. On the *ω*-limit set, the solution defines a semigroup. Classical theorems of dynamical systems recurrence apply. By Birkhoff’s theory, a minimal invariant set is a union of recurrent orbits.

## 4. On a system of partial differential equations of reaction–diffusion type with a forcing term

We now consider the PDE system 4.1 with which is of reaction–diffusion type with a spatial heterogeneity and a forcing term

The same inequalities are in order and again allow us to conclude that the orbit is precompact in *Y* . The single difference is that, as the system is no longer autonomous, the corresponding dynamical system is a process in the sense of Daffermos. It is, as usual, possible to associate a semigroup on the product (Haraux 1991). We conclude in both cases the existence of an attractor. We now proceed to numerical simulations. They tend to show that the attractor is not discrete and that it contains effectively recurrent trajectories. Some of these trajectories are perhaps periodic, although we cannot determine this yet. These recurrent trajectories are associated with a wave propagation. Our approach should be compared with previous contributions. In Yanagita *et al*. (2005, 2006), the authors considered a one-dimensional chain of coupled FitzHugh–Nagumo systems of the type
where *I* is a periodic input stimulus added to the first element and there are periodic boundary conditions. The authors also discuss the existence of a propagating pulse in terms of *ϵ*. Note that a simple change of variable *v*_{i}↦*v*_{i}−*I* yields the system
whose continuous limit displays
Hence, this PDE system is quite similar to the one studied here except that there is a single diffusion on *u* and that it is one-dimensional. We expect little change due to the diffusion term in *v* and accordingly our system may be seen as a two-dimensional version of theirs. As they explore the dependency in *ϵ*, we have focused here on the dependency on *c*(*x*,*t*) in two ways: the size of the central nucleus and the excitability of the neighbouring cells (*c*_{0}). Our interpretation in terms of propagation of bursting oscillations was not evidenced previously.

## 5. Numerical simulations of the two-dimensional case

Numerical simulations are based on a discretization and a finite difference method with five points in space and an explicit Euler scheme in time. We do not focus on the weak diffusion limit and set *α*=*β*=1. The domain is a square with a grid of 100×100 squares. Initial data are chosen around *u*(0,*x*)=−1.5,*v*(0,*x*)=0.1. The solution at time *t*_{k} is represented at each point of the grid as a vector . Given the solution at time *t*_{k}, the solution at time *t*_{k+1} is given by
5.1
5.2
where the operator *δ*_{ij} is defined as
for the points that are not on the boundary of the domain and, for instance, *δ*_{10} and *δ*_{00} are defined by

The operator *δ*_{ij} is defined similarly for other points of the boundary, and probably with *v*_{kh} instead of *u*_{kh}.

The two components of the reaction are *F*_{1}(*u*,*v*)=*f*(*u*)−*v*, and *F*_{2}(*u*,*v*,*i*,*j*)=*u*−*c*_{ij}(*t*_{k}) (*c*_{ij}(*t*_{k}) is the vector that represents *c*(*x*,*t*) on the grid at time *t*_{k}).

Parameter *ϵ* is fixed around *ϵ*=0.1; parameter *δ* has been varied around 0.01 and 0.1. The figures below are the result of the simulation for *δ*=0.01. The central nucleus *B*(0,*r*) is represented by four cells in the case of *δ*=0.01.

### (a) Numerical simulations of system (3.1)

The numerical simulations shows a propagation if the parameter *c*_{0}, which characterizes the excitability of the neighbouring cells, is above some threshold. More precisely, we include the result of the numerical simulations for the following fixed values of *c*_{0}:

*c*_{0}=−1.3 (figure 1)The solutions evolve to a stationary solution.

*c*_{0}=−1.195 (figure 2)The solutions still evolve to a stationary solution after several oscillations of small amplitude.

*c*_{0}=−1.19302 (figures 3 and 4)Above this value of

*c*_{0}, the solution is no longer asymptotically stationary. After a rather long time, there is a propagation of a wave. The time of appearance of this propagation depends on the initial condition.*c*_{0}=−1.19 (figure 5)After some oscillations of small amplitude, there is a propagation of a wave.

*c*_{0}=−1.15 (figure 6)There are no more transitory oscillations of small amplitude and there is a propagation of a wave.

Similar numerical simulations have been done with *δ*=0.1. For the same size of the central nucleus, for all values of *c*<−1 the solutions evolve asymptotically to a stationary solution. But if the size of the central nucleus (*B*(0,*r*)) is increased then there is again a propagation. This is related to the critical size of a pacemaker, which is discussed in Keener & Sneyd (1998, ch. 14).

### (b) Numerical simulations of system (4.1) and propagation of the bursting oscillations

The same result is observed numerically: there is a threshold such that, if the excitability *c*_{0} is below the threshold, the solution evolves to a stationary solution. If it is above, there is a propagation of the bursting oscillations. The number of spikes increases as *γ*→0 and/or as *c*_{0}→−1. Initial conditions are around the values *u*(0,*x*)=−1,*v*(0,*x*)=0.1, *δ*=0.01; propagation of bursting oscillations is seen in figure 7.

## 6. Conclusion

In conclusion, the numerical simulations suggest that the attractor whose existence is proved mathematically contains non-trivial recurrent solutions associated with a propagation phenomenon. The discretization techniques also suggest that a similar situation occurs for a large system of linearly coupled forced FitzHugh–Nagumo systems. We find previous interpretations of similar systems in terms of an external periodic forcing acting on a central nucleus in Yanagita *et al*. (2005, 2006). In this regard, our system proposes a natural extension to a two-dimensional system (with both diffusions on the potential and on the recovery variable). This relates to physiological experiments made in cardiodynamics (Chialvo & Jalife 1987; Michaels *et al*. 1989) in which a central nucleus of pacemaker is vagally driven. Our interpretation in terms of propagating bursting oscillations seems novel.

## Acknowledgements

This research was supported by ANR grant ‘Analyse non linéaire et Rythmes du vivant’.

## Footnotes

One contribution of 17 to a Theme Issue ‘From biological and clinical experiments to mathematical models’.

- © 2009 The Royal Society