Understanding the enhanced synchronization of delay-coupled networks with fluctuating topology

We study the dynamics of networks with coupling delay, from which the connectivity changes over time. The synchronization properties are shown to depend on the interplay of three time scales: the internal time scale of the dynamics, the coupling delay along the network links and time scale at which the topology changes. Concentrating on a linearized model, we develop an analytical theory for the stability of a synchronized solution. In two limit cases the system can be reduced to an"effective"topology: In the fast switching approximation, when the network fluctuations are much faster than the internal time scale and the coupling delay, the effective network topology is the arithmetic mean over the different topologies. In the slow network limit, when the network fluctuation time scale is equal to the coupling delay, the effective adjacency matrix is the geometric mean over the adjacency matrices of the different topologies. In the intermediate regime the system shows a sensitive dependence on the ratio of time scales, and specific topologies, reproduced as well by numerical simulations. Our results are shown to describe the synchronization properties of fluctuating networks of delay-coupled chaotic maps.


Introduction
In many interacting systems, the transmission time for information exceeds the time scale of the internal node dynamics. Hence, delay-coupled networks are relevant in a variety of fields, including coupled optical or opto-electronic systems, communication and transportation systems, social networks and biological networks as gene regulatory and neural systems. For example, in the brain a coupling delay between interacting neurons arises from the conduction time of an electric signal along the a e-mail: o.dhuys@aston.ac.uk we show that the dynamics depends sensitively on the ratio of time scales and the properties of the temporal topologies.
The paper is organized as follows: in Section 2, we present the modelling equations and apply the master stability formalism to the corresponding linearized model. In Section 3, we discuss the behavior of the system in the fast and slow network approximation, respectively. The interplay of time scales is largely discussed as well by paying attention to the parity effect in the relation between τ and T n . Sections 4 and 5 present our numerical simulations both for linear and nonlinear systems, respectively. The last section is devoted to conclusions. The more elaborate analytical calculations are presented in the Appendix.

Modelling equations
We start from a general model of N identical scalar elements coupled with interaction delays,ẋ with x i ∈ R. The coupling topology is modelled by a time-varying N × N adjacency matrix A(t), whose rows add up to one to ensure the existence of a permutation symmetric state. The coupling delay τ is constant over the links.
To determine the stability of a symmetric solution x 1 (t) = x 2 (t) = . . . = x N (t) ≡ x(t) (i.e. a symmetric fixed point, an in-phase oscillatory solution or a chaotic state in complete synchronization), the modelling equation is linearized around the symmetric solution x(t). A small perturbation δx i (t) = x i (t) − x(t) evolves according to, where f and g are the derivatives of the functions f and g, respectively, evaluated along the symmetric solution x(t).
We consider the simplest case with constant coefficients: the first term f (x(t)) ≡ −λ 0 represents an "internal" decay rate of the nodes; if the nodes are chaotic, it reduces to the (opposite) sub-Lyapunov exponent [14]. The term g (x(t − τ )) ≡ κ represents the coupling strength. After renaming the variables δx i ≡ x i , the linearized model can be rewritten asẋ where x(t) = (x 1 (t), . . . , x N (t)) T . In a constant network, it is straightforward to solve this system analytically by calculating the master stability function. Evaluating equation (3) along the eigenvectors v k of the adjacency matrix A, one findṡ with γ k the eigenvalue associated to the eigenvector v k . The system evolves then exponentially with a rate given by the master stability function λ(γ k ). In the limit of long delays τ λ −1 0 , and in the absence of a strongly unstable solution, λ 0 > 0, the set of exponential solutions of equation (4) {λ(γ k )} can be written as a pseudocontinuous spectrum [39] 1132 The European Physical Journal Special Topics This result applies for steady states, or for simple chaotic systems with a constant slope, as the Bernouilli map. However, it is also a first order approximation for chaotic systems [40] and reproduces the scaling properties of the spectrum of Lyapunov exponents of chaotic systems with time-delay τ and eigenvalues of magnitude |γ| [14,15].
Here, we consider a coupling matrix A(t) that is not constant: it changes discontinuously after a network time T n , running through a sequence of network topologies as A 1 , A 2 , . . . Thus, the system is non-autonomous, with time-dependent parameters [41]. However, in the following we show that, in certain limits, the synchronization properties of the system under a fluctuating topology can be described with a constant "effective" coupling topology A eff , allowing for to calculate a master stability function for nodes coupled with a time-varying topology.

Fast network approximation
In instantaneously coupled networks, it is well known that, if the network changes fast enough, the effective network is the average network over time [37]. This so-called "fast switching approximation" is valid as well in delay-coupled networks. Indeed, if T n λ −1 0 , one can approximate so that at t = t 0 + M T n , up to first order in T n , we retrievė This leads to an "effective" adjacency matrix We illustrate this result for a small network of ten coupled nodes, that alternates regularly between two topologies. The topologies A 1 and A 2 are randomly chosen, non-commuting small-world networks with 10 nodes, with normalized row sum. Their adjacency matrices are plotted in Figure 1.
According to the fast switching approximation, the nodes evolve exponentially. The decay or growth rate is given by the most unstable solution of equation (5), using the effective adjacency matrix A eff given by equation (8). For the exponential decay rate, we find Along the synchronization manifold, we have γ eff = 1 by construction, and hence the dynamics on the synchronization manifold is the same as in a constant network. The transverse stability, which determines the synchronization properties of the network, is determined by the second largest eigenvalue of the arithmetic mean network.
The transverse evolution is measured by the variance over the nodes. We have so that the transverse decay rate (TDR), λ, can be estimated from the evolution of the variance as We compare the theoretical decay rate (Eq. (9)) with the numerically calculated evolution of the variance (Eq. (11)) in Figure 2. As initial function, we used white, uniformly distributed noise. We find that the agreement between theory and simulations is excellent.

Slow network approximation
Let us consider the situation in which the network time T n is equal to the coupling delay τ , and both are larger than the instantaneous decay rate λ −1 0 of the nodes, λ −1 0 T n = τ . In this case, the coupling is constant during each delay interval, and it is straightforward to integrate equation (3). We consider an arbitrary initial function x 0 (t), t ∈ [−τ, 0]; decomposing into its Fourier components with ω n = 2πn/τ , one finds for the evolution of the nth mode during the first delay interval x which is solved by Since the terms proportional to e −λ0t become negligible after a short transient of order O(λ −1 0 ), we can approximate the general solution during the first delay interval as x 1 (t) = n x 1n e iωnt , with Note that, for a constant network, this corresponds to the decay rates given by pseudocontinuous spectrum equation (5).
Repeating this procedure for M time delays, and thus M alternations of the topology, one finds a general solution x M (t) = n x M n e iωnt , with Thus, we retrieve an "effective" adjacency matrix This prediction is verified numerically in Figure 3. Again, we simulated the model equation (3) for ten nodes, with the coupling configuration alternating regularly between the networks shown in Figure 1. Our theory predicts a TDR from equation (9), where in the slow network limit γ eff is given by the eigenvalues of the effective adjacency matrix equation (17). Figure 3 compares the variance of the nodes for white noise initial conditions (black curve) with the theoretical prediction. Also in this case the theoretical prediction provides an excellent approximation for the numerical decay rates. The difference between the best linear fit of the simulations and the theoretical decay rate is around 1%. Comparing Figures 2 and 3, it is clear that the arithmetic mean network synchronizes faster than the geometric mean network. This is the case for most pairs of stochastic matrices. In the fast switching approximation typically the network synchronizes faster than both topologies between which it alternates, and thus the network fluctuations can be said to enhance synchronization. In contrast, one usually finds a TDR that is in between the decay rates of both topologies in the slow network limit. Thus, an appropriate choice of topologies allows to control the synchronisation properties, also in the slow network limit.

Interplay of time scales
The linear model given in equation (3) is difficult to solve in the most general case. Nonetheless, it is possible to find exact solutions for the special case that the delay time τ is a multiple of the network time T n . The details of this calculation are outlined in the Appendix. We summarize the main results in the following. As a first simplification, we consider a regular alternation between two topologies, with respective adjacency matrices A 1 and A 2 . Thus, the system repeats the cycle exactly every 2T n . The results can be generalised to a periodic sequence of topologies A 1 , . . . , A M . In the direct integration, we have assumed a constant initial function x(t) = x 0 for t < 0, motivated by the fact that for constant coupling, and in the slow network limit (see Eqs. (5) and (16)) the zeroth Fourier mode is the least stable, and consequently, determines the stability. Similarly to the slow network limit, it is then possible to  integrate equation (3) for consecutive delay intervals, concatenating the output of an interval with the input for the next.
The evolution of the network depends on the interaction of the two coupling topologies, as, during an interval of duration T n with constant topology A 1 , the history term may depend on A 1 and A 2 or switch between both coupling configurations at a certain point. Thus, we expect a periodic (parity) effect in the overall TDR with respect to mod (τ, 2T n ). This is confirmed analytically, as we find a significantly different evolution depending on whether the coupling delay τ is an even or odd multiple of the network time T n .

Asymptotic behavior for τ = 2M T n
If the coupling delay τ is an even multiple of the network time T n , τ = 2M T n , M ∈ N, our analytic calculation (outlined in Appendix A.1) indicates that the arithmetic mean network, with an effective adjacency matrix A given by determines the asymptotic decay rate, independent of T n , even though the fast network limit does not apply. However, the time to reach this asymptotic decay rate increases with the network time. While this decay rate is reached quickly in fast switching networks, the transient time can be estimated as 4λ 0 T n τ for slowly changing topologies. Thus, transient dynamics can play an important role as the network time increases.
In particular, for relatively slow networks λ 0 T n 1, this initial transverse exponential evolution rate is approximated by that of the most unstable network, |γ eff | = max(|γ 1 |, |γ 2 |), with γ 1 and γ 2 the respective maximal transverse eigenvalues of A 1 and A 2 . Hence, it is possible that, although the synchronized state is asymptotically linearly stable, a perturbation may be initially amplified, and drive a nonlinear system out of the basin of attraction of the synchronized state.
To illustrate our analytical results, we have simulated a network with an even ratio of coupling delay τ and the network time T n , shown in Figure 5. In order to simplify comparison between the dynamics for even and odd ratios of time scales, we consider a regular alternation between topologies with commuting adjacency matrices. Our networks have 11 nodes, which evolution is modelled by equation (3), while the alternating networks are chosen to be a clockwise and a counterclockwise unidirectional ring, each coupled to its two nearest neighbors. A sketch of these topologies are shown in Figure 4. We chose τ = 20T n . For our choice of T n = 5, one observes that the asymptotic decay rate is only reached after ∼ 4λ 0 T n = 20 delay intervals. The theoretical TDR is given by equation (9), with γ eff in the long time given being the maximal transverse eigenvalue of the effective adjacency matrix (Eq. (18)). Once more, the agreement between theoretical approximation and numerical simulations is good.

Asymptotic behavior for τ = (2M + 1)T n
If the delay time τ is an odd multiple of the network time T n , the topologies interfere in a different way: the adjacency matrix A 1 enters the history term multiplying the coupling matrix A 2 .
If λ 0 T n is small, we retrieve the fast switching approximation. For slower networks λ 0 T n 1 it is only possible to compute the long-time evolution, under the condition that both adjacency matrices commute. This calculation is outlined in Appendix A.2. In this case, the asymptotic TDR is given by equation (9), with an effective eigenvalue γ eff given by with γ 1 and γ 2 the respective transverse eigenvalues of A 1 and A 2 along the same eigenvector, andγ = √ γ 1 γ 2 . Thus, the effective eigenvalue approaches the geometric mean eigenvalue γ eff = √ γ 1 γ 2 for large T n . Note that the first order correction in equation (19) is proportional to (λ 0 T n ) −1 for slow networks. In particular, in the case that |κγ| = λ 0 , this dependency explains the power law behavior λ ∝ T −1 n demonstrated in Figure 9. In contrast to the even ratio of time scales, there is no pronounced transient behavior, as also in the slow network limit the initial decay rate is approximated by the geometric mean network. We compare these theoretical results to numerical simulations in Figure 6. The network equation (3) of 11 nodes, alternates regularly between the commuting topologies shown in Figure 4, with τ = 25T n . We initialized our system again with white noise. We find again excellent agreement between theory and simulations. Note the considerably different TDRs in Figures 5 and 6, while the network times do not differ much (T n = 5 and T n = 4, respectively), illustrating the sensitive dependence of the system on the ratio of time scales.

Simulations for varying network time
To explore the synchronization properties in the full range of network times T n , we have performed numerical simulations of the aforementioned linear system equation (3) with delayed interactions: In order to easily compare with the analytic results, we consider A(t) to be a discontinuous matricial process which proceeds through the alternation of two matrices, A 1 and A 2 , every T n . The initial function in our simulations is constant x(t) = x 0 for all t < 0 (in contrast to the numerical results shown in Figs. 2-6).
Unless otherwise mentioned, our network consists of three nodes, N = 3. Our first choice for A 1 and A 2 is the cyclic choice, which is illustrated in Figure 7. The system alternates between the two different cyclic directed graphs with three nodes, i.e.: Note that these matrices commute and, moreover, they are inverse, i.e. In a second set of simulations, the commuting choice, we use the topologies A 1 and A 2 as follows: In this case, all theoretical results can be applied, but the matrices are no longer inverse.
A third possible choice is the random choice for A 1 and A 2 . The random adjacency matrices employed in our simulations are constructed the following way: the entries A ij with i = j are assigned a uniformly distributed random in the interval [0, 1]. Afterwards, the adjacency matrices are row-normalized in order to respect the stochasticity condition. For the employed examples, we check the matrices to be non-commutative, and compute the spectral gap of the arithmetic and geometric average.
Our simulations are based on a simple Euler integrator with a short enough timestep (∆t = 10 −2 ), endowed with a memory structure from where the required timedelayed data can be obtained in an efficient way. The initial condition is chosen randomly, each unit x i drawn from a uniform distribution in [0, 1], but, in contrast to the simulations presented in the previous section, constant in the interval [−τ, 0]. Figure 8 shows a few illustrative histories when two cyclic matrices are alternated along equation (3), using λ 0 = κ = 1, τ = 100 and two values for the network switching time: T n = 10 and T n = 4, and a fixed random initial condition. The top panel shows the evolution of the first node, x 1 (t). For T n = 10 the delay time is an even multiple of the network time, while for T n = 4 it corresponds to an odd multiple. As predicted by equation (18) for our choice of parameters, in all the cases shown, the different components approach a synchronized state.
The second panel of Figure 8 shows the evolution of the standard deviation between all components, along with our estimate of the TDR. The slower decay for T n = 4 results from the difference between τ being an odd and an even multiple of T n . Notice that the precise measurement of the decay rate is hindered by the strong oscillations, whose periodicity is given by the interaction delay τ .
Estimating the TDR can be a complicated numerical task. In order to perform this in a robust way, we have adopted the following algorithm: (a) a set of temporal windows are selected within the simulation range, with random lengths and random initial points, (b) for each of them, the TDR is estimated using a least-squares exponential fit of σ versus time and (c) the average for the different outcomes provides our best estimate for the TDR, while its deviation provides a natural approximation to the error. The main reason to use this technique is that it tends to discard the initial transient region, which is one of our biggest concerns. Yet, this technique is still amenable to future improvements.  (3) for the cyclic choice of A1 and A2, given by equation (20). Top: evolution of a single component. Bottom: evolution of the deviation. We use N = 3, λ0 = κ = 1, τ = 100, Tn = 4 and Tn = 10. The straight lines correspond to the best exponential fit, allowing us to estimate the transversal decay rate.
Using this method, we have computed the TDR in a variety of cases. Our theoretical estimates for the TDR are provided by equation (9), where the value of γ eff corresponds to the arithmetic mean of the matrices in the fast-switching regime λ 0 T n 1 and to the geometric mean in the slow-switching regime T n ∼ τ . The top panel of Figure 9 shows our estimates for the TDR as a function of network time T n , using τ = 100 and λ 0 = κ = 1. Each panel is devoted to a different choice for the A 1 and A 2 adjacency matrices: random (a), cyclic (b) and commuting (c). The horizontal lines mark our theoretical estimate in the fast-switching (continuous line) and slow-switching (dashed line) regimes. For the random case chosen (top panel), where matrices A 1 and A 2 do not commute, we observe a good correspondence to the fast switching limit for small values of T n and a correspondence to the slow network case for a range of network times T n ≈ τ , but there is no general trend visible. This may be due to the proximity of both limits.
In the commuting cases (central and bottom), however, there is a clear trend visible. Again, for fast fluctuations the transverse stability is well approximated by the fast switching approximation. The TDR evolves in both cases towards the slow network limit as T n ∼ τ , the general trend is well described by the decay rate for odd multiples (Eq. (19)). For the commuting and random choices, the estimate for the decay rate decreases further as T n τ . This can be explained by the fact that the dynamics operates on a time scale τ and thus adiabatically follows the temporary network, which evolves on a slower time scale in this limit.
For our choice of parameters, the slow-switching regime in the cyclic case provides a TDR ∼ 0, as shown in our data. In this case, the decay from the fast to the slow switching regimes corresponds approximately to a power-law, with TDR λ ∼ T −1 n , as shown in the inset of Figure 9 (middle). This power law is explained by the odd multiple decay rate. Indeed, combining equations (19) and (9), we can deduce,

1142
The European Physical Journal Special Topics where we usedγ = −1, our parameter choice κ = λ 0 = 1 and (λ 0 T n ) −1 is considered small. The bottom panel shows the same data, but with a restructured abscissa: the new independent variable is τ /T n ≡ M . In this case, the parity oscillations in the fixed and the cyclic case are clearly visible. In particular, we observe synchronizing resonances corresponding to the even values of M , in agreement to the theoretical predictions. In the commuting case, the TDR at these resonances is particularly well described by the arithmetic mean, also in agreement with the theory (Eq. (18)).
In Figure 10, we check whether this evolution from arithmetic mean effective topology in the fast network limit to geometric mean in the slow network limit holds for larger networks with random (non-commuting) adjacency matrices, A 1 and A 2 . We simulated networks of sizes N = 20, N = 40 and N = 60 (see Fig. 10 from bottom to top, respectively). Again, the fast and slow regime approximations to the TDR are marked with dashed and continuous horizontal lines, respectively. We observe a general, qualitative agreement between the theoretical prediction and the numerical data, even though the matrices are non-commuting in this case. Moreover, we observe the same strong parity oscillations.

Synchronization of delay-coupled chaotic maps
We illustrate our results for linear, time-continuous systems in a network of chaotic maps. Similar as in previous work [38], we consider a time-varying network of N delay-coupled Bernouilli maps with a discrete time evolution modelled by where x = (x 1 , . . . , x N ) is the vector of the Bernouilli unit states, f (x) = ax mod 1, with a > 1. The network adjacency matrix A(t) belongs to a sequence {A 1 , A 2 , · · · } of adjacency matrices randomly sampled from a Newmann-Watts small-world network ensemble [43]. Instances of this ensemble are directed small-world networks similar to the standard Watts-Strogatz networks, with a directed outside ring and each node has a probability p of establishing a new "shortcut" link with a randomly chosen node. The difference with standard small-world networks is that here shortcuts are added without removing the corresponding ring links, thus keeping the outside ring fixed and ensuring the connectivity of the ensemble networks. Like in the preceding sections, the connectivity switches instantly every T n timesteps. The system's evolution to our linear model can be compared to the timecontinous model equation (3), by identifying the (opposite) sub-Lyapunov exponent and the internal decay rate. The coupling strength can simply be mapped, Notice that for this system, the instantaneous decay rate λ 0 is dependent on the coupling strength . For chaotic maps, the TDR from the linear system coincides with the so-called synchronization, or transverse Lyapunov exponent (SLE). That is, the rate governing the evolution of small perturbations around the synchronized state, x i (t) = x(t), ∀i. In a fixed network of Bernoulli units, the SLE is, similar the linear system (Eq. (5)), computed as follows [44]: where γ 2 is the adjacency matrix' second largest eigenvalue. We simulated the dynamics, equation (23), for several values of the coupling strength , with a fixed p and for two network switching times: T n = 100 and T n = 10. The comparison between the resulting numerical SLE and the effective SLE corresponding to the fast and slow network approximations is plotted in Figure 11. The first case, T n = 100, corresponds to the slow regime conditions λ −1 0 T n = τ , and the numerical SLE values follow the slow network approximation closely. The second case, T n = 10, crosses over from the fast network regime λ −1 0 T n at low values to the slow network regime at values of ∼ 1. This is visible in the measured

Conclusion
We have developed an analytical formalism within the linearized limit in order to understand synchronization phenomena in the case of delay-coupled networks with a fluctuating topology, studied previously by us in [38]. Considering a regular alternation of the topology between different configurations, we have studied both the case of fast and slow fluctuations. Based on the master stability function approach, we derive the stability of a symmetric state based on the alternation of topologies and the interplay between the timescales involved: the internal timescale of the network nodes, the coupling delay, and the characteristic fluctuation time of the network. In two different regimes, we have derived an "effective" coupling topology: when the network fluctuations are faster than the internal time scale, the fast switching approximation can be extended to time-delay systems, and we find that the effective network is given by the arithmetic average. When the coupling delay and network fluctuation time are similarly large compared to the internal time scale of the network nodes, the effective network is given by the geometric average over the different topologies.
As the network time varies between the fast and slow limit, we find analytically a parity effect: the synchronization properties depend on the ratio of the delay time and the network time. In particular, if the adjacency matrices of the different topologies commute, one retrieves the fast network limit if the ratio between both time scales is even, while for odd ratios the behavior evolves from the fast to the slow network limit as the network time increases. Complementing these results with numerical simulations, we (broadly) recover this evolution from fast to slow network limit, from arithmetic to geometric mean network with increasing network time. This trend is visible for all network sizes, and for non-commuting and commuting topologies as well, but the agreement with the analytical theory is better in the latter case, and for larger networks. The parity effect is shown to be a universal feature for a regularly alternating topology.
Finally, we compare our theoretical results to the synchronization properties of an ensemble of delay-coupled chaotic maps. As the coupling strength, and thus the internal time scale varies, we show that the synchronization boundaries shift between the geometric and arithmetic means of the ensemble, confirming the linear theory. It is worth stressing however, that the effective synchronization properties of a nonlinear network might differ from this linear theory, due to long transient times and wild oscillatory behavior in the linear regime.
Future extensions of the research might be the study of other network ensembles, such as random Erdös-Rényi graphs, scale-free networks or even more complicated graphs of multiplex type with further application to real world problems concerning transport and energy issues or problems of supply networks in the general context of "Smart cities". Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Appendix A
We outline the analytical calculations leading to the results presented in Section 3.3. We directly integrate equation (3), using constant initial conditions x(t) = x 0 for t < 0, and a regular alternation of the coupling topology between adjacency matrices A 1 and A 2 after a network time T n . During the first delay interval, the dynamics then is modelled aṡ The European Physical Journal Special Topics After an initial transient of the order of O(λ −1 0 ), the resulting dynamics during the delay interval is 2T n -periodic: the system evolves exponentially between t = kT n and t = (k + 1)T n , when it changes direction. The values at these turning points are denoted κ λ0 x 1A and κ λ0 x 1B respectively. During the first half cycle of its periodic motion, the solution x A1 (t) of equation (A.1) reads During the second half cycle of its periodic motion, the solution x B1 (t) of equation (A.1) reads Assuming a continuous periodic motion x B1 (T n ) = κ λ0 x 1A and x A1 (T n ) = κ λ0 x 1B , we find a solution for the turning points Note that in the limit λ 0 T n 1, we find as a leading order approximation which corresponds to the fast switching approximation. However, in order to find the long term evolution, one needs to integrate over the next delay intervals.
A.1 Solution and convergence rate if τ = 2M T n Let us focus on the case when τ is an even multiple of the network time T n , τ = 2M T n , with M ∈ N arbitrarily large. By writing equation (3) for the second delay interval we obtainẋ A,B2 (t) = −λ 0 x A,B2 (t) + κA 1,2 x A,B1 (t) .

(A.21)
Thus, one finds the slowest decay rate when choosing the sign ofγ such that its argument is closest to the argument of γ 1 + γ 2 .