Synchronization of moving integrate and fire oscillators

We present a model of integrate and fire oscillators that move on a plane. The phase of the oscillators evolves linearly in time and when it reaches a threshold value they fire choosing their neighbors according to a certain interaction range. Depending on the velocity of the ballistic motion and the average number of neighbors each oscillator fires to, we identify different regimes shown in a phase diagram. We characterize these regimes by means of novel parameters as the accumulated number of contacted neighbors.


Introduction
Among the many emerging phenomena we observe in nature, synchronization is one of the most paradigmatic examples. It can be roughly understood as the collective dynamics of units whose internal state evolves periodically in time and when they interact tend to synchronize their internal variables [Pikovsky et al., 2001]. The achievement of the final synchronized state (if any) strongly depends on the interaction pattern of the system [Boccaletti et al., 2006;Arenas et al., 2008]. Up to now, synchronization has been mainly analyzed in fixed topologies but we are witnessing the first evidences that links between agents can evolve in time [Olfati-Saber et al., 2007;Onnela et al., 2007;Valencia et al., 2008]. The case in which this evolution of the network topology is an effect of the agents mobility is a particularly interesting case [Buscarino et al., 2006;Tanner et al., 2003;Buhl et al., 2006]. The effect of this changing patterns of interaction on synchronization features has been analyzed in different settings, for instance in chemotaxis [Tanaka, 2007], mobile ad hoc networks [Römer, 2001], wireless sensor networks [Sivrikaya & Yener, 2004], and the expression of segmentation clock genes [Uriu et al., 2010].
In the recent literature, studies on synchronization in dynamically evolving complex networks have been mainly concentrated in the case when the topology changes very fast. This is the so-called fastswitching approximation (FSA) [Belykh et al., 2004;Frasca et al., 2008;Porfiri et al., 2006;Stilwell et al., 2006], which replaces the real interaction between agents by the "mean field assumption" that all agents interact with an effective strength that corresponds to the probability that any pair of agents are connected.
Recently, it has been proposed a general framework of mobile oscillator networks where agents perform random walks in a two-dimensional (2D) plane [Fujiwara et al., 2011b]. It has been shown that FSA fails when the time scale of local synchronization is shorter than the time scale of the topology change due to the agent motion. New behaviors arise due to the interplay between instantaneous network topology, agent motion, and interaction rules. This framework, that reduces to FSA when velocity is high enough, is valid for models whose evolution can be well approximated by linear dynamics. This actually holds for models such as populations of Kuramoto oscillators [Kuramoto, 1984;Acebrón et al., 2005], whose evolution, after a short transient time, is very well described by a set of linear equations that can be solved in terms of spectral properties of the Laplacian matrix [Fujiwara et al., 2011a].
In the present paper, we focus on a dynamical system, a population of integrate and fire oscillators (IFO), where linearization is not a good approximation, since the evolution takes place in two different time scales. One for the slow evolution of the internal state variables (the phase and the orientation) and the other for the fast interaction between the units (pulse coupling). During the last years it has been shown that the interaction structure plays a fundamental role in the dynamics of IFO networks. Zillmer et al. [2009] observed different dynamical regimes due to network connectivity in a system formed by inhibitory integrate-and-fire neurons that were randomly connected. Also, the underlying network structure can affect the speed with which the system reaches the synchronized state, as studied by Grabow et al. [2011]. Usually, IFO have been used to model neural systems but we can also find some examples of applications in other fields, as for example in economy [Erola et al., 2011]. Models where the oscillators do not remain fixed, and the network of interactions changes with time can find a direct application in biological systems such as flashing fireflies, that interchange light signals while searching for potential mates [Mirollo & Strogatz, 1990;Ramirez-Avila et al., 2011].
In the present case we will show that the interplay between agents motion and phase evolution towards a synchronized state presents different asymptotic behaviors, reminiscent of the observation in Kuramoto oscillators [Fujiwara et al., 2011b] and agents using communication protocols [Baronchelli & Diaz-Guilera, 2011]. We identify, furthermore, the possible mechanisms in the different regions of the parameter space.
The organization of the paper is as follows. In the next section we introduce the model. Then we show the results for different regions of the parameter space, velocity of the agents and range of interaction, and later we identify the different microscopic mechanisms that lead the system to a globally synchronized state.

The model
We propose a setting in which a population of N integrate and fire oscillators (IFOS) [Mirollo & Strogatz, 1990] move at a constant velocity V in a bidimensional plane of size L with periodic boundary conditions. Each agent has two degrees of freedom corresponding to an internal phase φ ∈ [0, 1] and orientation θ i ∈ [0, 2π], both randomly set in an uniform manner at the initial configuration.
The evolution of the system takes place on two different timescales. The slow timescale sets the pace at which the phases of the agents increase uniformly with period τ , until they reach a maximum value of 1, when a firing event occurs. Then the phase is reset and the oscillator is randomly reoriented. Upon this event at time t, the firing oscillator influences its nearest neighbors (oscillators at minimal distance) producing an update in their phases by a factor ǫ: . (2) The phase and orientation resetting corresponds then to the fast time scale. If the neighbor's phase update overcomes the phase maximum, another firing event is triggered and this process goes on repeatedly until all shots have ceased. At this point, the time t runs again until the next firing event occurs. The system is synchronized when we encounter in the system a succession of consecutive firing events (avalanche) equal to the system size N , since after this fact all oscillators will remain synchronized forever because all of them will have the same period τ with or without interactions. For the sake of clarity we define the (discrete) time T , as the number of times a given oscillator (that we will identify with oscillator 1 in our computer simulations) has fired. This allows us to define T sync as the minimum number of (integer) cycles this reference oscillator takes to enter the synchronized state (i.e. the number of updates needed for an avalanche of size N to occur). We propose a geometric condition for neighbor selection upon a firing event as shown in figure 1: Every agent scans a circular area of radius R around it and shots the neighbors therein. We introduce a parameter r ∈ [0, 1] that indicates the fraction of the system available for interaction and relates both R, L variables and the average outgoing degree of the nodes of our evolving network Throughout this paper, we have used fixed parameters L = 100, τ = 1, N = 50 and ǫ = 0.02 ∼ O(1/N ), while analyzing the explicit dependence on the mobility parameters, r and V . Before proceeding to show the results of our simulations, we need to note that the type of proposed interaction range in this system has been reported to show statistical properties similar to a continuous percolation [Fujiwara et al., 2011b], that in the case of static oscillators occurs for approximately r c ≈ 4.51/(N − 1) = 0.09 [Dall & Christensen, 2002;Balister et al., 2005]. In our range of study, we hope to observe some traces of this percolation as well as saturation properties observed in other moving oscillator systems at high speeds [Fujiwara et al., 2011b].
It is also important to notice that we have kept the dynamical evolution of the units at its maximal simplicity since we are mostly interested in the interplay between motion (and hence construction of a dynamical network) and internal dynamics and how synchronization emerges as a collective property of the system.

Results
We present in figure 2 the results of our simulations. A preliminary observation points out that the roles and importances of V and r change throughout our map.
For high enough values of r (r r c ) , the synchronization time, T sync , is almost unaffected by the values of the speed V . Although this time is dependent on r, its range of possible values is much narrower (by orders of magnitude) than below the value, r c , that characterizes the static percolation transition. Actually, below this critical value the velocity plays the crucial role since for a fixed value of r the synchronization time changes by a factor of 10 4 . To bring further insight to the map, we show a profile of the efficiency of the system rT sync dependence on r for various velocities. This efficiency variable balances the range of interactions (which has an "energy cost" related with the number of shots to different neighbors and their strength) and the number of shots that the system needs to synchronize (related with T sync ).
It seems clear that mobility of agents helps to minimize energy consumption but as we approach the critical percolating value of r c , we observe that a range of behaviors emerge. On one hand, we observe for high velocities that the efficiency of the process remains roughly constant independently of r, since the extreme mobility of the agents compensate the reduced range of its interactions and successfully diffuses the synchronization process around the system, a process that is equivalent to the observation in other settings [Frasca et al., 2008;Fujiwara et al., 2011b;Baronchelli & Diaz-Guilera, 2011]. On the other hand, if the mobility of the agents is reduced, then the path to synchronization is more lengthly as well as more energy consuming. Synchronization is still possible (below the percolation static limit r c ) but the time to achieve it grows very fast, even resulting in an effectively infinite time 1 .
These results lead us to identify different regimes and the consequent transitions between them. On the one hand, we find a "topological" transition at the critical static value r c , since above it there is basically no velocity dependence whereas below r c the influence of V is determinant. On the other hand, below r c , where synchronization is made possible by the mobility of the agents, we can identify a transition (depending on both r and V ) that separates two dynamical regimes.

Mechanisms
Up to this point, we have studied the efficiency of our system and detected several regimes strongly dependent on two factors, the mobility of our agents and its range of interaction.
In this part of the paper we want to study the distinct mechanisms that the system uses in its path to synchronization. To this end, we introduce an order parameter η(T ) = cos (2πφ(T )) that is an increasing function that measures the overall synchrony of our system, ranging from a uniform phase distribution of our oscillators (η(T ) = 0) to complete synchronization 2 (η(T ) = 1).
To couple our V and r control parameters we also introduce the number of distinct interactions per oscillator N i c (accumulated encounters with different agents by a oscillator i). This value is bounded (on average) between a minimal starting value of N c 0 = (r−1)N and a maximal value of N c max = N −1 and it provides information about the evolution of our system's synchronizing mechanism. Since the bounding of this value depends on r, we introduce a normalized magnitude χ, that is a quantification of the mixing of our system. When the mixing is minimal ( N c = N c 0 ) χ is 0.
As the system mixes, i.e. the oscillators increase its average number of contacted neighbors, it can grow up to its maximum value χ = 1, i.e. N c = N c max . We have calculated the time evolution of both η and χ for different values of V ∈ {0.1, 1, 100} and r ∈ {0.05, 0.1, 0.2}. In figure 3 one can see the evolution of both parameters measured at the same time instants (and sufficiently averaged over enough realizations) together with the difference between them over time η(T ) − χ(T ) that give insight about the topological evolution of our system as it synchronizes. The figures show three clearly distinct patterns.
For high velocities we observe in fig. 3a) a gradual increase of the order parameter and a minor influence on r at fixed V , indicating that the synchronization emerges evenly on the system in a global fashion, due to the quickly changing topology of the network (neighbors of a given oscillator change rapidly). This regime (which we call diffusive) requires the inter-contact of the majority of the system, but this circumstance is rapidly achieved due to the strong mobility of the agents. In fact, this regime is optimal as far as the synchronization time is concerned, since the interactions are more effective. These conclusions were obtained for populations of Kuramoto oscillators [Fujiwara et al., 2011b], for which this regime corresponds to the region of validity of the FSA.
In the opposing case, when velocities are small enough, this behavior is completely lost and a step function appears, indicating that the slow mobility of the agents allows them to synchronize locally with their neighbors creating islands of synchrony. The sudden increase of the order parameter occurs at a regular pace, fact that points out that whenever the islands are disbanded (change of neighbors), they still transmit the local synchrony to the neighboring groups, mechanism that allows for system synchrony while keeping χ in small values. The initial height of the steps is dependent on r and decreasing as χ(T ) grows due to the limited range of η(T ) available states. Finally, as we decrease r approaching the critical value r c , we observe a transition from a local to a bounded regime, where the synchronizing time is so long that again allows for the interaction of the majority of the agents among themselves upon synchronization time (due to the bounded nature of the system). In this regime, the range of interaction is very reduced, and so is the size of the clusters, so an agreement between the multiple clusters created (if any) comes after almost all the system has interacted. Consequently, the increasing of η with χ is slower (many small steps, see fig. 3a) ) while the final value χ becomes larger ( fig. 3).
In fig. 3b) we provide an explicit time evolution of the difference η(T ) − χ(T ) in order to make the three regimes and the influence of r better identified.
It is interesting to study the final mixing of our system upon synchronization as shown in figure 4. This value χ(T = T sync ) ≡ χ sync together with T sync characterize the evolution of the system towards the synchronized final state. These features depend on both r and V . At fixed velocity the final mixing of the system decreases as the interaction range grows. This is caused by the fact that although an increased r induces more mixing (as oscillators find new neighbors more easily) it also drastically reduces (below the critical values r c ) the synchronizing time T sync thus reducing the chances of encounters between different oscillators. Above r c we observe a saturation of the values as the dependence of T sync in r and V is practically lost. Figure 4 b) shows a change between the relation pattern of χ sync and efficiency rT sync ranging from the diffusive regime (mixing independent of rT sync ) to the curves where both the local and bounded regimes are shown. It is important to note the similarity of the observed shapes (for a wide range of V values) where we find the local phase concentrated around the minimum value of χ sync that gradually grows in a power law fashion as the performance of the system decreases (it consumes more energy to synchronize).
The introduction of the new parameter χ sync allows us to present a phase diagram of our system relating the overall performance (in terms of efficiency) with the synchronizing mechanism used (figure 5). We identify the diffusive regime in the zone of high velocities V ∼ O(10) where both values of χ and rT sync are almost independent of r. This zone falters into the bounded one as V and specially r decrease, where both the efficiency and the mixing of the system is reduced. Finally for small enough velocities the local zone is clearly visible with low values of system mixing. From the map one clearly observes that the most   beneficial synchronizing mechanisms (in terms of energy consumption) are the diffusive and local ones.

Conclusions
Following the recent literature on complex systems, one of the hottest topics is the relation between dynamics and topology of interactions. In particular, there many evidences that the patterns of interaction change rapidly with time, thereby completely altering and conditioning the dynamical properties of the system.
Here we have proposed a framework in which agents move on a plane and are allowed to interact in a pulsing way. Each agent, representing a phase oscillator, moves at a common velocity and changes its internal phase with a common period. When this internal phase reaches a threshold value, the oscillator "fires", thus resetting its own phase as well as its orientation. At the same time, this oscillator interacts with the neighbors within a certain distance by changing their phases. Notably, this interaction setting makes the system reach a final synchronized state in which all oscillators fire at unison within the same fast scale.
Keeping all geometrical parameters constant, we analyze the system behavior by changing solely the velocity of the agents and the interaction range. Notice that the motion of the agents is what ensures that the system will be able to synchronize. If the agents are static, a minimal fixed topology will be required to connect them (what is called a giant component in network terminology). In contrast, a population of moving agents (even when the interaction range is small) will eventually synchronize.
We measure the time needed for the system to synchronize as a function of the two relevant agent parameters, the velocity and the fraction of population they interact with. Our model can be applied, for instance, to the field of wire-less communications by introducing the performance, which stands for the total number of signals emitted by the population to reach the synchronized state. Our simulations show that the time required for synchronization ranges widely, depending on the speed of the agents. As expected, the optimal performance is achieved when agents move very fast, irrespective of the interaction range; however, synchronization in systems with low mobilities dramatically depend on the interaction range.
Finally, we have focused on the mechanisms that the system follows along its path to synchronization. We have introduced a new order parameter, that stands for the fraction of different units each oscillator has interacted with. This order parameter complements performance, and a novel phase diagram enables us to identify three different characterizable mechanisms leading to synchronization: i) the diffusive mechanism, where the system very quickly reaches synchronization by means of extremely effective interactions of each oscillator with a large fraction of the population; ii) the local mechanism, where agents only interact with a small population fraction, but local synchronization is sufficient to lead the system to the globally synchronized state; and finally iii) the bounded mechanism, characterized by very slow motion and short range of interaction, which allow a high degree of accumulated mixing during a very long synchronizing time.
The identification of these mechanisms that relate mobility and interaction to synchronization will undoubtedly be crucial in similar models of populations of moving agents. Indeed, understanding these phenomena will help to design optimal protocols to dissect more realistic settings, as well as to predict their behavior, by defining their interaction rules.