Exact detection of direct links in networks of interacting dynamical units

The inference of an underlying network topology from local observations of a complex system composed of interacting units is usually attempted by using statistical similarity measures, such as Cross-Correlation (CC) and Mutual Information (MI). The possible existence of a direct link between different units is, however, hindered within the time-series measurements. Here we show that, for the class of systems studied, when an abrupt change in the ordered set of CC or MI values exists, it is possible to infer, without errors, the underlying network topology from the time-series measurements, even in the presence of observational noise, non-identical units, and coupling heterogeneity. We find that a necessary condition for the discontinuity to occur is that the dynamics of the coupled units is partially coherent, i.e., neither complete disorder nor globally synchronous patterns are present. We critically compare the inference methods based on CC and MI, in terms of how effective, robust, and reliable they are, and conclude that, in general, MI outperforms CC in robustness and reliability. Our findings could be relevant for the construction and interpretation of functional networks, such as those constructed from brain or climate data.


I. INTRODUCTION
Inferring the underlying connectivity topology of a complex network from observed data is nowadays the object of intense research. However, the limits for the exact inference of direct links in real-world systems composed by interacting dynamical units are still not fully understood. Understanding this limitations is often crucial in many applications in social and natural sciences. In order to infer the underlying network, usually, the observed data comes from time-series recorded at the different units. Then, a direct link between units is assumed depending on how interdependent these observations are. For example, by recording the activity of different brain regions, one wishes to infer which are the relevant structural or functional brain connections by comparing similarity patterns [14][15][16]. In general, the outcome is a complex network [12,13] that interconnects the individual units and allows for a better understanding of the overall system behavior.
Indeed, in undirected connected complex networks, every pair of units is joined by some path. Consequently, any pair of units will exchange some level of correlation or information due to the overall connecting topology. Therefore, the ability to detect a direct link between any two units is hindered within the similarity measured value. Nevertheless, when the links are homogeneous, one expects that directly connected units have larger values of the similarity measure than indirectly connected ones [25]. Then, the existence of a certain threshold, τ , that can split the similarity values into two sets can be expected. If a similarity value between two units is larger than τ , it is considered to be significant and a consequence of a direct link between the two units. Otherwise, it is less significant and a consequence of the lack of a direct link between the two units. When the strengths of the links are heterogeneous, the weak links are further hindered within the similarity measure values and a bivariate analysis can be insufficient [16]. Moreover, in the presence of strong coupling, global patterns in the system's behavior emerge, creating an effective topology which makes the underlying network inference process unfeasible. Similarly, inference fails for very weak coupling, where the system behaves as being composed of uncoupled units. Avoiding such fully coherent or incoherent behaviors is critical for the detection of the direct links.
Since different topologies are inferred for different τ values, the problem of finding the optimal τ value which recovers the largest portion of the underlying network is far from trivial. For example, in Ref. [7] the presence of dynamical noise in the individual units was shown to enable the identification of an optimal threshold giving an accurate prediction of a network topology, based solely on the measurement of dynamical correlations. However, the method requires computing the inverse matrix of the dynamical correlation matrix, which can be computationally demanding, and also the influence of non-additive noise and/or observational noise remains an open question. Other methods for link identification are based on perturbing the individual units. For example, the method proposed in Ref. [26] requires performing independent, simultaneous, and random phase resettings in all the units, which can be impractical in many real-world systems.
In this work, we show that, when the ordered values of CC (computed in absolute value, i.e., the Pearson coefficient) and/or MI (computed via ordinal pattern analysis [27][28][29][30]) exhibit a discontinuous curve, an adequate τ value allows to infer the underlying topology without errors. The exact link detection is demonstrated by considering various discrete-time dynamical units (logistic maps, circle maps, etc.) that mutually interact in different coupling topologies, including random networks. This means that both methods are able to infer the exact underlying network that interconnects the units from local time-series measurements. As a result, the structure of the interacting system is directly related to its function. We find that the existence of this τ occurs even when observational noise and heterogeneities (in the links and/or in the units) are present.
Our results are based on a critical comparison of the CC and MI inference methods effectiveness (what is the portion of the underlying network that is reconstructed correctly), robustness (how the effectiveness is affected when parameters are changed), and reliability (results yield consistent inferred networks). We conclude that MI outperforms CC as the most robust, in particular, is the least sensitive to the choice of τ value, and reliable measure. To the best of our knowledge, such reliable reconstruction of networks without errors, i.e., the most effective scenario, from time-series measurements has not been previously obtained.

II. MODEL AND METHODS
We consider logistic maps, tent maps, circle maps, and novel maps recently proposed [32] for representing nonlinear optical elements (in the following, referred to as optical maps). We let the units have a degree of heterogeneity by using non-identical parameters. For the underlying topology, we use random networks (RN) [33] and small-world networks (SW) [34] with homogeneous and heterogeneous weights. These networks are characterized by the number of nodes, N , the connectivity parameter, p, and the weights heterogeneity degree, g.

A. Network topologies
We use archetypal networks for coupling the discrete maps, namely, random [33] and small-world [34] topologies. Our Random Networks (RN) are characterized by the probability p of adding a link between every disjoint pair of nodes on a ring graph of N nodes [ Fig. 1(a)]. For p = 0 we have a ring graph and for p = 1 an all-to-all network, while random graphs are obtained for intermediate values of p [35] (where the Wigner semicircle distribution of eigenvalues is achieved and the node degrees are Poisson distributed). On the other hand, our Small-World networks (SW) are characterized by the probability p ′ of rewiring each link of a regular graph of degree k = N/4 [ Fig. 1(b)], as in Ref. [34]. These topologies define the underlying interconnections of our network of coupled maps, and our goal is to be able to infer them through the use of similarity measures. The choice of these two types of networks is due to the difference in their number of links, denoted by M . Our RN are mainly sparse networks for p ∼ 0 and N large. For finite p, the expected number of links is given by E[M p ] = p N (N − 3)/2 + N , which is the random component plus the fixed ring structure (left panel in Fig. 2). Our SW networks have a large number of links for large N , namely, E[M p ′ ] = k N/2 = N 2 /8 for every p ′ (right panel in Fig. 2). Consequently, the underlying number of connections in our system changes appreciably when using RNs or SWs topologies.
The effect of coupling heterogeneity between the maps on the similarity measure inference is dealt by using weighted networks. After an underlying topology is fixed, i.e., a particular adjacency matrix A ij is set, random weights are associated to each existing direct link. Specifically, we define a weighted network by where 1 > g ≥ 0 is the coupling heterogeneity degree parameter and ξ ij ∈ [−1, 1] is an uncorrelated zero-mean uniformly distributed random number. Symmetry in the weights is set by using W ij = W ji , which keeps the links undirected.

B. Map's equation of motion
The behavior of each map is governed by the equation where, for the logistic map, f (r, and for the optical map (namely, novel maps recently proposed [32] (2), r i is the i-th map parameter, ǫ is the coupling strength, W ij accounts for the weight of the link [Eq. (1)], and d i = N j=1 W ij is the weighted degree of node i. The behaviour of each map in the network depends on its parameter, r i , on the coupling strength, ǫ, to other maps, and on the underlying topology. For simplicity we set the map parameters to be 0 < r i = r ⋆ −ξ i δr, where r ⋆ is arbitrarily fixed to be in the chaotic regime of the map, ξ i is a random number uniformly distributed in [0, 1], and 0 ≤ δr < r ⋆ is the distribution width controlling the degree of heterogeneity in the maps. Heterogeneous coupling strengths are addressed by adding symmetric weights to the connections [Eq. (1)]. Hence, the dynamical behaviour of the units depends on the topology of the network (either RN or SW, characterized by the connectivity parameter p, number of nodes N , and weight heterogeneity degree g) and on ǫ, r ⋆ , and δr.
N time-series are obtained from the trajectories of the N maps, generated from random initial conditions. Unless otherwise stated, the length of the time-series is T = 5 × 10 4 . The CC and MI similarity measures are computed from these time-series. In particular, the MI is computed from symbolic sequences of ordinal patterns of length D = 4 [27][28][29][30].

C. Similarity Measures
The zero-lag Pearson CC (referred only by CC on the following) is defined by where σ i is the i-th map time-series standard deviation, n is the time average of node's i orbit, and T is the number of iterations that the orbit has. In particular, we use the absolute value of the CC as the similarity measure for the inference process.
The MI is defined by transforming the time-series {x (i) n } T n=0 into a symbolic sequence and then calculating the probability of appearance of each symbol in the sequence. The symbolic transformation we use is the ordinal analysis [27,28]. The ordinal analysis transforms a length D sliding window of each time-series, e.g., the vector {x The symbol is the number of permutations needed to order the components of the vector in a set of strictly increasing values. This means that each map's trajectory is encoded into a sequence of L ≃ T /D symbols if the sliding windows are non-overlapping (which is the choice we make to have equally probable symbols in the case where the time-series is random, e.g., for the surrogates of the maps trajectories).
MI is then calculated from where P (α i ) [P (β j )] is the probability of having a particular symbol α i = 1, . . . , D! [β j = 1, . . . , D!] in the encoded sequence of map i [j], namely, the frequency that α i [β j ] appears in the encoded i-th [j-th] map trajectory. Similarly, P (α i , β j ) is the joint probability that map i has a symbol α i and map j a symbol β j (possibly different than α i ) in each encoded trajectory at equal times. The choice of encoding is supported due to its simple implementation on experimental data and robustness under noisy time-series observations. Furthermore, the encoding is implemented without the need to define arbitrary partitions of the system's phase space. Also, the symbolic sequence is easily interpreted. For example, a periodic orbit of period P is transformed to a symbolic sequence with a unique symbol α if an embedding dimension D = P is used. In other words, only one of the possible D! symbols appears in the symbolic sequence of a periodic orbit of period P = D. In such a case, the symbolic entropy of the sequence is zero (H = − D! α=1 P (α) log 2 [P (α)]), hence, the MI is also zero. When the period of the orbit is different than D, the symbolic sequence has a non-null entropy, however, the joint entropy between two periodic orbits is the sum of the entropies for each orbit (independent orbits). Consequently, M I = 0 between two periodic orbits.
On the other hand, the CC value between periodic orbits depends on the phase-lag (φ) value between the two orbits. It is safe to say that for small phase-lags (φ ≪ 1) CC ∼ 1 and for large phase-lags (φ ∼ π) CC ∼ −1, hence, also close to unity in absolute value.
The threshold, τ , used to split the similarity values, is a control parameter that allows to define the inferred adjacency matrix, A τ : A τ, ij = 1 if the similarity measure between maps i and j is larger than τ , and A τ, ij = 0 otherwise. The error, ∆, between the inferred and the true adjacency matrices, is defined by The minimum value of ∆ is 0, which corresponds to an exact detection of the true underlying topology. The maximum value of ∆ is 1, and occurs only if all the links are inferred incorrectly. The effectiveness of a similarity measure (CC or MI) is quantified by the error ∆ between the true topology and the inferred topology, which is a function of the particular threshold τ chosen. We also consider the receiver operating characteristic (ROC) curve, which quantifies the true positive rate (TPR) and false positive rate (FPR), each measure being a function of τ [21,31]. We consider that a measure is effective when ∆ ≃ 0, the TPR is maximum, and the FPR is minimum. The robustness of the CC or MI is quantified in terms of how sensitive ∆ is to changes in the system's parameters (map's parameter, network topology, and heterogeneity degree) and choice of τ value. We also analyze how ∆ depends on the length of the time-series, the level of observational noise, and the size of the network. We consider that a measure is robust, when small changes to the system's parameters or optimal τ value keep ∆ ≃ 0. The reliability of a method is the ability to give consistently similar results from similar observations; hence, a measure's reliability gives an estimation of the reproducibility of the results.

III. RESULTS
In the following we present results for chaotic logistic maps (r = 4) and circle maps (r = 0.35) coupled in RNs with N = 16 and g = 0.1. Results for other maps, coupled with other network topologies and heterogeneity degrees, are found in the Supplemental Material [38].   Fig. 3(b)] measures for a particular RN (p = 0.3) of identical (δr = 0) logistic maps coupled with ǫ = 0.06. Discontinuous curves are found for both, CC and MI, though the gap for MI is larger than the one for CC. We observe that the direct connections (indicated by darker dots -black online) are found to have large similarity values, while the indirect connections (indicated by lighter dots -green online) have lower values. For comparison, the CC/MI values for ǫ = 0 are shown in light gray. We note that, in contrast to CC, the MI values are larger than the uncoupled case for all pairs of nodes, meaning that MI acknowledges the underlying topological paths between nodes in the network.
As a general result, we note that the effectiveness of a similarity measure to infer the underlying topology relies on the existence of a discontinuous jump in its ordered values. The abrupt change corresponds to a difference between the values of the similarity measure for direct connections and the values for indirect connections. We find that if a gap in the ordered sequence of CC or MI values exists, a value of τ within this gap infers the underlying topology without errors. When the gap is missing, the direct and indirect connection values are mixed within a continuous curve, hence ∆ > 0 for any τ value.
In Fig. 4(a) and (b) we indeed observe that the optimal choice of τ , namely, when ∆ = 0, is achieved for values of τ falling within the discontinuity gap of Fig. 3(a) and (b), respectively. Hence, Fig. 4 shows the effectiveness and τ -robustness of the two inference measures (in normalized units) for the same network of coupled maps as in Fig. 3  Most importantly, we found that, in general, MI is more robust than CC, in the sense that it is able to recover the underlying topology without errors for more threshold values. This is seen by comparing the lengths of the intervals where ∆ = 0 between Fig. 4(a) for CC and Fig. 4(b) for MI. The wider interval is explained by the existence of a larger gap in the values that the MI curve has in Fig. 3(b), in contrast to the smaller gap in Fig. 3(a) [however, for some parameters, CC can be more robust than MI in other aspects, as will be shown in Fig. 5(c)].
Next, we consider the receiver operating characteristic (ROC) curves, which quantify the true positive rate (TPR) and false positive rate (FPR) that each measure has as a function of τ [21,31]]. The ROC curves, shown in Fig. 4(c) and (d), provide further information about the type of errors made as a function of the threshold. When ǫ = 0, as the threshold increases from τ = 0 to τ = 1, both the TPR and FPR decrease (links are not inferred, regardless if they exist or not). On the contrary, when ǫ > 0, as τ increases, only the FPR decreases, while all the existing links are correctly inferred (the TPR remains constant). When τ is increased above the optimal range of values, then the TPR starts to decrease, as existing links are not inferred (the FPR remains constant).
In situations where the knowledge of the underlying topology is missing, the error ∆ [Eq. (5)], TPR, and FPR, cannot be computed. However, if the ordered values of the CC or MI exhibit a discontinuity gap, as in Fig. 3, the links can still be divided into two sets. The links at the left of the discontinuity correspond to indirect connections, while any value at the right reveals a direct link [25]. Every value of CC or MI inside the gap can be chosen as a possible threshold value τ and the width of the gap determines the sensitivity of the method. If the gap in the ordered CC or MI values is absent, the method is not capable to infer the correct topology.
For practical applications in real-word data, it is important to analyze how the measures depend on the length of the data set, and on how noisy the data is. In the following, observational noise is considered by adding uncorrelated zero-mean uniformly distributed noise η (i) ∈ [−1, 1] of strength Γ to each data set.  Furthermore, to make the results reliable, we average the min(∆) value that is found for each Γ and T among 5 RN realizations with equal statistical characteristics. For the sake of clarity, the white dots in the darker regions indicate where min(∆) = 0, i.e., the perfect reconstruction of all the RNs.
We can see that MI infers exactly all the network realizations for moderate noise strengths (Γ < 0.1) and orbits with T ≥ 3 × 10 4 for the logistic maps [white dots in Fig. 5(b)] and T ≥ 5 × 10 4 for the circle maps [ Fig. 5(d)]. However, we find that CC fails to provide reliable results for the logistic maps, as it only infers the underlying topology for some RN realizations [dark color in Fig. 5(a)]. On the other hand, CC outperforms MI for circle maps [white dots in Fig. 5(c)]. In general, we note that both methods are effective for moderate Γ and g when sufficient data is available.
Next we show that the exact detection of direct links is also possible when the individual units are heterogeneous (δr = 0.1), for a wide range of RN parameters and coupling strengths. This can be seen in Fig. 6, that displays the min(∆) (in color code) as a function of the connectivity (p) and the coupling strength (ǫ). Each value of min(∆) is averaged over 5 network realizations. Results for such parameter space for other maps and topologies are presented in the Supplemental Material [38]. Specifically, in Fig. 6 we see how the region where min(∆) = 0 for fixed N in the (ǫ, p) space changes depending on the units dynamical behavior. In particular, we note that with the exception of a robust window located for ǫ ∈ (0.02, 0.10) and p < 0.5 where fully inco-herent behavior is found [dark region in Fig. 6(a) and (b)], coupled chaotic logistic maps have periodic windows [light region in Fig. 6(b)] and synchronized behavior [triangular region in the upper corner of Fig. 6(a) and (b)] which make the inference impossible. On the other hand, no coherent behavior is found for the circle maps in the same (ǫ, p) space [ Fig. 6(c) and (d)]. Thus, the region where perfect inference is possible (white circles), is mainly limited by the amount of data available (disregarding ǫ ∼ 0). The ǫ-robustness of the CC or MI results, depends on the dynamic of the units composing the system and the topology (p). Although, these results are reliable and the methods are effective if dynamical coherence is avoided and sufficient data is available. In order to retrieve similar regions of perfect inference for larger (N ) networks, we find that larger data sets are needed (See Supplemental Material [38]).

IV. CONCLUSIONS
To conclude, we have shown that the Cross-Correlation coefficient (CC, calculated in absolute value) and the Mutual Information (MI, calculated from ordinal patterns) are able to infer, without errors, the underlying topology of different coupled discrete maps. In other words, CC and MI detect exactly the underlying direct links of the network of maps. We showed that, while both methods require weakly coupled units (avoiding the presence of global patterns, or strong desynchronization), the MI is in general more robust and reliable. To the best of our knowledge, reliable reconstruction of networks without errors from time-series measurements has not been previously obtained.
Our results are relevant to various fields, where complex networks of interactions are often inferred via a CC or MI statistical similarity analysis of observed timeseries. A careful consideration of the shape of the distribution of similarity values could allow for selecting optimal thresholds for the inference of direct links, as opposite to the often employed methods based on quantiles or on deviations from surrogate data.