# The Impact of Competing Time Delays in Coupled Stochastic Systems

###### Abstract

We study the impact of competing time delays in coupled stochastic synchronization and coordination problems. We consider two types of delays: transmission delays between interacting elements and processing, cognitive, or execution delays at each element. We establish the scaling theory for the phase boundary of synchronization and for the steady-state fluctuations in the synchronizable regime. Further, we provide the asymptotic behavior near the boundary of the synchronizable regime. Our results also imply the potential for optimization and trade-offs in synchronization problems with time delays.

###### keywords:

Stochastic synchronization and coordination, Time delays, Scaling###### Pacs:

05.40.-a, 05.45.Xtsort&compress

## 1 Introduction

Coordinating, distributing, and balancing resources in networks is a complex task as these operations are very sensitive to time delays Saber_IEEE2007 (); Hunt_PRL2010 (). To understand and manage the collective response in these coupled interacting systems, one must understand the interplay of stochastic effects, network connections, and time delays. In synchronization, coordination, and consensus problems in coupled interacting systems Arenas_PhysRep2008 (); GK_Science2003 (); Nishikawa_PNAS2010 (); Saber_IEEE2007 (); Hunt_PRL2010 (), individual units attempt to adjust their local state variables (e.g., pace, load, orientation) in a decentralized fashion. They interact or communicate only with their local neighbors in the network, often with explicit or implicit intention to improve global performance. Applications of the corresponding models range from physics, biology, computer science to control theory, including synchronization problems in distributed computing GK_Science2003 (); Guclu_PRE2004 (); Guclu_FNL2005 (), coordination and control in communication networks GK_PRE2007 (); Hunt_PRL2010 (); Johari_IEEE2001 (); Saber_IEEE2004 (); Saber_IEEE2007 (); Jadbabaie_IEEE2006 (); Scutari_IEEE2008 (), flocking animals Reynolds_ACM1987 (); Vicsek_PRL1995 (); Cucker_IEEE2007 (), bursting neurons Atay_PRL2003 (); Gassel_FNL2007 (); Chen_EPL2008 (); Chen_PRE2009 (), and cooperative control of vehicle formation Fax_IEEE2004 ().

In this Letter, we study the impact of competing, finite non-zero time delays in stochastic synchronization or coordination problems, which are present in most real communication, information Hunt_PRL2010 (); Johari_IEEE2001 (); Saber_IEEE2004 (); Saber_IEEE2007 (); Huberman_IEEE1991 (); Strogatz_PRE2003 (); Chen_PhysA2004 (), and biological systems Hutchinson_1948 (); May_Ecology1973 (); Ji_FNL2008 () including neurobiological networks Gassel_FNL2007 (); Chen_EPL2008 (); Chen_PRE2009 (). (Throughout this Letter, we use the terms coordination and synchronization synonymously.) Delays can be attributed to both non-zero transmission times between the nodes and to non-zero finite times it takes to process (possibly including cognitive delays) and execute the desired action at the nodes. Here, we investigate the importance and impact of these two types of delays in a simple synchronization problem in noisy environment with two linearly coupled nodes.

Singularities in critical phenomena and phase transitions Dorog_RMP2008 (), which are often present in coupled interacting systems consisting of a large number of nodes , are typically associated with progressively more eigenvalues of the coupling operator (e.g., the Laplacian) getting arbitrarily close to zero. Strictly speaking, these singularities are exhibited only by systems approaching the thermodynamic limit (where the density of eigenvalues does not vanish sufficiently fast, or itself becomes singular about zero in the limit). For example, in spatially-embedded physical systems these singularities are typically exhibited by the relevant response functions and fluctuations in the long-wavelength limit GK_Science2003 (). In complex networks Barab_sci (); Watts_Nature1998 (); BarabREV (); MendesREV () these singularities can be suppressed as a result of sufficient amount of randomness in the connectivity pattern GK_Science2003 (); Barahona_PRL2002 (); Kozma_UGA2004 (); Kozma_PRL2004 (); Kim_PRL2007 ().

In contrast, the instability governed by time delays is associated with a single mode exceeding a threshold value (in a simple case, associated with the eigenmode of the network Laplacian with the largest eigenvalue Hunt_PRL2010 (); Saber_IEEE2004 ()). Therefore, the underlying instability is present even in the simplest network with two nodes (=). Here we focus on a two-node network, which qualitatively captures the generic features of the coordination behavior when the delays are present due to both transmission between nodes and processing/execution at each node. For simplicity, we will refer to this instability as “critical”, even though this singularity does not require infinitely many degrees of freedom. In networks, consisting of a large number of nodes, the effect of time delays is not qualitatively different, but can be “amplified” by heterogeneous (or scale-free) Barab_sci (); BarabREV (); MendesREV () connectivity patterns: in the case of uniform time delays Hunt_PRL2010 (); Saber_IEEE2007 (); Saber_IEEE2004 (), the effective coupling of the most relevant singular mode is the largest eigenvalue of the coupling operator, which itself can diverge with the system size Arenas_PhysRep2008 (); Fiedler_1973 (); Anderson_1985 (); Mohar_1991 (), severely limiting synchronizability and coordination.

## 2 A Stochastic Model with Local and Transmission Delays

Differential equations with delays Bellman_1963 () describing complex systems have a long history, originally motivated by the emergence of business and economics cycles Kalecki_1935 (); Frisch_1935 (); Hayes_1950 (), and also naturally appearing in the context of stability of ecological systems, in models in population dynamics, and in game theory Hutchinson_1948 (); May_Ecology1973 (); MacDonald_1978 (); Ruan_2006 (); Yi_JTB1997 (); Miekisz_2007 (). There have been recent works combining stochastic differential equations with delays Kuchler_MCS2000 (); Amann_PhysA2007 (); Ibanez_FNL2008 () with applications ranging from population dynamics, epidemiology, and immunology to cell kinetics and finance Bocharov_JCAM2000 (); Tian_JCAM2007 ().

Here we consider a model for local coordination where time delays are attributed to two separate origins: one is the transmission between the two nodes, the other is processing the information and executing the action at each node, denoted by and , respectively. We investigate the simplest stochastic model where the coordination or synchronization attempt between the two nodes, in terms of the relevant state variables , is captured by linear relaxation

(1) |

Here, is delta correlated noise with and with noise intensity , . is the coupling strength between the two nodes. For initial conditions, we use for .

To simplify notation we introduce and (). Further, since we are interested in the synchronization (or coordination) between the two nodes, we focus on the relative difference which is governed by

(2) |

where and . The special case of the above equation has been investigated in our earlier work Hunt_PRL2010 (). Here, we study the impact of both types of delays, corresponding to the general case . Our quantity of interest is , capturing the relative deviation of the relevant state variables on the two nodes. By definition, the system is synchronizable if the fluctuations reach a finite steady state, . In the absence of time delays () one immediately finds Gardiner_1985 (), i.e., the system is synchronizable for any . Further, the stronger the coupling, the better the synchronization: is a monotonically decreasing function of .

Next we study and analyze the case with time delays. Employing standard Laplace transform Hunt_PRL2010 (); Bambi_JEDC2008 (), one can immediately write the formal solution for Eq. (2)

(3) |

where , , are the zeros of the characteristic equation

(4) |

on the complex plane. Then for the noise-averaged fluctuations one finds

(5) | |||||

where in the last expression of the above equation we introduced the scaled variables and . From Eq. (4) and from the definition of these scaled variables it is evident that are the solutions of the scaled characteristic equation

(6) |

and consequently, the solutions depend only on , i.e., . From the structure of the above characteristic equation it follows that if is a solution of Eq. (6) so is its complex conjugate . From Eq. (5) it is clear that synchronization can only be achieved if for all . To identify the boundary of the region of synchronizability, one has to find the solution(s) with a vanishing real part, i.e., with Frisch_1935 (); Hayes_1950 (); May_Ecology1973 (); Ruan_2006 (). Elementary analysis yields and

(7) |

Thus, for a fixed , the system is synchronizable if . Results obtained by numerically integrating Eq. (2) note_numerical () together with the analytic expression Eq. (7) are shown in Fig. 1.

We will discuss the phase diagram and some limiting cases in terms of the original variables, the local delay and the transmission delay , in the final section.

## 3 Scaling and Asymptotics in the Steady State

Now we turn to analyzing the steady-state fluctuations, in particular, their scaling behavior in the synchronizable regime, . Here, the fluctuations remain finite, and in the steady state () from Eq. (5) one obtains

(8) |

where

(9) |

is the scaling function for the steady-state fluctuations. [Recall that are the solutions of the scaled characteristic equation Eq. (6).] Thus for a given ,

(10) |

where in our notation, we suppressed the dependence to highlight the scaling behavior of the fluctuations, valid for each separately. Figure 2 shows the steady-state fluctuations before (a) and after (b) scaling, and demonstrates the data collapse for the scaled variables according to Eq. (10).

The scaling function is a non-monotonic function of its argument, diverging at and , and exhibiting a single minimum between these points [Fig. 2(b)]. This non-monotonic feature of the scaling function with a single minimum between is present for all [Fig. 3]. Thus, for fixed non-vanishing and finite delays, there is an optimal value of the coupling constant for which the steady-state fluctuation attains its minimum value. For stronger couplings, the overall coordination between the two nodes weakens, and for , it completely deteriorates.

Next, we briefly discuss the asymptotic behavior of the scaling function near the boundaries of the synchronizable regime. The fluctuations of diverge at the end points of this interval [as at least for one , ], indicating the breakdown of synchronization. Near these endpoints, the sum in Eq. (9) is dominated by the term(s) where Hod_2010 (). These are the solutions which have (negative) real parts with the smallest amplitude. As we show in Appendix A, to leading order,

(11) |

as , and

(12) |

as () with given in Appendix A [Eq. (24)]. From the numerical results [Fig. 2(b)] it is also apparent that the scaling function varies slowly between (and away from) the singular points, thus, can be reasonably well approximated Hod_2010 () throughout the full synchronizable regime by

(13) |

with also given in Appendix A [Eq. (28)].

Figure 2(b) and Fig. 3 show that the above approximate scaling function Eq. (13) (being asymptotically exact near the singular points) matches the numerical data very well. In particular, it captures the basic non-monotonic feature of the results obtained from numerical integration, exhibiting a single minimum

(14) |

in the interval.

As can also be seen in Fig. 3, the theoretical asymptotic behavior, captured by the approximate scaling function Eq. (13) becomes less accurate for small near . Other than lacking higher-order corrections to the asymptotic expressions, this is due in part to the time discretization in the numerical integration note_numerical (). For sufficiently small values, the condition will not hold, and deviations between the results of the time-discretized numerical scheme and those of the continuous system Eq. (2) will become more significant and noticeable.

## 4 Discussion and Summary

Having established the scaling theory for the phase boundary [Eq. (7)] and for the fluctuations [Eq. (13)], it is insightful to express our main findings explicitly in terms of the two types of delays appearing in the original formulation of the problem, Eq. (1). From Eq. (7), for the boundary of the synchronizable regime one immediately finds

(15) |

While explicitly expressing the critical line vs is prohibitive due to the implicit nature of Eq. (15), one can produce a plot for it numerically without any difficulties [Fig. 4].

We can gain further insight into the different impact of the two types of delays by considering two limiting cases. First, consider the case when , i.e., when the transmission delays are much larger than the local processing, cognitive, or execution delays. This is equivalent to the limit in our scaling expressions. From Eq. (7) one finds or . Thus, there is no singularity in the fluctuations for any finite provided that . Further, from Eq. (13) (with the coefficients given in Appendix A) for the steady-state fluctuations in the same limit we find

(16) |

In the other limiting case, , i.e., the transmission delays are much smaller than the local processing delays. This is equivalent to the limit in our scaling expressions. In this limit Eq. (7) reduces to or . The steady-state fluctuations approach

(17) |

provided that .

Figure 4 and Eqs. (16) and (17) highlight the subtle differences between the impacts of the two types of delays. We may regard the local delays as the dominant ones, in that as long as , there are no singularities for any finite , and increases linearly with as [Eq. (16)]. On the other hand, for every , there is a sufficiently large such that the fluctuations become singular. In particular, when the transmission delays are much smaller than the local processing delays, the fluctuations diverge as [Eq. (17)].

Inside the synchronizable regime, for fixed and (with the exception of note_min_limit ()), the steady-state fluctuations always exhibit a single local minimum as a function of the coupling constant [Eq. (14)]. This feature naturally presents scenarios for optimization and trade-offs. In particular, in real systems, the effective coupling constant between two interacting nodes can be controlled by the frequency of pairwise communications. This implies that too much communication can cause more harm than good, and further, there is an optimal rate of communications between the nodes which minimizes the size of the fluctuations. Also, as was already shown for the special case Hunt_PRL2010 (), in large network-coupled systems, decreasing connectivity may be beneficial if the system is too close or beyond the synchronization boundary.

In this Letter we only considered linear couplings between two nodes in the presence of noise and competing time delays. We established the phase diagram and the scaling theory for synchronizability, and provided the asymptotic behavior for the relevant scaling functions. Nonlinearities are undoubtedly important in real systems, and will likely further increase the complexity of the already rich phase diagram, such as recurring and alternating patterns of synchronizable and unsynchronizable regions Chen_EPL2008 (); Chen_PRE2009 (); Strogatz_PRE2003 ().

## Acknowledgements

We thank S. Hod for bringing our attention to his preprint and subsequently published work on the asymptotic method for the scaling function Hod_2010 (). This work was supported in part by DTRA Award No. HDTRA1-09-1-0049, by the Army Research Laboratory under Cooperative Agreement Number W911NF-09-2-0053, and by the Office of Naval Research Grant No. N00014-09-1-0607. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the Army Research Laboratory or the U.S. Government.

## Appendix A Asymptotic Behavior of the Scaling Function Near the Synchronization Thresholds

Here, we employ the method in Ref. Hod_2010 () to calculate the dominant contributions in Eq. (9) near the boundaries of the synchronizable regime. We assume that solutions of the characteristic equation

(18) |

change continuously with the parameter . Thus, if is a solution for , then for a small change in the parameter, , the corresponding solution can be written as . Substituting this into the characteristic equation, to lowest order we find

(19) |

For , there is a single solution with vanishing real part, , thus for small

(20) |

The dominant contribution for the scaling function as comes from the corresponding term in Eq. (9), to leading order yielding

(21) |

For [Eq. (7)], there is a pair of solutions (complex conjugates) with vanishing real parts . When is in the vicinity of (), to lowest order, these solutions behave as

(22) |

where . The dominant contributions for the scaling function as then come from the two terms in Eq. (9) when , yielding

(23) | |||||

where

(24) |

Finally, we obtain the approximate scaling function for the full interval, using some heuristics following Ref. Hod_2010 (). As observed from our numerical results [Fig. 2], the scaling function varies slowly between (and away from) the two singular points. Then, it can be approximated by

(25) |

In principle, the constant could be determined by matching the minimum value of the scaling function. Since it is not known analytically, instead we resort to the heuristics of Ref.Hod_2010 () where the constant is determined in such a way that it matches next-to-leading order corrections of the asymptotic behavior, e.g., near . To that end, we find the next-to-lowest order corrections to the solution of Eq. (18) in the vicinity of ,

(26) |

Keeping the relevant orders in the dominant term in Eq. (9), we obtain

(27) | |||||

In order to match this next-to-leading order correction as with the proposed approximate scaling function Eq. (25), one must have

(28) |

## References

- (1) R. Olfati-Saber, J.A. Fax, and R.M. Murray, Proc. IEEE 95, 215 (2007).
- (2) D. Hunt, G. Korniss, and B.K. Szymanski, Phys. Rev. Lett. 105, 068701 (2010).
- (3) A. Arenas et al., Phys. Rep.469, 93 (2008).
- (4) G. Korniss, M.A. Novotny, H.Guclu, Z. Toroczkai, P.A. Rikvold, Science 299, 677 (2003).
- (5) T. Nishikawa and A.E. Motter, Proc. Natl. Acad. Sci. U.S.A. 107, 10342 (2010).
- (6) H. Guclu and G. Korniss, Phys. Rev. E 69, 065104(R) (2004).
- (7) H. Guclu and G. Korniss, Fluctuation and Noise Letters 5, L43 (2005).
- (8) G. Korniss, Phys. Rev. E 75, 051121 (2007).
- (9) R. Johari and D. Kim Hong Tan, IEEE/ACM Trans. Networking 9, 818 (2001).
- (10) R. Olfati-Saber and R.M. Murray, IEEE Trans. Automat. Contr. 49, 1520 (2004).
- (11) A. Papachristodoulou and A. Jadbabaie, in Proc. of the 45th IEEE Conf. on Decision and Control, (IEEE, 2006) pp. 4307–4312.
- (12) G. Scutari, S. Barbarossa, and L. Pescosolido, IEEE Trans. Sig. Process. 56, 1667 (2008).
- (13) C.W. Reynolds, Computer Graphics 21, 25 (1987).
- (14) T. Vicsek et al., Phys. Rev. Lett. 75, 1226 (1995).
- (15) F. Cucker and S. Smale, IEEE Trans. Automat. Contr. 52, 852 (2007).
- (16) F.M. Atay, Phys. Rev. Lett. 91, 094101 (2003).
- (17) M. Gassel, E. Glatt, F. Kaiser, Fluctuation and Noise Letters 7, L225 (2007).
- (18) Q. Wang, Z. Duan, M. Perc, and G. Chen, Europhys. Lett. 83, 50008 (2008).
- (19) Q. Wang, M. Perc, Z. Duan, and G. Chen, Phys. Rev. E 80, 026206 (2009).
- (20) J.A. Fax and R.M. Murray, IEEE Trans. Automat. Contr. 49, 1465 (2004).
- (21) T. Hogg and B.A. Huberman, IEEE Trans. on Systems Science and Cybernetics 21, 1325 (1991).
- (22) M.G. Earl and S.H. Strogatz, Phys. Rev. E 67, 036204 (2003).
- (23) C. Li and G. Chen, Physica A 343, 263 (2004).
- (24) G.E. Hutchinson, Ann. N.Y. Acad. Sci. 50, 221 (1948).
- (25) R.M. May, Ecology 54, 315–325 (1973).
- (26) Lin Ji, Weiguo Xu, and Qianshu Li, Fluctuation and Noise Letters 8, L1 (2008).
- (27) S.N. Dorogovtsev, A.V. Goltsev, and J.F.F. Mendes, Rev. Mod. Phys. 80, 1275 (2008).
- (28) A.-L. Barabási and R. Albert, Science 286, 509 (1999).
- (29) D.J. Watts and S.H. Strogatz, Nature 393, 440 (1998).
- (30) R. Albert and A.-L. Barabási, Rev. Mod. Phys. 74, 47 (2002).
- (31) S.N. Dorogovtsev and J.F.F. Mendes, Adv. in Phys. 51, 1079 (2002).
- (32) M. Barahona and L.M. Pecora, Phys. Rev. Lett. 89, 054101 (2002).
- (33) B. Kozma and G. Korniss, in Computer Simulation Studies in Condensed Matter Physics XVI, edited by D.P. Landau, S.P. Lewis, and H.-B. Schüttler, Springer Proceedings in Physics Vol. 95 (Springer-Verlag, Berlin, 2004) pp. 29-33.
- (34) B. Kozma, M. B. Hastings, and G. Korniss, Phys. Rev. Lett. 92, 108701 (2004).
- (35) D.-H. Kim and A.E. Motter, Phys. Rev. Lett. 98, 248701 (2007).
- (36) M. Fiedler, Czech. Math. J. 23, 298 (1973).
- (37) W.N. Anderson and T.D. Morley, Lin. Multilin. Algebra 18, 141 (1985).
- (38) B. Mohar, in Graph Theory, Combinatorics, and Applications Vol. 2, edited by Y. Alavi, G. Chartrand, O.R. Oellermann, and A.J. Schwenk (Wiley, New York, 1991) pp. 871–898.
- (39) R. Bellman and K. Cooke Differential-Difference Equations (Academic Press, New York, 1963).
- (40) R. Frisch and H. Holme, Econometrica 3, 225 (1935).
- (41) M. Kalecki, Econometrica 3, 327 (1935).
- (42) N.D. Hayes, J. London Math. Soc. s1-25, 226 (1950).
- (43) M. MacDonald, Time Lags in Biological Models (Springer-Verlag, Berlin 1978).
- (44) S. Ruan, in Delay Differential Equations and Applications, edited by O. Arino, M.L. Hbid, and E.A. Dads, NATO Science Series II: Mathematics, Physics and Chemistry, Vol. 205 (Springer, Berlin, 2006) pp. 477–517.
- (45) Tao Yi and Wang Zuwang J. Theor. Biol. 187, 111 (1997).
- (46) J. Miekisz, e-print arXiv:q-bio/0703062 (2007).
- (47) U. Kuchler and E. Platen, Mathematics and Computers in Simulations 54,189 (2000).
- (48) A. Amann, E. Scholl, and W. Just, Physica A 373, 191 (2007).
- (49) S.A. Ibanez, A. Fendrik, P.I. Fierens, R.P.J. Perazzo, and D.F. Grosz, Fluctuation and Noise Letters 8, L315 (2008).
- (50) G.A. Bocharov and F.A. Rihan, J. Comput. Appl. Math. 125, 183 (2000).
- (51) T. Tian et al., J. Comput. Appl. Math. 205, 696 (2007).
- (52) C.W. Gardiner, Handbook of Stochastic Methods 2nd ed. (Springer-Verlag, New York, 1985).
- (53) M. Bambi, J. Econ. Dynamics and Control 32, 1015 (2008).
- (54) The time discretization of Eq. (2) naturally has its own effects on the stability through the numerical scheme. Choosing and will yield only small corrections to the behavior of the underlying continuous-time system.
- (55) S. Hod, Phys. Rev. Lett. 105, 208701 (2010).
- (56) In the limit of , , as can be obtained from Eq. (14).