Synchronization to extreme events in moving agents

Interactions amongst agents frequently exist only at particular moments in time, depending on their closeness in space and movement parameters. Here we propose a minimal model of moving agents where the network of contacts changes over time due to their motion. In particular, agents interact based on their proximity in a two-dimensional space, but only if they belong to the same fixed interaction zones. Our research reveals the emergence of global synchronization if all the interaction zones are attractive. However, if some of the interaction zones are repulsive, they deflect synchrony and lead to short-lasting but recurrent deviations that constitute extreme events in the network. We use two paradigmatic oscillators for the description of the agent dynamics to demonstrate our findings numerically, and we also provide an analytical formulation to describe the emergence of complete synchrony and the thresholds that distinguish extreme events from other intermittent states based on the peak-over-threshold approach.


Introduction
Research to understand the interplay between complex networks and the dynamical properties of coupled oscillators has been a hotspot for the last few decades and the developing phenomenon of synchronization [1][2][3][4] is one of the most important dynamical processes that has been in the center of these researches. Cooperation [5,6] and time series analysis [7,8] in complex network have been studied in the past few years. From the perspective of synchronization among coupled oscillators placed into a complex network [9,10], the correlation between the network's topology and local dynamics is quite decisive. Here, synchronization signifies a process of adaptation to a common collective behavior of oscillators due to their interaction. In most of the previous studies of such systems, the network topology is assumed to be invariant over time and thus the system is controlled by a deterministic static formation for all the course of time. But such a crude assumption regarding the network connectivity inhibits one to model and study most of the practical instances.
Recently, time-varying networks have grabbed the attention of the researchers due to their enormous applications in various fields like functional brain network [11], epidemic modeling [12], communication systems [13,14] and many more. Time-varying networks, also known as temporal networks [15] indicate those networks in which links get activated for a certain course of time. On the assumption of time-invariant nodes which are static over time, many network architecture is studied, e.g. power transmission elements are considered as such nodes among which haphazard links are treated as the coupling between elements of the power transmission system [16]. Even for functional brain networks [11], these types of nodes are considered to characterize the dynamical evolution. Particularly, the scenario of time-varying networks owing to the mobility in the nodes is really a significant platform to study several dynamical processes over them in which nodes can move in the space and interact with each other based on their physical proximity. Collective behavior in random geometric graphs and random walkers has been studied [17,18]. For instance, in case of coordinated motion of

Mathematical model for the moving oscillators
We start with describing the mathematical framework of the network model of N moving agents which are distributed in a two dimensional (2D) XY-plane g g g g , , sq. units that contain m number of restricted non-overlapping circular zones, each of radius r units. The interactions among the agents are activated only when they enter into the same pre-defined circular zone. In a similar fashion, other types of interaction zones (closed) of square, rectangle or even ellipse shape could have been considered. But, we are considering here simple circular zones, without any loss of generalization (later we show that the analytical condition for synchronization only depends on the area of interaction coupling zone). For illustration, let us consider the upper left interaction coupling zone in figure 1, which contains only one agent (smaller yellow circle). So, this agent cannot interact with any other agents for that particular instant of time. Similarly, the remaining 9 moving agents (outside any coupling zone) in the space å are deprived of mutual interaction with that agent and among each other. But for the coupling zone containing 2 agents in the upper right of the figure 1, those two agents are only interacting between themselves for that occasion. The same logic is applicable for the interaction among the 3 moving agents in the lower left coupling zone in the figure. At some time, it may also happen that there will be a zone in which no moving agent is present, e.g. the empty zone in the bottom right. Now we discuss the scheme of random movement of N agents in the two-dimensional plane. Initially, N identical agents are randomly placed in the region å . Any kind of bumping collision among those random moving agents is forbidden. The agents are allowed to move in any direction t i q ( ) with velocity where v is the modulus of the agent velocity and t i N , 1,2 , , For better visualization, an exemplary trajectory of a moving agent is depicted graphically in figure 2. To ensure the fact that the moving agents always remain confined within the region å , the periodic boundary condition is applied. We can relate these moving agents with vehicles equipped with global positioning system device, which enables to track location. The movement of each agent is governed by the dynamical upgrading rule t where t y i ( ) is the position of the ith agent in the plane at any time t and ΔT is the motion integrating step-size. At each integration step, t i q ( ) and therefore Furthermore, a d-dimensional dynamical system is allied with each moving agent. As per the above stated assumptions, the dynamical equation of each agent can be written as , ..., , : given by the system's dynamics, K is the coupling constant, H(x j ) is the vector coupling function for diffusive type of interaction and the time-varying matrix G t g t is the Laplacian matrix describing the location of all agents and hence the connectivity pattern of the network at any time t with respect to the pre-specified coupling zones. The Laplacian matrix G(t) is defined as a symmetric, positive semi-definite and M-matrix and deliberately calculated as the difference between the degree matrix and the adjacency matrix of the network. Particularly, g ij (t)=g ji (t)=−1 if both ith and jth agents enter the same coupling zone at time t, and otherwise zero. Since any kind of self-interaction is intolerable, so we assign g ii (t)=h, where h is the number of agents in the same coupling zone at time t, to expect synchronization as one of the possible states of the whole dynamical system. The system given in equation (1) is Figure 1. Schematic diagram of a two dimensional space where N=15 agents (smaller yellow circles) are moving independently without any influence of any other agents' motion, in the periodic planar space g g g g , , ]with g=10. Four disjoint coupling zones (larger circles) are pre-defined in å . Whenever, the moving individuals enter the same zone, they interact with each other. Note that the agents from one coupling zone do not possess any kind of interaction with the agents belonging to any other coupling zones, they only bother about the members within the same coupling zone, not about the other members of å . integrated using the Runge-Kutta-Fehlberg method with a fixed time step defined as 7 ΔT=0.01. The initial conditions are chosen randomly from the phase space where the chaotic attractor resides.

Analytical condition for synchronization
In this section, we analytically derive the necessary condition for stable synchrony among moving agents in the two-dimensional plane with two attractive interaction zones. As mentioned earlier, synchronization refers to a process in which coupled systems regulate their dynamics to achieve a coherent evolution, and time-averaged Laplacian matrix G g ij = [ ]plays a crucial role to understand the synchronization behavior of a time-dependent networked dynamical system. According to the [56], under fast switching among all the possible network configurations, if the system of coupled oscillators given by with time-static coupling topology G g ij = [¯] goes through stable synchrony and there is a constant λ such that ( )¯, then there exists ò 0 such that for all fixed 0<ò<ò 0 , the system described as (with the time-varying Laplacian G(t/ò)) also supports stable synchronization. So if the time-average of the ( ) exhibits synchrony of the system, then the time varying network will also possess synchronization. The synchronization of the time-evolving network can thus be predicted and assured if the time-averaged network G supports synchronization. Henceforth, we now focus to calculate the time-averaged coupling matrix G .
To simplify the calculations, first we consider two disjoint confined zones with the same area in the physical space å . Let, A={1, 2, K, N} be the set of all moving agents in the physical space å and B A Í be the set of all interacting agents lying within any one of the pre-defined coupling zone and C A Í be another set of all interacting agents lying in any other disjoint pre-defined coupling zone, such that there is no common agents lying in both sets B and C. In other words, B C Ç = AE. Let G B denotes the corresponding Laplacian matrix for the set B and B i be the resulting sum over all the coupling matrices generated by the ith agent. Let G A be the zerorow sum Laplacian matrix with ) agents except the agents j and k, we have to choose i 2 -( )elements, which can be done in N i . Then the ( j, k)th element of the matrix B i is where k ä B and with the same assumption that B contains only i agents, we have to deal with (i−1) elements excluding the ith agent, which to be selected out of remaining (N−1) elements. Hence, the diagonal elements of B i are Thus, from equations (4) and (5), we can write Let, out of N moving agents, i individuals go to one restricted zone, while j moving individuals go to another cramp zone. Note that those i agents are completely different from the other j agents and their motion is completely independent from others. Let, p i be the probability that i agents are coupled with each other. Then, where P is the probability of an agent for lying within a coupling zone. In our case, P where r is the radius of the circular zones. Then the time-averaged matrix becomes G pB pB for m disjoint restricted zones with same area. Since, G A is a zero-sum matrix so that the manifold of synchronous states is neutrally stable, G also possesses the same property. Hence, the eigen values of G are λ 1 =0 and mP N, , the coupling term provides sufficient amount of coherence among the oscillators to reach synchronization for the whole system.

Numerical findings
In order to look at the dynamical behavior of the network numerically, we first define the synchronization error ] is taken into consideration as the coupling function. In the following subsections, our attention will be to identify the synchronization region by changing the network parameters, namely the coupling strength K, modulus of velocity v and the number of attractive interaction zones m. We fix the parameters associated with the local dynamics of each individual agent in the chaotic regime.

Lorenz system
We first consider the chaotic Lorenz system [42], where the state dynamics of each moving individuals is represented by Without loss of any generality, we now meditate on the appearance of synchronization for N=100 identical Lorenz oscillators moving with modulus of velocity v=5.0 and m number of interaction zones of radius r=4.0. We perceive that the synchronization error E descends to zero at K=1.0 and K=0.5 respectively for m=2 and m=4 number of attractive coupling zones, which convey well agreement with the critical interaction strengths K c =0.5 and K c =1.0 derived from the relation (8) reflecting the effect of m, for α=3. 15. We detect that with the increasing number of interaction zones, the agents can communicate in more space now, as a result of that E decreased to zero more rapidly for larger m.
Here we note that the obtained results remain valid for any other value of the network size N. Particularly, higher N would need lower values of the interaction strength K to achieve synchronization for fixed values of other network parameters (results not shown here). Even that inverse proportional relation between N and K is obvious from equation (8). Since, interaction zones are fixed in plane, mobility of nodes have huge significance in obtaining complete synchronization. On the other hand, the modulus of velocity v does not appear explicitly in the relation (8). So, in order to understand the simultaneous influence of K and v on synchrony, we plot the phase diagram in the K-v parameter space, in figure 3. We affix m=2 interaction zones of radius r=4 in the plane and observe variation in between the parameters K and v for randomly moving oscillators initially in the phase space. We notice that there exists an optimal interval of v for which complete synchronization can be assured depending on the values of the interaction strength K.

Rössler system
Next we take chaotic Rössler oscillator [43] into consideration which is described as With the same number of agents as before, and v=2.5 with m=2 and r=4.0, the synchronization error E with respect to K, drops from a non-zero value to zero at K c =0.039 and continues to be zero.  figure 4(a), we consider two coupling zones of radius r=4.0 as before, centered at g g 2, 2 ( )and g g 2, 2 --( ) and plot the K-v parameter space in terms of the emergence of synchronization. As can be seen, for the whole range of K ä [0, 0.1] and v ä [0, 2.5], the higher the K, the lower the v is needed to achieve synchrony. Next, we add two more attracting coupling zones centered at g g 2, 2 -( ) and g g 2, 2 -( ) with all the other parameters fixed as before and plot the K-v parameter space in figure 4(b). A similar sort of behavior as in figure 4(a), is observed here but with a significant enhancement in the synchronization region, as the agents are now able to interact inside two more attracting coupling zones.

Inclusion of repulsive coupling zone and emergence of extreme events
Let us now move on to examine the impact of the presence of repulsive zones in the movement space. For this, we will investigate both Lorenz and Rössler oscillators as before.

Lorenz system
To scrutinize the scenario, two interaction coupling zones centered at g g 2, 2 ( )and g g 2, 2 -( ) , each of radius r=4.0, are considered on which N=100 identical Lorenz oscillators are randomly moving in any direction with uniform modulus velocity v=5.0. We place the coupling zone centered at g g 2, 2 ( )as the attractive one and treat the other one as repulsive coupling zone, with K a and K r as the attractive and repulsive coupling strengths respectively. Under this arrangement, we spot that the error trajectory continues to be on the invariant synchronous manifold (i.e. E=0) for most of the time, but due to the presence of the repulsive coupling zone with adequate strength, error trajectory experiences non-uniform, uncertain extensive expeditions from the synchronous manifold (see figure 5(a)). This is because the attractor moves near a saddle orbit from where the trajectories get repelled resulting in large excursions. After such deviation, those trajectories come back to the invariant manifold due to stretching and folding nature of chaotic attractors. Such ). For (b), we place m=4 coupling zones centered at locally repulsive deflection designates on-off intermittent [53,59] like behavior of the invariant manifold containing nonlinear chaotic orbits. Whenever the error at a certain instant of time crosses the threshold H S (details discussed later), that event can be considered as an appropriate candidate for extreme events.
To confirm whether the perceived spikes in figure 5(a) reveal the existence of extreme events or not, a sufficiently long time series of the variable E is accumulated and we have constructed a probability density function (PDF) for the event sizes E n defined as the local maximum values of E. The collected time series is sufficient enough so that due to statistical regularity, inclusion of new sample does not affect the structure of the observed L-shaped PDF [52]. From figure 5(b), it is quite evident that the occurrence of such event sizes is much more than that according to Gaussian distribution and the corresponding histogram possesses a long tail, which guarantees the existence of extreme events.
The presence of events with diverse sizes necessitates to prescribe a criterion that discriminate extreme events from other small-sized events. So, for further characterization of the extremeness, using peak over threshold (POT) approach, we make an attempt based on the specification of a threshold value H S (say) indicating a significant height that should be crossed by those rare events in order to be counted as an extreme event.
According to the Fisher-Tippet theorem [60] in extreme value theory [61], three types of distributions are mainly used to model the maximum or minimum of the collection of random observations from the same distribution. Specifically, they are the Gumbel, Fréchet, and Weibull distribution [62]. In general, Gumbel and Fréchet distributions correlate with largest extreme value, whereas Weibull model deals with the smallest extreme value [63][64][65]. Moreover, the Gumbel distribution is unbounded defined on the entire real axis and the Fréchet distribution is bounded for the lower side, x>0 but has a heavy upper tail. So, in seek of the threshold of extreme value indicator, we consider the PDF of a Weibull random variable as follows where p>0 is the shape parameter and λ>0 is the scale parameter of the distribution [65]. Next we proceed by adopting the method developed by Massel [66]. In the asymptotic limit of infinite time observational window, the distribution function can be expressed in terms of Weibull random distribution with where σ is the the standard deviation of peaks from a time series and under this assumption, the error height H is twice the envelope, that is H x x 2 for 0  = . Thus,   [68][69][70].
Now, for our network model, we calculate the mean and standard deviation of the error and we have the threshold which is equal to mean plus eight times of the standard deviation calculated over the time interval as in figure 5(a). We set this threshold as a horizontal line (red line) in the figure 5(a) plotted over the variation of error with respect to time.
In figure 6, we identify different dynamical states, namely synchronization, intermittent spikes, extreme events as a consequence of on-off intermittent states and desynchronization states in the K a -K r coupling parameter space within the ranges K a ä [5,10] and K r ä [−0.55, 0.0]. We particularly demonstrate the sensitivity of the negative coupling parameter K r on the network dynamics. For smaller values of K r , the attractive interactions dominate and the networked system maintains complete synchronization for suitable optimal values of K a . However, as soon as K r becomes higher, we need much higher K a to retain synchrony, otherwise, with increasing K r depending on K a , the synchronization becomes intermittent and gradually generate extreme events. This scenario is depicted through the EE region in the phase diagram. If we continuously increase K r , the synchrony gets completely deflected producing desynchronization.

Rössler oscillators
Again, we inspect the influence of the presence of a repulsive zone in the movement space by keeping fixed Rössler oscillators in each of these moving agents. We fix one attracting coupling zone centered at g g 2, 2 ( ) and one repulsive coupling zone centered at g g 2, 2 -( ) , each of same radius r=4.0, having the same area. Under this framework, N=100 Rössler oscillators, with modulus velocity v=2.5, display high amplitude fluctuation from the synchronous manifold (see figure 7(a)) for attractive coupling strength K a =1.0 and repulsive coupling strength K r =−0.1. As discussed earlier, the chaotic dynamics of the network was confined to the original synchronization manifold as long as the coupling zones provide attractive interactions. But whenever along with attractive zone the movement space consists of a zone bringing appropriate repulsive force, the synchronization becomes intermittent and the dynamical network occasionally blows out of the manifold.
On-off intermittency in the form of aperiodic switching among the synchronous state (i.e. when the error is zero) and chaotic eruption of the oscillation is observed in the error dynamics (see figure 7(a)). Whenever the error at a certain instant of time crosses the threshold H S , that event can be considered as an appropriate candidate for extreme events. The spikes observed here are recurrent, aperiodic and of different amplitudes indicating the origination of extreme events. For a better understanding of this scenario, in figure 7(b), the histogram associated to the PDF of the event sizes is plotted from a sufficiently long time series of the dynamical units where E n stands for local maximum values of E (see equation (9)). We observe non-Gaussian (L-shaped) distribution of the events in the semi-log scale [52], that clarifies the occurrence of events with heterogeneously distributed sizes and hence of the extreme events.
In a similar fashion, we observed extreme events while considered different numbers of attractive or repulsive coupling zones (figures not shown here) and the inclusion of repulsive coupling zones are identified as necessary in order to realize extreme events.

Discussions and conclusions
In most of the previous works [50][51][52][53][54][55], the studies on extreme events are done either by considering only single isolated dynamical system or by considering static dynamical network formalisms and the possible temporal evolution of the network itself are rather ignored. But, dynamical systems in nature are rarely isolated and more importantly, in most of the cases interacting dynamical units in physical, biological and social networks undergo through time-varying connections instead of being coupled via static regular topologies. In contrast, we here have explored the phenomenon of synchronization together with its deflection towards the emergence of extreme events, based on fundamental principles of motion and network science, in a time-evolving model of interaction among moving agents, where each individual is moving in a two-dimensional plane with some preexisting coupling zones. The agents are allowed to interact with each other only when they move into the same coupling zone. The study of extreme events [71] have received less attention, although discovery of definitive extreme events indicating threshold is highly desirable for its scientific understanding, particularly for the problems, such as earthquakes, epileptic seizures and social dynamics, where the governing equations are unknown. Extreme events on complex networks play an important role due to its wide possible applications on dynamical processes. The analysis of extreme events can also help in the management of disease, for example to prevent sudden outbreaks. For instance, our study may be beneficial to avoid extreme events in many humanmade systems, e.g. reduction of overusage of antibiotics, which causes worsening medical cost and mortality especially for life-threatening bacteria infections [72]. Chen et al already showed that hysteresis can appear as an unprecedented roadblock for the recovery of vaccination uptake [73]. Our perceived study may be also useful to understand this puzzling phenomenon too.
In order to authorize our results, we successfully dealt with two well known chaotic systems, Lorenz and Rössler oscillators. The emergence of synchronization in this described time-varying complex network is thoroughly investigated in presence of only attractive coupling zones. Parameter regions for the synchronization state have been figured out by varying modulus velocity v of the moving oscillators for different number of coupling zones m. Numerically, it is shown that the modulus velocity v of the moving agents plays a crucial role in order to obtain synchronization of the system. The interval of interaction strength that induces synchrony has also been estimated analytically. Besides, extreme event like behavior deflecting the synchronization manifold is observed in the network whenever the repulsive coupling zone with effective strength has been introduced in the physical space. Non-Gaussian distribution of the event size is identified that explains the emergence of the extreme events.
In this paper, extreme events are considered as those short lasting events characterized by • unpredictable appearance with respect to time due to random mechanism with very low probability, • having amplitudes higher than H S =M+8σ, and • they appear much more often than they would according to Gaussian statistics [74].
The extreme value distribution obtained through block extrema approach demands an extensive sequence of data to conclude robust interpretation. This complication is somehow controlled to some extent using peak over threshold approach by considering a suitable extreme event indicating threshold depending on their probabilities of occurrence or on the point where they have potential consequences. Next, we provide a proper justification for the choice of the threshold value that distinguishes extreme events from the usual intermittency in synchrony, and with the help of this threshold value, the occurrence of the extreme events are established. We have also been able to map the dynamical states of synchronization and extreme events in the parameter plane of the attractive and repulsive interaction strengths, which points out to the fact that mere inclusion of repulsive zones may preserve stable synchrony region or induce intermittent behavior or even purely desynchronization region, but not necessarily extreme events. To obtain extreme events like phenomena, the proper choice of both coupling parameters is crucial. This also sheds light to the fact that proper choice of both coupling parameters can help the system push away from extreme events, which helps to mitigate extreme events.