# Pattern selection in a lattice of pulse-coupled oscillators

## I Introduction

The study of the collective behavior of populations of interacting nonlinear oscillators has attracted the interest of physicists and mathematicians for many years since they can be used to modelize several chemical, biological and physical systems. Among them, we should mention cardiac pacemakers cells, integrate and fire neurons and other systems made of excitable units. Most of the theoretical papers that have appeared in the scientific literature deal with oscillators interacting through continuous-time couplings, allowing them to describe the system by means of coupled differential equations and apply most of the modern nonlinear dynamics techniques. More challenging from a theoretical point of view is to consider a pulse-coupling, or in other words, oscillators coupled through instantaneous interactions that take place in a very specific moment of its period. The richness of behavior of these pulse-coupled oscillatory systems includes synchronization phenomena, spatio-temporal pattern formation (we could mention, for instance, traveling waves, chessboard structures, and periodic waves ), rhythm anihilation, self-organized criticality,…

Most of the work on pattern formation has been done in mean-field models or populations of just a few oscillators. However, such restrictions do not allow to consider the effect of certain variables whose effect can be crucial for realistic systems. The specific topology of the connections or geometry of the system are some typical examples which usually induce important changes in the collective behavior of these models. Pattern formation usually takes place when oscillatory units interact in an inhibitory way, although it has also been shown that the shape of the interacting pulse, when the spike lasts for a certain amount of time, or time delays in the interactions can lead to spatio-temporal pattern formation also in the case of excitatory couplings. Only recently, general solutions for the general case, where the patterns existence and stability is proved, have been worked out. The aim of this paper is to study some pattern properties and get a quantitative estimation of the probability of pattern selection under arbitrary initial conditions or, in the language of dynamical systems, the volume of the basin of attraction of each pattern. Keeping this goal in mind, we will use the general results given in  where assuming a system defined on a ring the authors developed a mathematical formalism powerful enough to get analytic information of the system. Not only about the mechanisms which are responsible for synchronization and formation of spatio-temporal structures, but also, as a complement, to proof under which conditions they are stable solutions of the dynamical equations.

Despite the apparent simplicity of the model, some ring lattices of pulse-coupled oscillators are currently used to modelize certain types of cardiac arhythmias where there is an abnormally rapid heartbeat whose period is set by the time that an excitation takes to travel the circuit . Moreover, there are experiments where rings of a few R15 neurons from Aplysia are constructed and stable patterns are reported . Our 1d model allows us to study analytically the most simple patterns and understand their mechanisms of selection.

The structure of this paper is as follows. In Sec II we review the model introduced in  as well as set the notation used throughout the paper. In Section III we study some pattern properties which will be useful for, in Section IV, propose an estimation of the probability of selection of each pattern. In the last section we present our conclusions.

## II The model

Our system consists in a ring of $`(N+1)`$ pulse-coupled oscillators. The phase of each oscillator $`\varphi _i`$ evolves linearly in time

$$\frac{d\varphi _i}{dt}=1i=0,\mathrm{},N$$
(1)   
until one of them reaches the threshold value $`\varphi _{th}=1`$. When this happens the oscillator fires and changes the state of its rightmost neighbor according to

$`\varphi _i1\{\begin{array}{c}\varphi _i0\hfill \\ \varphi _{i+1}\varphi _{i+1}+\epsilon \varphi _{i+1}\mu \varphi _{i+1}\hfill \end{array}`$

subjected to periodic boundary conditions, i.e. $`N+10`$, and where $`\epsilon `$ denotes the strength of the coupling and $`\mu =1+\epsilon `$. Where we have assumed that, from an effective point of view, the pulse-interaction between oscillators, as well as the state of each unit of the system, can be described in terms of changes in the phase, or in other words, in terms of the so called phase response curve (PRC), $`\epsilon \varphi `$ in our case. A PRC for a given oscillator represents the phase advance or delay as a result of receiving an external stimuli (the pulse) at different moments in the cycle of the oscillator. We will assume $`\epsilon <0`$ througout the paper, as we are only interested in spatio-temporal pattern formation and $`\epsilon >0`$ always leads to the globally synchronized state. This linear PRC has physical sense in some situations. For instance, it shows up when we expand the non-linear PRC for the Peskin model of pacemaker cardiac cells  in powers of the convexity of the driving or in neuronal modelling. In practice, however, this condition can be relaxed since a nonlinear PRC does not change the qualitative behavior of the model provided the number of fixed points of the dynamics is not altered. Moreover, a linear PRC has the advantage of making the system tractable from an analytical point of view.

Let us describe the notation used in the paper. The population is ordered according to the following criterion: The oscillator which fires will be always labeled as unit 0 and the rest of the population will be ordered from this unit clockwise. After the firing, the system is driven until another oscillator reaches the threshold. Then, we relabel the units such that the oscillator at $`\varphi =1`$ is now unit number 0, and so on. This firing + driving (FD) process for $`N+1`$ oscillators can be described through a suitable transformation

$$\stackrel{}{\varphi }^{}=T_k(\stackrel{}{\varphi })\stackrel{}{1}+\text{𝕄}_k\stackrel{}{\varphi },$$
(2)   
where $`\text{𝕄}_k`$ is a $`N\times N`$ matrix, $`\stackrel{}{\varphi }`$ is a vector with $`N`$ components, $`\stackrel{}{1}`$ is a vector with all its components equal to one and $`k`$ stands for the index of the oscillator which will fire next. We call this kind of transformation a firing map, and we have to define as many firing maps as oscillators could fire, that is, index $`k`$ must run from $`k=1`$ ($`\varphi _1`$ fires) to $`N`$ ($`\varphi _N`$ fires). For example, in the $`N+1=4`$ oscillators case we have that the firing map corresponding to the FD process where $`\varphi _2`$ is the next oscillator which do fire,

| $`\varphi _0=1`$ | $`\stackrel{\text{firing}}{}`$ | $`0`$ | $`\stackrel{\text{driving}}{}`$ | $`1\varphi _2=\varphi _2^{}`$ |
| --- | --- | --- | --- | --- |
| $`\varphi _1`$ | $``$ | $`\mu \varphi _1`$ | $``$ | $`\mu \varphi _1+1\varphi _2=\varphi _3^{}`$ |
| $`\varphi _2`$ | $``$ | $`\varphi _2`$ | $``$ | $`1=\varphi _0^{}`$ |
| $`\varphi _3`$ | $``$ | $`\varphi _3`$ | $``$ | $`\varphi _3+1\varphi _2=\varphi _1^{}`$ |

would be $`T_2(\stackrel{}{\varphi })`$

$$\left(\begin{array}{c}\varphi _1^{}\\ \varphi _2^{}\\ \varphi _3^{}\end{array}\right)=\left(\begin{array}{c}1\\ 1\\ 1\end{array}\right)+\underset{\text{𝕄}_2}{\underset{}{\left(\begin{array}{ccc}0& 1& 1\\ 0& 1& 0\\ \mu & 1& 0\end{array}\right)}}\left(\begin{array}{c}\varphi _1\\ \varphi _2\\ \varphi _3\end{array}\right)$$
(3)   
and so on. Once we have defined all possible firing maps for a given number of oscillators we can proceed to deal with the attractors or fixed points of the system dynamics. As has been proved in  these fixed points must be cycles of $`N+1`$ firings. We define a cycle as a sequence of consecutive firings where each oscillator fires once and only once. Mathematically, each cycle is described by means of a return map. The return map is the transformation that gives the evolution of $`\stackrel{}{\varphi }`$ during a cycle and is the composition of all firing maps involved in the firing sequence of that cycle

$$\stackrel{}{\varphi }^{}=T_{c_1}T_{c_2}\mathrm{}T_{c_{N+1}}(\stackrel{}{\varphi })\stackrel{}{R}_c+\text{𝕄}_c\stackrel{}{\varphi },$$
(4)   
where $`T_{c_i}T_{c_j}(\varphi )`$ is the usual composition operation $`T_{c_i}(T_{c_j}(\varphi ))`$ and

$`\stackrel{}{R}_c=\stackrel{}{1}+{\displaystyle \underset{i=c_1}{\overset{c_N}{}}}({\displaystyle \underset{j=c_1}{\overset{i}{}}}\text{𝕄}_j)\stackrel{}{1}\text{and}\text{𝕄}_c={\displaystyle \underset{j=c_1}{\overset{c_{N+1}}{}}}\text{𝕄}_j.`$

Note that not all possible combinations of firing maps are allowed, just those ones whose indices $`c_i`$ sum $`p(N+1)`$ without any partial sum equal to $`q(N+1)`$, where $`p>q`$ are positive integers.

As all firing maps are linear transformations, return maps are also linear. There are $`N!`$ possible cycles in the $`N+1`$ oscillators case (all permutations of firing sequences with the initial firing oscillator $`\varphi _0`$ fixed). Following our previous example, for the four oscillators case all possible firing sequences and their associated return maps are

$`A:0,1,2,3T_1T_1T_1T_1`$   
$`B:0,1,3,2T_2T_3T_2T_1`$   
$`C:0,2,1,3T_1T_2T_3T_2`$   
$`D:0,2,3,1T_3T_2T_1T_2`$   
$`E:0,3,1,2T_2T_1T_2T_3`$   
$`F:0,3,2,1T_3T_3T_3T_3`$   
Now, in order to find the attractors of the dynamics, we must solve the fixed point equation

$$\stackrel{}{\varphi }_c^{}=\stackrel{}{R}_c+\text{𝕄}_c\stackrel{}{\varphi }_c^{},$$
(5)   
for every cycle $`c`$. Formally,

$$\stackrel{}{\varphi }_c^{}=\stackrel{}{R}_c(\text{𝕀}\text{𝕄}_c)^1.$$
(6)   
As was shown in , there are $`N`$ different stable solutions to the whole set of fixed point equations. Their stability is assured by the fact that $`\epsilon <0`$, since it guarantees that all eigenvalues of $`\text{𝕄}_c`$ lie inside the unit circle for all cycles $`c`$. In our four oscillators example these solutions are

$`\stackrel{}{\varphi }_A^{}=`$ $`(1,{\displaystyle \frac{3}{4+3\epsilon }},{\displaystyle \frac{2}{4+3\epsilon }},{\displaystyle \frac{1}{4+3\epsilon }})`$   
$`\stackrel{}{\varphi }_B^{}=\stackrel{}{\varphi }_C^{}=\stackrel{}{\varphi }_D^{}=\stackrel{}{\varphi }_E^{}=`$ $`(1,{\displaystyle \frac{1}{2+\epsilon }},1,{\displaystyle \frac{1}{2+\epsilon }})`$   
$`\stackrel{}{\varphi }_F^{}=`$ $`(1,{\displaystyle \frac{1}{4+\epsilon }},{\displaystyle \frac{2+\epsilon }{4+\epsilon }},{\displaystyle \frac{3+\epsilon }{4+\epsilon }})`$   
Which are a kind of four-oscillators traveling wave, chessboard and inverse traveling wave structures.

From now on we will label such solutions with index $`m`$ ($`m=1\mathrm{}N`$) since their first component always satisfy

$$\varphi _1^{}=\frac{m}{N+1+m\epsilon }.$$
(7)   
Therefore, in the example, we relabel patterns $`\stackrel{}{\varphi }_A^{}`$ as $`m=3`$, $`\stackrel{}{\varphi }_B^{}`$, $`\stackrel{}{\varphi }_C^{}`$, $`\stackrel{}{\varphi }_D^{}`$, $`\stackrel{}{\varphi }_E^{}`$ as $`m=2`$ and $`\stackrel{}{\varphi }_F^{}`$ as $`m=1`$.

Since there are $`N!`$ possible cycles and $`N`$ solutions to Eq. (7) there will be some fixed points or patterns which will appear more than once, so, we shall use $`C(N+1,m)`$ to characterize these degeneracies. In the example, the values of the degeneracies are $`C(4,1)=C(4,3)=1`$ and $`C(4,2)=4`$. In general, patterns which are solutions of cycle consisting in the iterative application of the same firing map (like A and F in our example) have no periodicities whereas the ones solution of mixtures of differents firing maps (B,C,D and E) have some periodic structure that are also solution of Eq. (7) for a case with less oscillators. In Fig. 1 we can visualize the solutions for $`N+1=2,3`$ and $`4`$ oscillators and realize that solution $`m=2`$ for the four oscillators case is a periodic composition of solution $`m=1`$ for the two oscillators case.

## III Pattern properties

As we have seen, the stability of all patterns solution of Eq. (6) is guaranteed by the fact that $`\epsilon <0`$, but the existence of such solutions is not ensured. In fact, for small values of the coupling strength $`|\epsilon |`$ all patterns do exist, but, as we increase it, some patterns disappear. The reason is that the solution loses its physical meaning because $`\varphi _1^{}>1`$. Their first component is always the one that becomes larger than unity earlier and this happens, for each $`m`$ and according to Eq. (9), when

$$\epsilon <\epsilon _m^{}=1\frac{N+1}{m}.$$
(8)   
Our coupling strength range of interest ends at $`\epsilon =1`$, since at $`\epsilon 1`$ we always find the same pathological dynamics which does not have any physical or biological sense. Realistic couplings never reach such higher values. Therefore, as $`\epsilon `$ runs from $`0`$ to $`1`$, all patterns whose $`m`$ satisfy $`m>\frac{N+1}{2}`$, disappear.

There is another interesting pattern property which has to do with the calculation of the pattern degeneracy $`C(N+1,m)`$. In principle, to calculate such degeneration, we should solve fixed point Eq. (6) for all possible cycles and count how many of them lead to the same pattern. Although for few oscillators the problem is quite straightforward, as we deal with higher and higher number of oscillators, the number of cycles increases (it grows as $`N!`$) and solving Eq. (6) becomes more difficult. Fortunately, there is another way of calculating $`C(N+1,m)`$ which reduces the problem to a combinatorial question. Lets show it through an example, in the previous four oscillators case, if we count, for each firing sequence, the number of oscillators which have received the pulse before firing, we can easily realize that this number is the same as its value of $`m`$

$`A:`$ $`0,\overline{1},\overline{2},\overline{3}`$ $`m=3`$   
$`B:`$ $`0,\overline{1},3,\overline{2}`$ $`m=2`$   
$`C:`$ $`0,2,\overline{1},\overline{3}`$ $`m=2`$   
$`D:`$ $`0,2,\overline{3},\overline{1}`$ $`m=2`$   
$`E:`$ $`0,3,\overline{1},\overline{2}`$ $`m=2`$   
$`F:`$ $`0,3,2,\overline{1}`$ $`m=1`$   
Here an upper bar means that the oscillator has already received a pulse during the cycle. The point is that it turns out that every pattern $`m`$ corresponds to a sequences of firings involving exactly $`m`$ oscillators that, when they do fire, had already received a pulse from their leftmost neighbor. Therefore, this property (we have checked for several values of $`N+1`$) allows us to associate every cycle with the pattern it leads to, just by counting these kind of firings. Now, calculating $`C(N+1,m)`$ becomes a straightforward matter. In Table I we have computed $`C(N+1,m)`$ for several values of $`N+1`$.

Apart from brute force counting, degeneracy distribution $`C(N+1,m)`$ can also be determined from the following relation

$`C(N+1,m)`$ $`=`$ $`mC(N,m)+`$ (10)   
$`(N+1m)C(N,m1),`$   
for $`2mN1`$. This recursion relation is closed by

$$C(N+1,1)=C(N+1,N)=1,$$
(11)   
which correspond to the firing sequences

$`0,N,(N1)\mathrm{}2,\overline{1}\text{and}0,\overline{1},\overline{2}\mathrm{}\overline{(N1)},\overline{N},`$

respectively.

From the previous relations one can deduce by induction the symmetry of the distribution with respect to its extremes at $`m=1`$ and $`m=N`$

$$C(N+1,m)=C(N+1,N+1m),$$
(12)   
and

$$\underset{m}{}C(N+1,m)=N!.$$
(13)   
Another interesting property is the period $`\mathrm{\Delta }_m^{N+1}`$ of each spatio-temporal pattern $`m`$. Since all oscillators are in a phase-locked state, they must oscillate with the same period. Then, as the intrinsic period of each oscillator is one, and when any oscillator receives the delaying pulse from its neighbor it has a phase equal to $`\varphi _1^{}`$, one can easily realize that the effective period is

$$\mathrm{\Delta }_m^{N+1}=1+\epsilon \varphi _1^{}=\frac{N+1+2m\epsilon }{N+1+m\epsilon }.$$
(14)   
Therefore, the larger the value of $`m`$, the longer the period of its associated pattern. It is important to notice that we have not fixed the value of such periods (each pattern has its own which is different from the others), since there are some authors who fix all periods equal to some constant, and use it as a condition to find the structures.

## IV Pattern selection

Once we have characterized all spatio-temporal patterns, we proceed to find some general formula which give us some estimation of the probability of each pattern to be selected, or in other words, an estimation of the volume of its basin of attraction. In order to achieve this objective, we should understand the mechanism which lead to the selection of a certain spatio-temporal structure and how is it modified as the parameters of the model ($`\epsilon `$ in our case) change.

There is an easy and straightforward way to get the essential features of this mechanism assuming that the probability of one oscillator to fire next is, basically, proportional to its phase (that is, if it has a phase slightly below $`1`$ it has a higher probability to be the next firing oscillator, whereas if it has a smaller phase, it will rarely fire next). Imagine the phases of all oscillators randomly distributed over the interval $`(0,1)`$. Then we let the system evolve till one of the oscillators reaches a phase $`\varphi _i=1`$ and emits a pulse that is received by its rightmost neighbor which lows its phase by an amount $`\epsilon \varphi _{i+1}`$. Now we assume that all phases are again randomly distributed over $`(0,1)`$ except the one which received the pulse whose phase is distributed over $`(0,1+\epsilon )`$. So, we get rid of memory effects (we know the oscillator that has fired should, now, have a phase equal to zero) and just keep in mind if each oscillator has received a pulse or has not. Therefore, the point is that under this conditions, the probability that one oscillator which has still not received a pulse do fire is some constant and, on the other hand, for the ones which had, is this constant times the factor $`(1+\epsilon )`$. Then, we can characterize the probability of having some cycle just by recalling how many oscillators do fire having previously received a pulse during that cycle. Basically, this probability is proportional to $`(1+\epsilon )^n`$ where $`n`$ stands for the number of oscillators which do fire having already received a pulse (the product of all constant terms will be absorbed in a normalization factor). This approach, where we assume all firings as almost-independent events, can be viewed as a kind of mean-field approximation. Then,as has been shown before, since cycles leading to the same pattern $`m`$ always exactly have $`m`$ oscillators that do fire having received the interacting pulse, we can give an estimation of the probability for pattern $`m`$ selection in the $`N+1`$ oscillators case

$$p_m^{N+1}(\epsilon )𝒩(\epsilon )C(N+1,m)(1+\epsilon )^m.$$
(15)   
Here $`𝒩(\epsilon )`$ is chosen so that summation of the probabilities over m gives 1

$$\underset{m}{}p_m^{N+1}(\epsilon )=1.$$
(16)   
In the limit of small coupling strength $`\epsilon 0`$, which is the more interesting case for the majority of physical and biological systems, one can assume that interaction plays almost no role when pattern selection takes place. That is, the fact that one oscillator has received the pulse from its neighbor does not low its probability to fire as the pulse does not modify appreciably its phase. Then, we can consider that all cycles have approximately the same probability to be selected, $`(1+\epsilon )^m1`$, and only pattern degeneracy has to be considered to get a good estimation of $`p_m^{N+1}`$

$$p_m^{N+1}\frac{C(N+1,m)}{N!}.$$
(17)   
The dominant pattern, that is, the one which has the larger probability to be selected coincides with the mean value of $`m`$ (due to the symmetric behavior of $`C(N+1,m)`$).

$$<m>_{N+1}=\underset{m}{}m\frac{C(N+1,m)}{N!}=\frac{N+1}{2}.$$
(18)   
For an odd number of oscillators $`<m>_{N+1}`$ does not exist and we have a competition between the two closest patterns $`m=N/2`$ and $`m=(N+2)/2`$. Recall that the most probable patterns turn out to be the ones with ”shortest wavelengths”, a fact that was already reported in simulations of these sort of systems. In Figs. 2 and 3 we check this new approximation for the $`N+1=10`$ and $`9`$ case and realize that expected results are in good agreement with simulations data.

There also is the interesting question of how does this probability distribution modifies when the number of oscillators increases. In Fig. 4 we show $`p_m`$ for different values of $`N+1`$. Since there are more possible values of $`m`$ available, as we increase $`N+1`$, $`p_m^{N+1}`$ diminishes. The distribution also gets narrower as we increase $`N+1`$ and this becomes clear when one studies the variance of $`p_m`$. It can be found that

$`<m^2>_{N+1}`$ $`=`$ $`{\displaystyle \underset{m}{}}m^2{\displaystyle \frac{C(N+1,m)}{N!}}`$ (19)   
$`=`$ $`{\displaystyle \frac{(N+1)^2}{4}}+{\displaystyle \frac{N+1}{12}}.`$ (20)   
We could not prove this without an explicit expression for $`C(N+1,m)`$ but we have checked it N up to $`170`$. Therefore

$$\sigma _{N+1}^2=\frac{N+1}{12}=\frac{<m>_{N+1}}{6}.$$
(21)   
It turns out that for a large number of oscillators almost all initial conditions lead to a pattern whose $`m`$ approximately falls in the interval $`<m>_{N+1}\pm \sqrt{<m>_{N+1}}`$. In order to compare it for different number of oscillators we have to normalize $`m`$ dividing by $`N+1`$. In that case, one observes that $`\sigma _{N+1}^21/\sqrt{N+1}`$ so that as we increase $`N+1`$, the spread of $`p_m^{N+1}`$ diminishes getting the distribution sharpened.

As Eq. (14) does not take into account the disappearance of the different patterns $`m`$ at the different values of $`\epsilon _m^{}`$ predicted by Eq. (9), it can not give a good quantitative estimation of pattern selection for higher coupling values. Nevertheless we can expand Eq. (14) to the leading order in $`\epsilon `$. For small $`\epsilon `$, $`p_m^{N+1}`$ are approximated by

$$p_m^{N+1}\frac{C(N+1,m)}{N!}(1+(m\frac{N+1}{2})\epsilon ).$$
(22)   
In Fig. 5 we compare this approximation with simulated data. The slopes near $`\epsilon =0`$ do agree with Eq. (21).

In our simulations we calculate the probability of each pattern to be selected just by counting how many realizations (with $`\varphi _0=1`$ and the rest of oscillators with random initial conditions) lead to each pattern $`m`$ and divide over the total number of realizations.

Although we only have a good quantitative estimation of $`p_m^{N+1}`$ for small values of $`\epsilon `$, Eq. (15) catches the two basic mechanisms responsible of pattern selection. On the one hand, it is clear that for higher values of the coupling strength $`|\epsilon |`$, when one oscillator receives a pulse, it lows its phase to almost zero and, consequently, its firing probability also does. Therefore pattern selection probability $`p_m^{N+1}(\epsilon )`$ is strongly controlled by the number of oscillators which have to fire having already received a pulse, that is, the probabilistic factor $`(1+\epsilon )^m`$. As a consequence, $`p_m^{N+1}`$ begin to decrease sooner when $`|\epsilon |`$ increases, the larger $`m`$ is. On the other hand, for small values of the coupling strength, interaction plays almost no role and $`p_m^{N+1}(\epsilon )`$ is dominated by the degeneracy factor $`C(N+1,m)`$. Therefore $`p_m^{N+1}(\epsilon )`$ for the different values of $`m`$ are basically ordered as $`C(N+1,m)`$. In Fig. $`6`$, $`7`$ and $`8`$ we show results from simulations of $`p_m^{N+1}(\epsilon )`$ for different number of oscillators.

## V Conclusions

In this paper we have studied some properties of the spatio-temporal patterns that appear in a ring of pulse-coupled oscillators with inhibitory interactions. We have focused our attention in estimating the probability of selecting a certain pattern under arbitrary initial conditions and have shown the two basic mechanisms responsible of that: the degeneracy distribution $`C(N+1,m)`$, for small values of $`\epsilon `$, and $`m`$, the number of oscillators that do fire having already received a pulse, for higher values of $`\epsilon `$. According to this, the different probabilities of selecting pattern $`m`$ start being distributed following the degeneracy distribution $`C(N+1,m)`$, and, as $`\epsilon `$ decreases, these probabilities diminish in a hierarchical way: the larger the value of $`m`$, the sooner its selection probability is going to decrease, so that only patterns with smaller m will survive for higher values of $`\epsilon `$. Moreover, some of the structures disappear, at the different values of $`\epsilon _m^{}`$, during this process. We have found out an approximation formula for $`p_m^{N+1}(\epsilon )`$ which takes into account all these mechanisms and gives us a quantitative estimation of the different selection probabilities for small $`\epsilon `$.

The estimation of the volume of the basin of attraction of each spatio-temporal pattern $`m`$ also gives us an idea of the stability of the different structures with respect to additive noise fluctuations (for instance, we can add some random quantity $`\eta `$ to all phases after each firing event or a continuous-time $`\eta (t)`$ in the driving). Simulations of arrays of noisy pulse coupled oscillators showed that our most probable patterns were also the most stable. The present paper only concerns spatio-temporal pattern formation in a ring of oscillators, nevertheless, all results are trivially generalized to bidirectional couplings. Although the question of what happens when dealing with higher dimension lattices remains opened, some simulations results in 2d  showed that almost all realizations lead to a chessboard pattern in analogy with our results in the ring. That makes us believe we have caught the basic features of the problem in our 1d model.

## Acknowledgments

The authors are indebted to C.J.Pérez and A.Arenas for very fruitful discussions. They also acknowledge extremely constructive suggestions from an anonymous referee. This work has been supported by DGICYT of the Spanish Government through grant PB96-0168 and EU TMR Grant ERBFMRXCT980183. X.G. also acknowledges financial support from the Generalitat de Catalunya.