Inferring the connectivity of coupled oscillators and anticipating their transition to synchrony through lag-time analysis

UPCommons Portal del coneixement obert de la UPC http://upcommons.upc.edu/e-prints Aquesta és una còpia de la versió author’s final draft d'un article publicat a la revista Chaos Solitons & Fractals http://hdl.handle.net/2117/182505 Article publicat / Published paper: Leyva, I.; Masoller, C. Inferring the connectivity of coupled oscillators and anticipating their transition to synchrony through lag-time analysis. Chaos solitons and fractals , 13 Abril 2020, vol. 133, 109604. DOI <10.1016/j.chaos.2020.109604> © <2020>. Aquesta versió està disponible sota la llicència CC-BYNCND 4.0 http://creativecommons.org/licenses/by-nc-nd/4.0/ Inferring the connectivity of coupled oscillators and anticipating their transition to synchrony through lag-time analysis Inmaculada Leyva 1,2and Cristina Masoller 3∗ 1Complex Systems Group & GISC, Universidad Rey Juan Carlos, Madrid, Spain 2Center for Biomedical Technology, Universidad Polit´ecnica de Madrid, Madrid, Spain 3Departament de Fisica, Universitat Politecnica de Catalunya, Barcelona, Spain Abstract The synchronization phenomenon is ubiquitous in nature. In ensembles of coupled oscillators, explosive synchronization is a particular type of transition to phase synchrony that is first-order as the coupling strength increases. Explosive sychronization has been observed in several natural systems, and recent evidence suggests that it might also occur in the brain. A natural system to study this phenomenon is the Kuramoto model that describes an ensemble of coupled phase oscillators. Here we calculate bi-variate similarity measures (the crosscorrelation, ρij , and the phase locking value, PLV ij ) between the phases, φi ( t ) and φj ( t ), of pairs of oscillators and determine the lag time between them as the time-shift, τij , which gives maximum similarity (i.e., the maximum of ρij ( τ ) or PLV ij ( τ )). We find that, as the transition to synchrony is approached, changes in the distribution of lag times provide an earlier warning of the synchronization transition (either gradual or explosive). The analysis of experimental data, recorded from Rossler-like electronic chaotic oscillators, suggests that these findings are not limited to phase oscillators, as the lag times display qualitatively similar behavior with increasing coupling strength, as in the Kuramoto oscillators. We also analyze the statistical relationship between the lag times between pairs of oscillators and the existence of a direct connection between them. We find that depending on the strength of the coupling, the lags can be informative of the network connectivity. Preprint submitted to Journal of L A T EX Templates January 23, 2020 arXiv:2001.08195v1 [nlin.AO] 22 Jan 2020 Keywords: synchronization, Kuramoto model, Rossler system, networks 1. Introduction Inferring the underlying connectivity of a complex system from the observed output signals is an important problem in nonlinear science with applications across disciplines. For example, functional brain networks are generated using statistical similarity measures applied to electroencephalogram (EEG) signals. In this approach, pairs of brain regions are considered functionally linked if there is high similarity between the time series recorded in the regions [ 1 , 2 ]. However, the ability of the functional brain networks methodology to infer the underling structure is still an open problem [ 3 ]. To demonstrate the potential of any approach as a clinical diagnostic tool it is necessary to test the inference methodology over datasets where the underlying connectivity is known. Many studies have been done in this direction, specially in gene networks [ 4 , 5 ]. However, network inference remains an open problem in the general case of dynamic units, where the interaction of the coupling structure and the dynamics of the units can hinder the detailed inference of the network connectivity (at the level single link). In parallel to the inference problem, in the analysis of functional networks is specially important to be able to detect and to predict changes in the structure and in the dynamics of the system, which can reveal the proximity to a functional transition, such as variations on the synchronization levels between brain areas. The sudden and nearly unannounced appearance of an epileptic crisis [ 6 ] is a good example, where it is a challenge to obtain information from the often subtle changes in the system dynamics [ 7 ]. In some cases the system presents a very low synchronization level just before reaching full synchronization in a sudden, irreversible way. This phenomenon, known as explosive synchronization (ES), has been observed in experimental systems [ 8 , 9 ], and it has been very recently hinted in the anesthetic-induced transition to unconsciousness [ 10 , 11 ], epilepsy [12] or chronic pain [13]. 2 In this work we analyze numerical and experimental datasets obtained from Kuramoto oscillators [ 14 ] and from R¨ossler-like electronic oscillators [ 15 ], respectively. We show how, taking into account the optimal lags in the computation of the similarity measure between each pair of time series, the resulting functional network allows to predict the proximity of the transition to synchrony, even when the early signs of the development of the process are hidden, as happens in the case of explosive synchronization. 2. Datasets 2.1. Kuramoto oscillators We consider a network of N Kuramoto phase oscillators, which is described by [14]: ˙ φi=ωi+d N X i=1 Aij sin (φj−φi) (1) where φi is the phase of the i th oscillator ( i = 1, , N ), ωi is its associated natural frequency, drawn from the homogeneous frequency distribution ω =[0,1]. The parameter d is the coupling strength, and A = {Aij} is the symmetric adjacency matrix: Aij = 1 if iand jare connected, else, Aij = 0. We use an extra parameter that allows changing the nature of the synchronization transition from a soft, gradual transition to an abrupt, first-order-like transition. As explained in Ref. [ 16 ], this can be done by explicitly imposing afrequency dissasortativity in the network, i. e., a constraint in the frequency differences between each pair of oscillators that are linked. With this procedure, we construct our network (i.e., the adjacency matrix) by randomly connecting the nodes up to obtain a preset average degree hki as in the usual Erd¨os-Reyni model, but with the condition that a pair of oscillators i and j can be connected only if |ωiωj|> γ , with γ a parameter that is refered to as frequency gap. This condition avoids that connected pairs of oscillators that have similar frequencies act as synchronization seeds. A consequence of this condition is that for high enough γthe synchronization transition is explosive [16]. 3 We use this procedure to construct networks with N =50 oscillators, with mean degree hki =5, and different values of the frequency gap γ . Over each network, the dynamics described by Eq. 1 is simulated for a wide range of coupling values d . For each value of d , the time series of the phases of the 50 oscillators are analyzed. 2.2. Electronic oscillators We also perform the network inference study over a set of experimental time series. The freely-available datasets, described in Ref. [ 15 ], come from an ensemble of N =28 networked R¨ossler chaotic electronic oscillators diffusely coupled through one of the variables. The structure of physical connections between oscillators (i.e., the adjacency matrix) is also provided, being a random matrix with 42 links (therefore 336 links do not exist). The dynamics of the variable describing the evolution of the each oscillator is given for a wide range of coupling strengths, capturing the transition from the unsynchronized behaviour to the synchronized one. As the oscillators are randomly coupled (i.e., there is no frustration parameter), the synchronization transition is not explosive. The time series have 30000 data points for every oscillator and coupling value. 3. Methods 3.1. Bi-variate similarity measures The lagged Pearson coefficient, ρij ( τ ), that is the absolute value of the crosscorrelation coefficient, is used to quantify the similarity between the time series xi ( t ) and xj ( t ) in oscillators i and j . When xi ( t ) and xj ( t ) with t = 1 . . . T are each normalized to zero-mean and unit variance, ρij(τ) can be calculated as ρij(τ) = 1 T−τ T−τ X t=1 xi(t)xj(t+τ) .(2) For each pair of oscillators i and j we calculate the lag, τij , that maximizes ρij ( τ ). We search for the maximum in the interval −τmax ≤τ≤τmax with τmax = T/ 5. 4 If several equal maximum values are found, the one with the smallest lag is selected. We measure the dynamical similarity between oscillators iand jas Sij =ρij(τij).(3) For the Kuramoto oscillators we calculate the cross-correlation between cos ( φi ( t )) and cos ( φj ( t )), while for the R¨ossler electronic oscillators, we calculate the cross-correlation between the experimental signals, i.e., the time series of the observed variables. The use of the cosine in the Kuramoto oscillators is motivated by the fact that we can consider cos ( φ ( t )) an “observed” variable (phases are usually not experimentally directly accessible, but they are calculated by using a suitable transformation, for example, the Hilbert transform). On the other hand, this measure of similarity has the drawback that oscillators with similar frequencies will have large cross-correlation values, regardless of the existence of a direct connection between them. To demonstrate that indeed the cross-correlation computed from the cosines of the phases contains useful information, we also quantify the dynamical similarity between Kuramoto oscillators i and j using the lagged phase locking value (PLV): PLVij(τ) = 1 T−τ T−τ X t=1 ei(φi(t)−φj(t−τ)).(4) As with the cross-correlation, we search for the lag that maximizes PLVij ( τ ) in the interval −τmax ≤τ≤τmax with τmax =T/5. 3.2. Global synchronization measures For the Kuramoto oscillators, the level of phase synchronization can be monitored by the order parameter given by: K=*1 N N X i=1 eiφi(t)+T (5) where h...iT denotes a time average. In the general case, as the coupling strength d increases, system represented by Eq. 1 undergoes a soft, second-order like phase 5 transition from an incoherent, desynchronized ( K∼ 0) state to a synchronous (K∼1) state, where all oscillators ultimately acquire the same frequency. For the R¨ossler chaotic electronic oscillators we use, as a global measure of synchronization, the bivariate similarity measure, Sij , averaged over all the network, hSijiij. 3.3. Network inference Using time series recorded from a set of 12 randomly coupled R¨ossler chaotic oscillators [ 17 ] it has been recently reported [ 18 ] that the mutual lags between connected pairs of oscillators are, on average, smaller than the mutual lags between unconnected pairs. Here we test if we can use this property to infer the underlying physical connectivity of the R¨ossler or of the Kuramoto oscillators, i.e., to reconstruct the network by classifying the links in two categories, existing and non-existing. To do this we first calculate, for each coupling strength, the set of N ( N− 1) / 2 similarity values, Sij = ρij or Sij = PLVij , and the corresponding set of lags, τij , that maximize ρij ( τ ) or PLVij ( τ ). Then, for each coupling strength we define two thresholds, one for the lags, τth , and one for the similarity values, Sth , and use the following criteria to classify a link as existing or not-existing: 1. SIM: the link between iand jexists if Sij > Sth, else, it does not exist. 2. AND: the link between i and j exists if τij < τth and Sij > Sth , else, it does not exist. 3. OR: the link between i and j exists if τij < τth or Sij > Sth , else, it does not exist. 3.4. Diagnostic ability of the classification criteria To quantify the efficiency of these criteria for uncovering the underlying network structure (i.e., the existing, Aij = 1, and the non-existing , Aij = 0, links) we calculate, for each criteria and for each value of the coupling, the area under the receiver operating characteristic (ROC) curve [ 19 ]. The ROC curve is obtained by varying the detection threshold and calculating the 6 1. True positives (TP): number of links that are correctly detected as existing. 2. True negatives (TN): number of links that are correctly detected as nonexisting. 3. False positives (FP): number of links that are incorrectly detected as existing. 4. False negatives (FN): number of links that are incorrectly detected as non-existing. Then, the true positive rate (TPR, also known as recall) is TP/(TP+FN)=TP/(# of existing links), and the false positive rate (FPR) is FP/(TN+FP)=FP/(# of non existing links). A ROC curve is obtained by ploting the TPR vs. FPR. The area under the ROC curve (AUC) is a measure of the goodness of a binary classifier: while random guessing gives a diagonal line, a perfect classifier has one (or more) thresholds that perfectly separate the existing and the non-existing links. In this situation, AUC=1. While the AUC is routinely used to quantify the performance of a binary classifier, in class imbalance scenarios (for example, when there are a lot of patients without a decease) it has been shown that the precision-recall curve is more informative because it does not depend on the number of true negatives [ 20 ]. The precision is the ratio of correct positive detections over all positive detections, TP/(TP+FP), and the precision-recall curve is obtained by plotting the precision vs. the recall (i.e., the TPR). Because our networks are sparce (only about 10% of the possible links exist) we also use as a performance measure the average precision that is the area under the precision-recall curve. 4. Results We begin by characterizing the synchronization transition as a function of the frequency gap, γ . Figure 1 displays the order parameter, Eq. 5, when γ = 0, 0.4 and 0.6. For γ =0 the network is constructed without restrictions, i.e., the links are distributed randomly and therefore, the synchronization transition is 7 second order. For γ =0.4 the transition to synchronization is more abrupt and it is explosive for γ=0.6. Figure 1: Kuramoto order parameter as a function of the coupling strength, for different values of the frequency gap, γ . We note that as γ increases, the synchronization transition becomes more abrupt. Figures 2 and 3 present the analysis of the Kuramoto network when the frequency gap is γ =0, using the lagged Pearson coefficient (Fig. 2) and the PLV (Fig. 3) respectively. For each procedure, panels (a) displays the correspondent average similarity value hSiji (where Sij = ρij in Fig. 2 and Sij = PLVij in Fig. 3), and panel (b) displays hτiji , averaged over all the pairs of oscillators and over time. We note that in both cases (using ρij or PLVij ), when the i−j link exists Sij tends to be larger than when it does not exist. We also note that, as the coupling strength increases, on average, for the existing links τij starts decreasing much earlier than for the non-existing links. We note that the difference is clear even if the order parameter, K , and the average similarity value, hSiji, don’t show any sign of the approaching transition. Regarding the use of the lag-time information for inferring the network structure, the panels (c) and (d) (that display the area under the ROC curve 8 [3] D. S. Bassett, O. Sporns, Network neuroscience, Nature neuroscience 20 (2017) 353. doi:https://doi.org/10.1038/nn.4502. [4] V. A. Smith, E. D. Jarvis, A. J. Hartemink, Evaluating functional network inference using simulations of complex biological systems, Bioinformatics 18 (suppl 1) (2002) S216–S224. [5] D. Marbach, J. C. Costello, R. K¨uffner, N. M. Vega, R. J. Prill, D. M. Camacho, K. R. Allison, A. Aderhold, R. Bonneau, Y. Chen, et al., Wisdom of crowds for robust gene network inference, Nature methods 9 (8) (2012) 796. [6] T. Wilkat, T. Rings, K. Lehnertz, No evidence for critical slowing down priorto human epileptic seizures, Chaos 29 (2019) 091104. doi:10.1063/1. 5122759. [7] E. van Diessen, S. J. Diederen, K. P. Braun, F. E. Jansen, C. J. Stam, Functional and structural brain networks in epilepsy: what have we learned?, Epilepsia 54 (11) (2013) 1855–1865. [8] I. Leyva, R. Sevilla-Escoboza, J. Buld´u, I. Sendina-Nadal, J. G´omezGarde˜nes, A. Arenas, Y. Moreno, S. G´omez, R. Jaimes-Re´ategui, S. Boccaletti, Phys. Rev. Lett. 108 (16) (2012) 168702. [9] S. Boccaletti, J. Almendral, S. Guan, I. Leyva, Z. Liu, I. Sendi˜na-Nadal, Z. Wang, Y. Zou, Explosive transitions in complex networks structure and dynamics: Percolation and synchronization, Physics Reports 660 (2016) 1–94. [10] M. Kim, G. A. Mashour, S.-B. Moraes, G. Vanini, V. Tarnal, E. Janke, A. G. Hudetz, U. Lee, Functional and Topological Conditions for Explosive Synchronization Develop in Human Brain Networks with the Onset of Anesthetic-Induced Unconsciousness, Frontiers in Computational Neuroscience 10 (2016) 1–15. 15 [11] M. Kim, S. Kim, G. A. Mashour, U. Lee, Relationship of topology, multiscale phase synchronization, and state transitions in human brain networks, Frontiers in computational neuroscience 11 (2017) 55. [12] Z. Wang, C. Tian, M. Dhamala, Z. Liu, A small change in neuronal network topology can induce explosive synchronization transition and activity propagation in the entire network, Scientific Reports 7 (1) (2017) 561. [13] U. Lee, M. Kim, K. Lee, C. M. Kaplan, D. J. Clauw, S. Kim, G. A. Mashour, R. E. Harris, Functional brain network mechanism of hypersensitivity in chronic pain, Scientific reports 8 (1) (2018) 243. [14] Y. Kuramoto, Chemical oscillations, waves and turbulence, Berlin: Springer, 1984. [15] R. Sevilla-Escoboza, J. M. Buldu, Synchronization of networks of chaotic oscillators: Structural and dynamical data sets, Data in Brief 7 (2016) 1185–1189. doi:10.1016/j.dib.2016.03.097. [16] I. Leyva, I. Sendina-Nadal, J. Almendral, A. Navas, M. Zanin, D. Papo, J. Buld´u, S. Boccaletti, Explosive transitions to synchronization in networked phase oscillators, Scientific Reports 3 (2013) 1281. [17] G. Tirabassi, R. Sevilla-Escoboza, J. Buld, C. Masoller, Inferring the connectivity of coupled oscillators from time-series statistical similarity analysis, Sci. Rep. 5 (2015) 10829. [18] N. Rubido, C. Masoller, Impact of lag information on network inference, Eur. Phys. J. Special Topics 227 (2018) 1243–1250. doi:10.1140/epjst/ e2018-800070-1. [19] T. Fawcett, An introduction to roc analysis, Pattern recognition letters 27 (8) (2006) 861–874. [20] D. S. Bassett, O. Sporns, The precision-recall plot is more informative than the roc plot when evaluating binary classifiers on imbalanced datasets, PLoS ONE 10 (2015) e0118432. doi:10.1371/journal.pone.0118432. 16