Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Unraveling Spurious Properties of Interaction Networks with Tailored Random Networks

  • Stephan Bialonski ,

    bialonski@gmx.net

    Affiliations Department of Epileptology, University of Bonn, Bonn, Germany, Helmholtz Institute for Radiation and Nuclear Physics, University of Bonn, Bonn, Germany, Interdisciplinary Center for Complex Systems, University of Bonn, Bonn, Germany

  • Martin Wendler,

    Affiliation Fakultät für Mathematik, Ruhr-Universität Bochum, Bochum, Germany

  • Klaus Lehnertz

    Affiliations Department of Epileptology, University of Bonn, Bonn, Germany, Helmholtz Institute for Radiation and Nuclear Physics, University of Bonn, Bonn, Germany, Interdisciplinary Center for Complex Systems, University of Bonn, Bonn, Germany

Abstract

We investigate interaction networks that we derive from multivariate time series with methods frequently employed in diverse scientific fields such as biology, quantitative finance, physics, earth and climate sciences, and the neurosciences. Mimicking experimental situations, we generate time series with finite length and varying frequency content but from independent stochastic processes. Using the correlation coefficient and the maximum cross-correlation, we estimate interdependencies between these time series. With clustering coefficient and average shortest path length, we observe unweighted interaction networks, derived via thresholding the values of interdependence, to possess non-trivial topologies as compared to Erdös-Rényi networks, which would indicate small-world characteristics. These topologies reflect the mostly unavoidable finiteness of the data, which limits the reliability of typically used estimators of signal interdependence. We propose random networks that are tailored to the way interaction networks are derived from empirical data. Through an exemplary investigation of multichannel electroencephalographic recordings of epileptic seizures – known for their complex spatial and temporal dynamics – we show that such random networks help to distinguish network properties of interdependence structures related to seizure dynamics from those spuriously induced by the applied methods of analysis.

Introduction

The last years have seen an extraordinary success of network theory and its applications in diverse disciplines, ranging from sociology, biology, earth and climate sciences, quantitative finance, to physics and the neurosciences [1][4]. There is now growing evidence that research into the dynamics of complex systems profits from a network perspective. Within this framework, complex systems are considered to be composed of interacting subsystems. This view has been adopted in a large number of modeling studies and empirical studies. It is usually assumed that the complex system under study can be described by an interaction network, whose nodes represent subsystems and whose links represent interactions between them. Interaction networks derived from empirical data (multivariate time series) have been repeatedly studied in climate science (climate networks, see [5][9] and references therein), in seismology (earthquake networks, see, e.g., [10][13]), in quantitative finance (financial networks, see e.g. [14][18] and references therein), and in the neurosciences (brain functional networks, see [19], [20] for an overview). Many interaction networks have been reported to possess non-trivial properties such as small-world architectures, community structures, or hubs (nodes with high centrality), all of which have been considered to be characteristics of the dynamics of the complex system.

When analyzing empirical data one is faced with the challenge of defining nodes and inferring links from multivariate noisy time series with only a limited number of data points due to stationarity requirements. Different approaches varying to some degree across disciplines have been proposed. For most approaches, each single time series is associated with a node and inference of links is based on time series analysis techniques. Bivariate time series analysis methods, such as the correlation coefficient, are used as estimators of signal interdependence which is assumed to be indicative of an interaction between different subsystems. Inferring links from estimates of signal interdependence can be achieved in different ways. Weighted interaction networks can be derived by considering estimated values of signal interdependence (sometimes mapped via some function) as link weights. Since methods characterizing unweighted networks are well-established and readily available, such networks are more frequently derived from empirical data. Besides approaches based on constructing minimum spanning trees (see, e.g., reference [14]), on significance testing [21][23], or on rank-ordered network growth (see, e.g., reference [15]), a common practice pursued in many disciplines is to choose a threshold above which an estimated value of signal interdependence is converted into a link (“thresholding”, see, e.g., references [5], [12], [16], [20]). Following this approach, the resulting unweighted interaction networks have been repeatedly investigated employing various networks characteristics, among which we mention the widely-used clustering coefficient and average shortest path length to assess a potential small-world characteristic, and the node degrees in order to identify hubs.

As studies employing the network approach grow in numbers, the question arises as to how informative reported results are with respect to the investigated dynamical systems. To address this issue, properties of interaction networks are typically compared to those obtained from network null models. Most frequently, Erdös-Rényi random networks [24] or random networks with a predefined degree distribution [25], [26] serve as null models; network properties that deviate from those obtained from the null model are considered to be characteristic of the investigated dynamical system. Only in a few recent studies, results obtained from network analyses have been questioned in relation to various assumptions underlying the network analysis approach. Problems pointed out include: incomplete data sets and observational errors in animal social network studies [27]; representation issues and questionable use of statistics in biological networks (see [28] and references therein); challenging node and link identification in the neurosciences [29][31]; the issue of spatial sampling of complex systems [31][33]. This calls not only for a careful interpretation of results but also for the development of appropriate null models that incorporate knowledge about the way networks are derived from empirical data.

We study – from the perspective of field data analysis – a fundamental assumption underlying the network approach, namely that the multivariate time series are obtained from interacting dynamical processes and are thus well represented by a model of mutual relationships (i.e., an interaction network). Visual inspection of empirical time series typically reveals a perplexing variety of characteristics ranging from fluctuations on different time scales to quasi-periodicity suggestive of different types of dynamics. Moreover, empirical time series are inevitably noisy and finite leading to a limited reliability of estimators of signal interdependencies. This is aggravated with the advent of time-resolved network analyses where a good temporal resolution often comes at the cost of diminished statistics. Taken together, it is not surprising that the suitability of the network approach is notoriously difficult to judge prior to analysis.

We here employ the above-mentioned thresholding-approach to construct interaction networks for which we estimate signal interdependence with the frequently used correlation coefficient and the maximum cross correlation. We derive these networks, however, from multivariate time series of finite length that are generated by independent (non-interacting) processes which would a priori not advocate the notion of a representation by a model of mutual relationships. In simulation studies we investigate often used network properties (clustering coefficient, average shortest path length, number of connected components). We observe that network properties can deviate pronouncedly from those observed in Erdös-Rényi networks depending on the length and the spectral content of the multivariate time series. We address the question whether similar dependencies can also be observed in empirical data by investigating multichannel electroencephalographic (EEG) recordings of epileptic seizures that are known for their complex spatial and temporal dynamics. Finally, we propose random networks that are tailored to the way interaction networks are derived from multivariate empirical time series.

Methods

Interaction networks are typically derived from multivariate time series () in two steps. First, by employing some bivariate time series analysis method, interdependence between two time series and () is estimated as an indicator for the strength of interaction between the underlying systems. A multitude of estimators [34][40], which differ in concepts, robustness (e.g., against noise contaminations), and statistical efficiency (i.e., the amount of data required), is available. Studies that aim at deriving interaction networks from field data frequently employ the absolute value of the linear correlation coefficient to estimate interdependence between two time series. The entries of the correlation matrix then read(1)where and denote mean value and the estimated standard deviation of time series . Another well established method to characterize interdependencies is the cross correlation function. Here we use the maximum value of the absolute cross correlation between two time series,(2)with(3)to define the entries of the cross correlation matrix . As practiced in field data analysis, we normalize the time series to zero mean before pursuing subsequent steps of analysis. Note that is then the maximum value of the absolute cross covariance function. Both interdependence estimators are symmetric ( and ) and are confined to the interval . High values indicate strongly interdependent time series while dissimilar time series result in values close to zero for sufficiently large.

Second, the adjacency matrix representing an unweighted undirected interaction network is derived from (or ) by thresholding. For a threshold entries and of are set to (representing an undirected link between nodes and ) for all entries (, respectively) with , and to zero (no link) otherwise. In many studies is not chosen directly but determined such that the derived network possesses a previously specified mean degree , where denotes the degree of , i.e., the number of links connected to node . More frequently, is chosen such that the network possesses a previously specified link density . We will follow the latter approach and derive networks for fixed values of .

To characterize a network as defined by , a plethora of methods have been developed. Among them, the clustering coefficient and the average shortest path length are frequently used in field studies. The local clustering coefficient is defined as(4) represents the fraction of the number of existing links between neighbors of node among all possible links between these neighbors [1], [2], [41]. The clustering coefficient of the network is defined as the mean of the local clustering coefficients,(5) quantifies the local interconnectedness of the network and .

The average shortest path length is defined as the average shortest distance between any two nodes,(6)and characterizes the overall connectedness of the network. denotes the length of the shortest path between nodes and . The definition of the average shortest path length varies across the literature. Like some authors, we here include the distance from each node to itself in the average (). Exclusion will, however, just change the value by a constant factor of .

If a network disintegrates into a number of different connected components, there will be pairs of nodes , for which no connecting path exists, in which case one usually sets and thus . In order to avoid this situation, in some studies in equation (6) is replaced by . The quantity defined this way is called efficiency [42], [43]. Another approach, which we will follow here and which is frequently used in field studies, is to exclude infinite values of from the average. The average shortest path length then reads(7)where(8)denotes the set of all pairs of nodes with finite shortest path. The number of such pairs is given by . Note that for .

In field studies, values of and obtained for interaction networks are typically compared with average values obtained from an ensemble of random Erdös-Rényi (ER) networks [24]. Between every pair of nodes is a link with probability , and links for different pairs exist independently from each other. The expectation value of the clustering coefficient of ER networks is [2]. The dependence of the average shortest path length of ER networks on and is more complicated (see references [2], [44]). Almost all ER networks are connected, if . ER networks with a predefined number of links (and thus link density) can also be generated by successively adding links between randomly chosen pairs of nodes until the predefined number of links is reached. During this process, multiple links between nodes are avoided.

Results

Simulation studies

We consider time series , , whose entries are drawn independently from the uniform probability distribution on the interval . We will later study the impact of different lengths of these random time series on network properties. In order to enable us to study the effects of different spectral contents on network properties, we add the possibility to low-pass filter by considering(9)where , and . By definition . With the size of the moving average we control the spectral contents of time series. We here chose this ansatz for the sake of simplicity, for its mathematical treatability, and because the random time series with different spectral contents produced this way show all properties we want to illustrate.

In the following we will study the influence of the length of time series on network properties by considering for different . For a chosen value of we determine realizations of and we denote each realization with . When studying the influence of the spectral content we will consider with different and with . We chose this value of because we are interested in investigating time series of short length as typically considered in field studies. For a chosen value of we determine realizations of and we denote realization with .

In order to keep the line of reasoning short and clear, we will present supporting and more rigorous mathematical results in Appendix S1 and refer to them in places where needed. In addition, since we observed most simulation studies based on to yield qualitatively the same results as those based on , we will present results based on only and report results of our studies based on whenever we observed qualitative differences.

Clustering coefficient.

Let denote the absolute value of the empirical correlation coefficient estimated for time series and , and let us consider realizations, . Because of the independence of processes generating the time series and because of the symmetry of the correlation coefficient, we expect the values of the empirical correlation coefficient calculated for finite and fixed to be distributed around the mean value . The variance of this distribution will be higher the lower we choose . If we sample one value out of the values it is almost surely that . Thus there are thresholds for which we would establish a link between the corresponding nodes and when deriving a network. Let us now consider a network of nodes whose links are derived from time series as before. For some the network will possess links and . We expect to observe for some fixed to be higher the larger the variance of the distribution of . Likewise, for fixed values of we expect to find to be higher the lower we choose a value of .

As a first check of this intuition we derive an approximation for the edge density by making use of the asymptotic limit (, see Appendix S1, Lemma 2 for details),(10)where denotes the cumulative distribution function of a standard normal distribution. In figure 1 (top left) we show the dependence of on for exemplary values of . Indeed, is decreasing in and .

thumbnail
Figure 1. Simulation results for the edge density, the clustering coefficient, and the effective length.

Top row: Dependence of edge density (left) and of clustering coefficient (right) on the threshold for different values of the size of the moving average and of the length of time series. Values of edge density obtained with the asymptotic limit (equation (10)) are shown as lines (top left). Bottom left: Dependence of the ratio on edge density . Note, that we omitted values of estimated quantities obtained for since the accuracy of the statistics is no longer guaranteed. Bottom right: Dependence of effective length as determined by equation (14) (black line) and its numerical estimate (red markers) on .

https://doi.org/10.1371/journal.pone.0022826.g001

The concession of taking the asymptotic limit when deriving equation (10) may limit its validity in the case of small values of in which we are especially interested. Thus, we approach this case by simulation studies. Let us consider values of obtained for realizations of two time series , . We estimate the edge density by(11)where for , and else. Note that does not depend on . This is because represents the (numerically determined) probability that there is a link between two vertices. The dependence of on for different values of shown in figure 1 (top left) indicates a good agreement between and for larger values of but an increasing difference for .

We proceed by estimating the clustering coefficient for our model using realizations of three time series , by(12)The dependence of on for various is shown in the top right part of figure 1. For fixed , decreases from with increasing values of which one might expect due to the decrease of . However, we also observe for that takes on higher values the lower .

In order to investigate whether the clustering coefficients of our networks differ from those of Erdös-Rényi networks we use equation (11) and obtain with . This allows the comparison with by considering the ratio . Remarkably, for a range of values of and (cf. lower left part of figure 1). even increases for small and . This indicates that there is a relevant dependence between the three random variables , , and for different indices and small . For and fixed edge density, converges to because the dependence between the random variables , , vanishes (i.e., the random variables will converge in distribution to independent normal random variables).

In order to gain deeper insights into the influence of the spectral contents of random time series on the clustering coefficient, we repeat the steps of analysis with time series , where is kept fix, and we choose different values of . Figure 1 (top panels and lower left) shows that the higher the amount of low-frequency contributions (large ) the higher and (for ), and the higher (for ). The difference between Erdös-Rényi networks and our time series networks becomes more pronounced () the smaller and the higher .

Given the similar dependence of , , and on and , we hypothesize that the similarity can be traced back to similar variances of and for some values of and . By making use of the asymptotic variance of the limit distributions of , we derive an expression relating Var and Var to each other (see Appendix S1, Lemma 1),(13)We are now able to define an effective length of time series,(14)for which . In the lower right part of figure 1 we show in dependence on . To investigate whether equation (14) also holds for small values of , we determine numerically, for different values of , for as well as for some chosen values of . Eventually, we determine for each value of a value of , for which and curves match in a least-squares sense, and denote this value with (see the lower right part of figure 1). We observe a maximum deviation and conclude that equation (14) indeed holds for small length of time series. Moreover, numerically determined dependencies of on , on , as well as on for pairs of values show a remarkable similarity to those dependencies obtained for pairs of values .

Thus, the clustering coefficient of networks derived from random time series of finite length and/or with a large amount of low-frequency contributions is higher than the one of Erdös-Rényi (ER) networks – independently of the network size (cf. equation (12)). This difference becomes more pronounced the lower the edge density , the lower the length of time series, and the larger the amount of low-frequency contributions. These results point us to an important difference between ER networks and our model networks: possible edges in ER networks are not only (1) equally likely but also (2) independently chosen to become edges. While property (1) is fulfilled in our model networks, property (2) is not.

Average shortest path length.

Next we study the impact of the length of time series and of the amount of low-frequency contributions on the average shortest path length of our model networks by employing a similar but different simulation approach as used in the previous section. To estimate , we consider networks with a fixed number of nodes (). We derive our model networks by thresholding , , . Let denote the average shortest path length for network with and different values of , and the average shortest path length for network with fixed value of () and different values of . With we refer to the average shortest path length obtained for the -th ER network of size and edge density . Mean values over realizations will be denoted as , , and respectively. Finally, we define and .

In figure 2 we show the dependence of and on for various values of and . All quantities decrease as increases which can be expected due to additional edges reducing the average distances between pairs of nodes of the networks. For fixed , takes on higher values the higher or the lower . With equation (14) we have which resembles the results obtained for the clustering coefficient. Differences between the average shortest path lengths of our model networks and ER networks (as characterized by ) become more pronounced the higher and the lower . For edge densities typically reported in field studies (), however, differences are less pronounced (, cf. figure 2 right) than the ones observed for the clustering coefficient ( for selected values of and , cf. figure 1 bottom left). We obtained qualitatively similar results for small () and large numbers of nodes ().

thumbnail
Figure 2. Simulation results for the average shortest path length.

Dependence of the average shortest path length (left) and of the ratio (right) on edge density for different values of the size of the moving average and of the length of time series. Lines are for eye-guidance only.

https://doi.org/10.1371/journal.pone.0022826.g002

Number of connected components and degree distribution.

Since the number of connected components of a given network might affect network characteristics such as the average shortest path length (see equation (7)), we investigate the impact of different length of time series and of the amount of low-frequency contributions on the average number of connected components of the networks derived from and . We determine as the mean of realizations of the corresponding networks and with we denote the mean value of the number of connected components in realization of ER networks. For the edge densities considered here we observe ER networks to be connected (cf. figure 3), , which is in agreement with the connectivity condition for ER networks, (for ). Similarly, we observe , even for small values of (cf. figure 3 right). In contrast, takes on higher values the lower and the higher (cf. figure 3 left). In order to achieve a better understanding of these findings, we determine degree probability distributions of our model networks. Let denote the estimated probability of a node to possess a degree , i.e., . With we will denote the estimated degree distribution for networks derived from . We briefly recall that the degree distribution of ER networks follows a binomial distribution,(15)which we show in figure 4 for and various (top panels and lower left panel). In the same figure we present our findings for for various values of and . We observe to be equal to within the error to be expected due to the limited sample size used for the estimation. For , however, we observe striking differences in comparison to the previous degree distributions. In particular, for decreasing and higher , the probability of nodes with degree increases, which leads to networks with disconnected single nodes, thereby increasing the number of connected components of the network.

thumbnail
Figure 3. Simulation results for the number of connected components.

Dependence of the number of connected components on the edge density for different values of the size of the moving average (left, for ) and of the length of time series (right, for ). Lines are for eye-guidance only.

https://doi.org/10.1371/journal.pone.0022826.g003

thumbnail
Figure 4. Simulation results for the degree distribution.

(a–c) Degree distributions estimated for realizations of networks derived from time series () via thresholding using various edge densities and for selected values of the size of the moving average and of the length of time series. The symbol legend in (a) also holds for (b) and (c). (d) Dependence of correlation () between node degrees and power content in the lower frequency range on the size of the moving average. Mean values of correlations obtained for realizations of networks for each value of are shown as crosses and standard deviations as error bars. Stars indicate significant differences in comparison to (Bonferroni corrected pair-wise Wilcoxon rank sum tests for equal medians, ). Lines are for eye-guidance only.

https://doi.org/10.1371/journal.pone.0022826.g004

We hypothesize that the observed differences in the number of connected components as well as in the degree distributions are related to differences in the spectral content of different realizations of time series for . In particular, a node with a low degree might be associated with a time series , which possesses, by chance, a small relative amount of low frequency contributions (or, equivalently, a large relative amount of high frequency contributions).

In order to test this hypothesis, we generate realizations of time series and estimate their periodograms for frequencies using a discrete Fourier transform [45]. denotes the Nyquist frequency, and periodograms are normalized such that . From the same time series, we then derive the networks using and determine the degrees . For some fixed we define the total power above (upper frequency range) as , and the total power below (lower frequency range) as . For each realization we estimate the correlation coefficients between the degrees and the corresponding total power contents in upper and lower frequency range, and , respectively, and determine their mean values, and . Note that by construction. We choose such that 40% of the total power of the filter function associated with the moving average is contained within the frequency range .

For increasing we observe in the lower right panel of figure 4 the degrees to be increasingly correlated with , which corresponds to an anti-correlation of degrees with . Thus, as hypothesized above, the observed differences in the degree distributions can indeed be related to the differences in the power content of the time series. We mention that the exact choice of does not sensitively affect the observed qualitative relationships as long as is fulfilled.

We briefly summarize the results obtained so far, which indicate a striking difference between networks derived from independent random time series using or (cf. equations (1) and (2)) and corresponding ER networks. First, we observed the clustering coefficient and the average shortest path length of our networks to be higher the lower the length of the time series (cf. figures 1 and 2). Second, for some fixed we observed and to be higher the larger the amount of low frequency components (as parametrized by ) in the time series. In addition, these contributions led to an increasing number of connected components in our networks and to degree distributions that differed strongly from those of the corresponding ER networks (cf. figures 3 and 4). We mention that as defined here (cf. equation (7)) tends to decrease for networks with an increasing number of connected components, and for . thus depends non-trivially on the amount of low frequency components in the time series. Third, for small edge densities and for short time series lengths or for a large amount of low frequency components, the clustering coefficient deviates more strongly from the one of corresponding ER networks () than the average shortest path length (; cf. figure 2 right and figure 1 (bottom left)). Networks derived from independent random time series can thus be classified as small world networks if one uses and as practical criterion, which is often employed in various field studies (cf. [31] and references therein).

Field data analysis

The findings obtained in the previous section indicate that strong low frequency contributions affect network properties and in a non-trivial way. We now investigate this influence in electroencephalographic (EEG) recordings of epileptic seizures that are known for their complex spatial and temporal changes in frequency content [46][49]. We analyze the multichannel ( channels) EEGs from 60 patients capturing 100 epileptic seizures reported in reference [50]. All patients had signed informed consent that their clinical data might be used and published for research purposes. The study protocol had previously been approved by the ethics committee of the University of Bonn. During the presurgical evaluation of drug-resistant epilepsy, EEG data were recorded with chronically implanted strip, grid, or depth electrodes from the cortex and from within relevant structures of the brain. The data were sampled at 200 Hz within the frequency band Hz using a 16-bit analog-to-digital converter. Electroencephalographic seizure onsets and seizure ends were automatically detected [51], and EEGs were split into consecutive non-overlapping windows of 2.5 s duration ( sampling points). Time series of each window were normalized to zero mean and unit variance. We determined and for all combinations of EEG time series from each window and derived networks with a fixed edge density in order to exclude possible edge density effects. With and as well as and we denote characteristics of networks based on and , respectively. In order to simplify matters, we omit the window indexing in the following.

We investigate a possible influence of the power content of EEG time series on the clustering coefficient and the average shortest path length by comparing their values to those obtained from ensembles of random networks that are based on properties of the EEG time series at two different levels of detail. For the first ensemble and for each patient we derive networks from random time series with a power content that approximately equals the mean power content of all EEG time series within a window. Let denote the estimated periodogram of each EEG time series , and with we denote the mean power for each frequency component over all time series. We normalize such that . We generate random time series of length whose entries are independently drawn from a uniform probability distribution, and we filter these time series in the Fourier domain using as filter function. We normalize the filtered time series to zero mean and unit variance and derive a network based on or (). We use 20 realizations of such networks per window in order to determine the mean values of network characteristics and as well as and based on or , respectively. Since the power spectra of all time series equal each other, these random networks resemble the ones investigated in the previous section.

With the second ensemble, we take into account that the power content of EEG time series recorded from different brain regions may differ substantially. For this purpose we make use of a well established method for generating univariate time series surrogates [52], [53], which have power spectral contents and amplitude distributions that are practically indistinguishable from those of EEG time series but are otherwise random. Amplitudes are iteratively permuted while the power spectrum of each EEG time series is approximately preserved. Since this randomization scheme destroys any significant linear or non-linear dependencies between time series, it has been successfully applied to test the null hypothesis of independent linear stochastic processes. For each patient, we generated 20 surrogate time series for each EEG time series from each recording site and each window, and derived networks based on either or (). Mean values of characteristics of these random networks are denoted as and as well as and , respectively.

We begin with an exemplary recording of a seizure of which we show in figure 5 (left) the temporal evolution of the relative amount of power in the - (0–4 Hz, ), - (4–8 Hz, ), - (8–12 Hz, ), and - (12–20 Hz, ) frequency bands. Prior to the seizure the -band contains more than 50% of the total power which is then shifted towards higher frequencies and back towards low frequencies at seizure end. is even higher after the seizure than prior to the seizure.

thumbnail
Figure 5. Evolving relative amount of power during epileptic seizures.

(Left) Relative amount of power contained in the - (, black), - (, blue), - (, green), and - (, red) frequency bands during an exemplary seizure. Profiles are smoothed using a four-point moving average. Grey-shaded area marks the seizure. (Right) Mean values (, , , ) of the relative amount of power averaged separately for pre-seizure, discretized seizure, and post-seizure time periods of 100 epileptic seizures. Lines are for eye-guidance only.

https://doi.org/10.1371/journal.pone.0022826.g005

In figure 6 we show the temporal evolution of network properties obtained for this recording based on (top panels) and (bottom panels). During the seizure both the clustering coefficients and and the average shortest path lengths and show pronounced differences to the respective properties obtained from the random networks. These differences are less pronounced prior to and after the seizure, where and even approach the values of and , respectively. and decrease during the seizure and already increase prior to seizure end, resembling the changes of (cf. figure 5 (left)). This is in accordance with results of our simulation studies: there we observed the clustering coefficient to be higher the larger the amount of low frequency components in the time series; this could also be observed, but to a much lesser extent, for the average shortest path length. Indeed, and vary little over time, and is only slightly increased after the seizure, reflecting the high amount of power in the -band.

thumbnail
Figure 6. Evolving network properties during an exemplary epileptic seizure.

Network properties and (top row, black lines) as well as and (bottom row, black lines) during an exemplary seizure (cf. figure 5 (left)). Mean values and standard deviations of network properties obtained from surrogate time series (, , , ) are shown as blue lines and blue shaded areas, respectively, and mean values and standard deviations of network properties obtained from the overall power content model (, , , ) are shown as red lines and red shaded areas, respectively. Profiles are smoothed using a four-point moving average. Grey-shaded area marks the seizure. For corresponding Erdös-Rényi networks and for all time windows.

https://doi.org/10.1371/journal.pone.0022826.g006

We only observe small deviations between and as well as between and , which appear to be systematic (for many windows and ). These suggest that for interaction networks derived from , both random network ensembles appear appropriate to characterize the influence of power in low frequency bands on clustering coefficient and the average shortest path length. In contrast, we observed differences between and , as well as between and . These differences were most pronounced during the seizure and for and also after the seizure. This finding indicates that the clustering coefficient and average shortest path length of interaction networks derived from depend sensitively on the power contents of EEG time series recorded from different brain regions. Thus, for these interaction networks only the random networks that account for the complex changes in frequency content of different brain regions prior to, during, and after seizures appear appropriate to characterize the influence of power in low frequency bands on clustering coefficient and the average shortest path length.

We continue by studying properties of networks derived from the EEG recordings of all 100 focal onset seizures. Due to the different durations of seizures (mean seizure duration: s) we partitioned each seizure into 10 equidistant time bins (see reference [50] for details) and assigned the time-dependent network properties to the respective time bins. For each seizure we included the same number of pre-seizure and post-seizure windows in our analysis and assigned the corresponding time-dependent network properties to one pre-seizure and one post-seizure time bin. Within each time bin we determined the mean value (e.g., ) and its standard error for each property. In figure 5 (right), we show for each time bin the mean values of the relative amount of power in different frequency bands of all seizure recordings (, , , ). Similar to the exemplary recording (cf. figure 5 (left)), we observe a shift in the relative amount of power in low frequencies prior to seizures towards higher frequencies during seizures and back to low frequencies at seizure end. The amount of power in the -band is on average higher in the post-seizure bin than in the pre-seizure bin.

In figure 7 we show the mean values of properties of networks in each time bin for all seizures. We observe , , , , , and to decrease during seizures and to increase prior to seizure end thereby roughly reflecting the amount of power contained in low frequencies (cf. figure 5 (right), ). and and to a lesser extent also and roughly follow the same course in time, however, with a slight shift in the range of values as already observed in the exemplary recording of a seizure (cf. figure 6). Differences between both random network ensembles are most pronounced in network properties based on , i.e., between and as well as between and . This corroborates the observation that the clustering coefficient and the average shortest path length of the random networks based on depend more sensitively on the power contents of EEG time series recorded from different brain regions than the respective quantities derived from . While and show a similar course in time, reaching a maximum in the middle of the seizures, we observe a remarkable difference between and prior to end of the seizures, where the amount of power in low frequencies is large. While decreases at the end of the seizures, does not and remains elevated after seizures. Interestingly, considering the corresponding quantities obtained from the second random network ensemble, fluctuates around and does not increase at the end of seizures, while, in contrast, increases at the end of the seizures, traversing an interval of values roughly three times larger than the interval containing values of . Taken together these findings suggest that the pronounced changes of the frequency content of EEG time series seen during epileptic seizures influence the values of the clustering coefficient and the average shortest path length.

thumbnail
Figure 7. Evolving network properties averaged over 100 epileptic seizures.

Mean values (black) of network properties (top left), (top right), (bottom left), and (bottom right) averaged separately for pre-seizure, discretized seizure, and post-seizure time periods of 100 epileptic seizures. Mean values of corresponding network properties obtained from the first and the second ensemble of random networks are shown as red and blue lines, respectively. All error bars indicate standard error of the mean. Lines are for eye-guidance only.

https://doi.org/10.1371/journal.pone.0022826.g007

A comparison of some value of a network property with the one obtained for a random network with the same edge density and number of nodes is typically achieved by calculating their ratio. If ER networks are used for comparison, the value of a network property is rescaled by a constant factor. In this case, the time-dependent changes of network properties shown in figure 7 will be shifted along the ordinate only. In order to take into account the varying power content of EEG time series recorded from different brain regions we instead normalize the clustering coefficients and the average shortest path lengths with the corresponding quantities from the second random network ensemble , , , and (cf. figure 8). We observe the normalized network properties to describe a concave-like movement over time indicating a reconfiguration of networks from more random (before seizures) towards a more regular (during seizures) and back towards more random network topologies. This is in agreement with previous observations using a different and seldom used thresholding method [50].

thumbnail
Figure 8. Evolving normalized network properties averaged over 100 epileptic seizures.

Mean values of and (left) as well as and (right) averaged separately for pre-seizure, discretized seizure, and post-seizure time periods of 100 epileptic seizures. All error bars indicate standard error of the mean. Lines are for eye-guidance only.

https://doi.org/10.1371/journal.pone.0022826.g008

Discussion

The network approach towards the analysis of empirical multivariate time series is based on the assumption that the data is well represented by a model of mutual relationships (i.e., a network). We studied interaction networks derived from finite time series generated by independent processes that would not advocate a representation by a model of mutual relationships. We observed the derived interaction networks to show non-trivial network topologies. These are induced by the finiteness of data, which limits reliability of estimators of signal interdependence, together with the use of a frequently employed thresholding technique. Since the analysis methodology alone can already introduce non-trivial structure in the derived networks, the question arises as to how informative network analysis results obtained from finite empirical data are with respect to the studied dynamics. This question may be addressed by defining and making use of appropriate null models. In the following, we briefly discuss two null models that are frequently employed in field studies.

Erdös-Rényi (ER) networks represent one of the earliest and best studied network models in mathematical literature and can be easily generated. They can be used to test whether the network under consideration complies with the notion of a random network in which possible edges are equally likely and independently chosen to become edges. We observed that clustering coefficient and average shortest path length for interaction networks derived from finite random time series differed pronouncedly from those obtained from corresponding ER networks, which would likely lead to a classification of interaction networks as small-world networks. Since the influence of the analysis methodology is not taken into account with ER networks, they may not be well suited for serving as null models in studies of interaction networks derived from finite time series.

Another null model is based on randomizing the network topology while preserving the degrees of nodes [26], [54], [55]. It is used to evaluate whether the network under consideration is random under the constraint of a given degree sequence. Results of our simulation studies point out that the structures induced in the network topology by the way how networks are derived from empirical time series cannot be related to the degree sequence only. We observed that and from interaction networks remarkably depended on the finiteness of the data, while the degree distribution did not (cf. figure 4 (a–c), ). The usefulness of degree-preserving randomized networks has also been subject of debate since they do not take into account different characteristics of the data and its acquisition [56], [57]. Moreover, the link-switching algorithm frequently employed for generating such networks has been shown to non-uniformly sample the space of networks with predefined degree sequence (see, e.g., references [25], [58]). This deficiency can be addressed by using alternative randomization schemes (see, e.g., [58][60] and references therein).

We propose to take into account the finite length and the frequency contents of time series when defining null models. For this purpose we applied the same methodological steps as in field data analysis (estimation of signal interdependence and thresholding of interdependence values to define links) but used surrogate time series [53] to derive random networks (second ensemble). These surrogate time series comply with the null hypothesis of independent linear stochastic processes and preserve length, frequency content, and amplitude distribution of the original time series. For these random networks, we observed (in our simulation studies) dependencies between properties of networks and properties of time series: the clustering coefficient , and, to a lesser extent, the average shortest path length are higher the higher the relative amount of low frequency components, the shorter the length of time series, and the smaller the edge density of the network. Results obtained from an analysis of interaction networks derived from multichannel EEG recordings of one hundred epileptic seizures confirm that the pronounced changes of the frequency content seen during seizures influence the values of and . Comparing these network characteristics with those obtained from our random networks allowed us to distinguish aspects of global network dynamics during seizures from those spuriously induced by the applied methods of analysis.

Our random networks will likely be classified as small-world networks when compared to ER networks which might indicate that small-world topologies in networks derived from empirical data as reported in an ever increasing number of studies can partly or solely be related to the finite length and frequency content of time series. If so, small-world topologies would be an overly complicated description of the simple finding of finite time series with a large amount of low frequency components. In this context, our approach could be of particular interest for studies that deal with short time series and low frequency contents, as, for example, is the case in resting state functional magnetic resonance imaging studies (see, e.g., references [61][65]). In such studies, taking into account potential frequency effects could help to unravel information on the network level that would be otherwise masked.

We observed the degrees of nodes of our random networks to be correlated with the relative amount of power in low-frequencies in the respective time series (cf. figure 4). The degree of a node has been used in field studies as an indicator of its centrality in the network (see, e.g., [2], [66] and references therein). Particular interest has been devoted to nodes which are highly central (hubs). In this context it would be interesting to study whether findings of hubs in interaction networks can partly or solely be explained by the various frequency contents of time series entering the analysis. In such a case, hubs would be a complicated representation of features already present on a single time series level. We are confident that our random networks can help to clarify this issue.

Our simulation studies were based on the simplified assumption that power spectra of all time series from which a network is derived are approximately equal. The dependencies of and on the power content could also be observed qualitatively for networks derived from EEG time series – that were recorded from different brain regions and whose power spectra may differ substantially among each other – but only if link definition was based on thresholding the values of the correlation coefficient (). Thus, estimating mean power spectra of multivariate time series can provide the experimentalist with a rule of thumb for the potential relative increase of and in different networks based on the correlation coefficient. This rule of thumb, however, might not be helpful if the maximum value of the absolute cross correlation () is used to estimate signal interdependencies. In this case, and depended sensitively on the heterogeneity of power spectra (see the second random network ensemble). It would be interesting to investigate in future studies, which particular properties of and can be accounted for these differences.

We close the discussion with two remarks, the first being of interest for experimentalists. Our findings also shed light on a network construction technique that relies on significance testing in order to decide upon defining a link or not [21]. For this purpose, a null distribution of a chosen estimator of signal interdependence () is generated for each pair of time series and a link is established if the null hypothesis of independent processes generating the time series can be rejected at a predefined significance level. It was suggested in Ref. [21] to use a limited subset of time series in order to minimize computational burden when generating null distributions. Our findings indicate that networks constructed this way will yield an artificially increased number of false positive or of false negative links which will depend on the frequency contents of time series being part or not part of the subset. Our second remark is related to network modeling. By choosing some threshold and generating time series that satisfy the relation between the size of the moving average and the length of time series, networks can be generated which differ in their degree distributions but approximately equal in their clustering coefficient and average shortest path length. This property could be of value for future modeling studies.

To summarize, we have demonstrated that interaction networks, derived from finite time series via thresholding an estimate of signal interdependence, can exhibit non-trivial properties that solely reflect the mostly unavoidable finiteness of empirical data, which limits the reliability of signal interdependence estimators. Addressing these influences, we proposed random network models that take into account the way interaction networks are derived from the data. With an exemplary time-resolved analysis of the clustering coefficient and the average shortest path length of interaction networks derived from multichannel electroencephalographic recordings of one hundred epileptic seizures, we demonstrated that our random networks allow one to gain deeper insights into the global network dynamics during seizures. Here we concentrated on and but we also expect other network characteristics to be influenced by the methodologies used to derive interaction networks from empirical data. Analytical investigations of properties of our random networks and the development of formal tests for deviations from these networks may be regarded as promising topics for further studies. Other research directions are related to the framework we proposed to generate random networks from time series. For example, parts of the framework may be exchanged in order to study network construction methodologies other than thresholding (e.g., based on minimum spanning trees [14] or based on allowing weighted links) or other widely used linear and nonlinear methods for estimating signal interdependence [34], [35], [39]. Other surrogate concepts [67][72] may allow for defining different random networks tailored to various purposes. We believe that research into network inference from time series and into random network models that incorporate knowledge about the way networks are derived from empirical data can decisively advance applied network science. This line of research can contribute to gain a better understanding of complex dynamical systems studied in various scientific fields.

Supporting Information

Acknowledgments

We thank Marie-Therese Kuhnert, Gerrit Ansmann, and Alexander Rothkegel for their helpful comments and Paula Daniliuc for proofreading the manuscript.

Author Contributions

Conceived and designed the experiments: SB MW KL. Performed the experiments: SB. Analyzed the data: SB. Wrote the paper: SB MW KL. Mathematical proofs: MW.

References

  1. 1. Newman MEJ (2003) The structure and function of complex networks. SIAM Rev 45: 167–256.
  2. 2. Boccaletti S, Latora V, Moreno Y, Chavez M, Hwang DU (2006) Complex networks: Structure and dynamics. Phys Rep 424: 175–308.
  3. 3. Arenas A, Díaz-Guilera A, Kurths J, Moreno Y, Zhou C (2008) Synchronization in complex networks. Phys Rep 469: 93–153.
  4. 4. Barrat A, Barthélemy M, Vespignani A (2008) Dynamical Processes on Complex Networks. New York, USA: Cambridge University Press.
  5. 5. Tsonis AA, Roebber PJ (2004) The architecture of the climate network. Physica A 333: 497–504.
  6. 6. Yamasaki K, Gozolchiani A, Havlin S (2008) Climate networks around the globe are significantly affected by El Niño. Phys Rev Lett 100: 228501.
  7. 7. Donges JF, Zou Y, Marwan N, Kurths J (2009) Complex networks in climate dynamics. Eur Phys J–Spec Top 174: 157–179.
  8. 8. Tsonis AA, Wang G, Swanson KL, Rodrigues FA, da Fontura Costa L (2010) Community structure and dynamics in climate networks. Clim Dynam. in press.
  9. 9. Steinhaeuser K, Chawla NV, Ganguly AR (2011) Complex networks as a unified framework for descriptive analysis and predictive modeling in climate science. Statistical Analysis and Data Mining 4: in press.
  10. 10. Abe S, Suzuki N (2004) Small-world structure of earthquake network. Physica A 337: 357–362. 23.
  11. 11. Abe S, Suzuki N (2006) Complex-network description of seismicity. Nonlinear Proc Geoph 13: 145–150.
  12. 12. Jiménez A, Tiampo KF, Posadas AM (2008) Small world in a seismic network: the California case. Nonlinear Proc Geoph 15: 389–395.
  13. 13. Krishna Mohan TR, Revathi PG (2011) Network of earthquakes and recurrences therein. J Seismol 15: 71–80.
  14. 14. Mantegna RN (1999) Hierarchical structure in financial markets. Eur Phys J B 11: 193–197.
  15. 15. Onnela JP, Kaski K, Kertesz J (2004) Clustering and information in correlation based financial networks. Eur Phys J B 38: 353–362.
  16. 16. Boginski V, Butenko S, Pardalos PM (2005) Statistical analysis of financial networks. Comput Stat An 48: 431–443.
  17. 17. Qiu T, Zheng B, Chen G (2010) Financial networks with static and dynamic thresholds. New J Phys 12: 043057.
  18. 18. Emmert-Streib F, Dehmer M (2010) Influence of the time scale on the construction of financial networks. PLoS ONE 5: e12884.
  19. 19. Reijneveld JC, Ponten SC, Berendse HW, Stam CJ (2007) The application of graph theoretical analysis to complex networks in the brain. Clin Neurophysiol 118: 2317–2331.
  20. 20. Bullmore E, Sporns O (2009) Complex brain networks: graph theoretical analysis of structural and functional systems. Nat Rev Neurosci 10: 186–198.
  21. 21. Kramer MA, Eden UT, Cash SS, Kolaczyk ED (2009) Network inference with confidence from multivariate time series. Phys Rev E 79: 061916.
  22. 22. Donges JF, Zou Y, Marwan N, Kurths J (2009) The backbone of the climate network. Europhys Lett 87: 48007.
  23. 23. Emmert-Streib F, Dehmer M (2010) Identifying critical financial networks of the DJIA: Toward a network-based index. Complexity 16: 24–33.
  24. 24. Erdős P, Rényi A (1959) On random graphs I. Publ Math Debrecen 6: 290–297.
  25. 25. Rao AR, Jana R, Bandyopadhyay S (1996) A Markov chain Monte Carlo method for generating random (0,1)-matrices with given marginals. Sankhya Ser A 58: 225–242.
  26. 26. Maslov S, Sneppen K (2002) Specificity and stability in topology of protein networks. Science 296: 910–913.
  27. 27. James R, Croft DP, Krause J (2009) Potential banana skins in animal social network analysis. Behav Ecol Sociobiol 63: 989–997.
  28. 28. Lima-Mendez G, van Helden J (2009) The powerful law of the power law and other myths in network biology. Mol Biosyst 5: 1482–1493.
  29. 29. Ioannides AA (2007) Dynamic functional connectivity. Curr Opin Neurobiol 17: 161–170.
  30. 30. Butts CT (2009) Revisiting the foundations of network analysis. Science 325: 414–416.
  31. 31. Bialonski S, Horstmann MT, Lehnertz K (2010) From brain to earth and climate systems: Small-world interaction networks or not? Chaos 20: 013134.
  32. 32. Antiqueira L, Rodrigues FA, van Wijk BCM, da F Costa L, Daffertshofer A (2010) Estimating complex cortical networks via surface recordings–a critical note. Neuroimage 53: 439–449.
  33. 33. Gerhard F, Pipa G, Lima B, Neuenschwander S, Gerstner W (2011) Extraction of network topology from multi-electrode recordings: Is there a small-world effect? Front Comp Neuroscience 5: 4.
  34. 34. Brillinger D (1981) Time Series: Data Analysis and Theory. San Francisco, USA: Holden-Day.
  35. 35. Pikovsky AS, Rosenblum MG, Kurths J (2001) Synchronization: A universal concept in nonlinear sciences. Cambridge, UK: Cambridge University Press.
  36. 36. Boccaletti S, Kurths J, Osipov G, Valladares DL, Zhou CS (2002) The synchronization of chaotic systems. Phys Rep 366: 1–101.
  37. 37. Kantz H, Schreiber T (2003) Nonlinear Time Series Analysis. Cambridge, UK: Cambridge Univ. Press, 2nd edition.
  38. 38. Pereda E, Quian Quiroga R, Bhattacharya J (2005) Nonlinear multivariate analysis of neurophysiological signals. Prog Neurobiol 77: 1–37.
  39. 39. Hlaváčková-Schindler K, Paluš M, Vejmelka M, Bhattacharya J (2007) Causality detection based on information-theoretic approaches in time series analysis. Phys Rep 441: 1–46.
  40. 40. Lehnertz K, Bialonski S, Horstmann MT, Krug D, Rothkegel A, et al. (2009) Synchronization phenomena in human epileptic brain networks. J Neurosci Methods 183: 42–48.
  41. 41. Watts DJ, Strogatz SH (1998) Collective dynamics of ‘small-world’ networks. Nature 393: 440–442.
  42. 42. Latora V, Marchiori M (2001) Efficient behavior of small-world networks. Phys Rev Lett 87: 198701.
  43. 43. Latora V, Marchiori M (2003) Economic small-world behavior in weighted networks. Eur Phys J B 32: 249–263.
  44. 44. Chung F, Lu L (2001) The diameter of sparse random graphs. Adv Appl Math 26: 257–279.
  45. 45. Press WH, Teukolsky SA, Vetterling WT, Flannery BP (2002) Numerical Recipes in C. Cambridge, UK: Cambridge University Press, 2nd edition.
  46. 46. Franaszczuk PJ, Bergey GK, Durka PJ, Eisenberg HM (1998) Time-frequency analysis using the matching pursuit algorithm applied to seizures originating from the mesial temporal lobe. Electroencephalogr Clin Neurophysiol 106: 513–521.
  47. 47. Schiff SJ, Colella D, Jacyna GM, Hughes E, Creekmore JW, et al. (2000) Brain chirps: spectrographic signatures of epileptic seizures. Clin Neurophysiol 111: 953–958.
  48. 48. Jouny CC, Franaszczuk PJ, Bergey GK (2003) Characterization of epileptic seizure dynamics using Gabor atom density. Clin Neurophysiol 114: 426–437.
  49. 49. Bartolomei F, Cosandier-Rimele D, McGonigal A, Aubert S, Regis J, et al. (2010) From mesial temporal lobe to temporoperisylvian seizures: A quantified study of temporal lobe seizure networks. Epilepsia 51: 2147–2158.
  50. 50. Schindler K, Bialonski S, Horstmann MT, Elger CE, Lehnertz K (2008) Evolving functional network properties and synchronizability during human epileptic seizures. Chaos 18: 033119.
  51. 51. Schindler K, Leung H, Elger CE, Lehnertz K (2007) Assessing seizure dynamics by analysing the correlation structure of multichannel intracranial EEG. Brain 130: 65–77.
  52. 52. Schreiber T, Schmitz A (1996) Improved surrogate data for nonlinearity tests. Phys Rev Lett 77: 635–638.
  53. 53. Schreiber T, Schmitz A (2000) Surrogate time series. Physica D 142: 346–382.
  54. 54. Roberts JM (2000) Simple methods for simulating sociomatrices with given marginal totals. Soc Networks 22: 273–283.
  55. 55. Maslov S, Sneppen K, Zaliznyak A (2004) Detection of topological patterns in complex networks: correlation profile of the internet. Physica A 333: 529–540.
  56. 56. Artzy-Randrup Y, Fleishman SJ, Ben-Tal N, Stone L (2004) Comment on “Network Motifs: Simple building blocks of complex networks” and “superfamilies of evolved and designed networks”. Science 305: 1107.
  57. 57. Milo R, Itzkovitz S, Kashtan N, Levitt R, Alon U (2004) Response to comment on “Network Motifs: Simple building blocks of complex networks” and “superfamilies of evolved and designed networks”. Science 305: 1107.
  58. 58. Artzy-Randrup Y, Stone L (2005) Generating uniformly distributed random networks. Phys Rev E 72: 056708.
  59. 59. Del Genio CI, Kim H, Toroczkai Z, Bassler KE (2010) Efficient and exact sampling of simple graphs with given arbitrary degree sequence. PLoS ONE 5: e10012. 27.
  60. 60. Blitzstein J, Diaconis P (2010) A sequential importance sampling algorithm for generating random graphs with prescribed degrees. Internet Mathematics 6: 489–522.
  61. 61. Eguiluz VM, Chialvo DR, Cecchi GA, Baliki M, Apkarian AV (2005) Scale-free brain functional networks. Phys Rev Lett 94: 018102.
  62. 62. van den Heuvel MP, Stam CJ, Boersma M, Hulshoff Pol HE (2008) Small-world and scale-free organization of voxel-based resting-state functional connectivity in the human brain. Neuroimage 43: 528–539.
  63. 63. Hayasaka S, Laurienti PJ (2010) Comparison of characteristics between region-and voxel-based network analyses in resting-state fMRI data. Neuroimage 50: 499–508.
  64. 64. Fransson P, Åden U, Blennow M, Lagercrantz H (2011) The functional architecture of the infant brain as revealed by resting-state fMRI. Cereb Cortex 21: 145–154.
  65. 65. Tian L, Wang J, Yan C, He Y (2011) Hemisphere- and gender-related differences in small-world brain networks: A resting-state functional MRI study. Neuroimage 54: 191–202.
  66. 66. Guye M, Bettus G, Bartolomei F, Cozzone PJ (2010) Graph theoretical analysis of structural and functional connectivity MRI in normal and pathological brain networks. Magn Reson Mater Phy 23: 409–421.
  67. 67. Small M, Yu D, Harrison RG (2001) Surrogate test for pseudoperiodic time series data. Phys Rev Lett 87: 188101.
  68. 68. Breakspear M, Brammer M, Robinson PA (2003) Construction of multivariate surrogate sets from nonlinear data using the wavelet transform. Physica D 182: 1–22.
  69. 69. Nakamura T, Small M (2005) Small-shuffle surrogate data: Testing for dynamics in fluctuating data with trends. Phys Rev E 72: 056216.
  70. 70. Keylock CJ (2006) Constrained surrogate time series with preservation of the mean and variance structure. Phys Rev E 73: 036707.
  71. 71. Suzuki T, Ikeguchi T, Suzuki M (2007) Algorithms for generating surrogate data for sparsely quantized time series. Physica D 231: 108–115.
  72. 72. Romano MC, Thiel M, Kurths J, Mergenthaler K, Engbert R (2009) Hypothesis test for synchronization: Twin surrogates revisited. Chaos 19: 015108.