II. Denitions and examples

The basic steps that have been used [1{4] to dene the point process in brain signals are summarized in Fig. 1. The raw data consist of time series recorded from the brain using functional magnetic resonance imaging (fMRI) corresponding to the activity of one of many thousands of small brain regions. It is accepted that this imaging technique measures a \blood oxygen level-dependent" signal (i.e., \BOLD") in each small region, giving an estimation of the blood oxygen saturation, which itself is proportional to local neuronal activity. The point process can be dened in diferent ways, but for the reasons discussed later, the end results are equivalent. As shown in Fig. 1, time points can be selected at the upward threshold crossings (here at unity) of the signal (lled circles). A second approach is to construct the point process by selecting the local peaks and/or valleys of the BOLD time series. For ease of discussion, we will deal with the second option in this paper. The temporal co-occurrence of the points denes the co-activation matrix (bottom graphs), which can be further averaged to estimate the correlation matrix of the system under study. Figure 2 illustrates, for those unfamiliar to the subject, three examples of typical BOLD time series that are usually recorded from the brain. From visual inspection it is already apparent that they are smooth traces, exhibiting temporal correlations, as will be further discussed later. There are also spatial correlations, for instance between the top two traces, which is evidenced in the images' heat map between counter-lateral regions.

A qualitative comparison of how well it works:

It has already been established, in diferent circumstances [2{4], that the co-activation matrix otained with the point process method is very similar to the correlation matrix computed from the full (i.e., continuous) BOLD signal. Figure 3 shows an example of a correlation matrix constructed from the point process computed from a subject while resting (data fully described in [4]). The results demonstrate that as few as 4 points are already sufcient to dene clusters of co-activation, as demonstrated previously in [2{4]. In addition, the results here show how de-activations (i.e., blueish colors) can also be evidenced by the PP approach.

Figure 2: Right: Examples of typical BOLD time series from three selected sites in the brain. The dashed box in the middle time series indicates a portion of the signal used later in the analysis presented in Fig. 4. Left: Images correspond to a snapshot of activity amplitudes (top), the correlations between the three selected seeds (red dotted circles) and the rest of the brain (middle three images), and the corresponding brain structural slice at MNI coordinate z=10 (bottom). The MNI x, y, z coordinates are -31 -95 10, respectively, for the top trace, 43 -73 10 for the middle trace and -10 -31 10 for the bottom trace. 

III. Why does it work? A simple theory

As discussed previously, the example in Fig. 3 implies a large compression; the question then is why a few points are enough to compute results similar to those obtained with the full signal. A simple visual inspection of the BOLD traces reveals that the type of signals we are dealing with are temporally correlated. This is very well known; the neuronal activity is temporally and spatially correlated, and furthermore, the activity is convoluted by the hemodynamic transfer function which in itself introduces additional temporal correlations.

Therefore, for any time series with these properties, it seems natural to think that the most informative points are those in which its derivative changes sign. The other points are redundant since they can be predicted, to a certain degree, by a linear estimator. This is illustrated in Fig. 4, using as an example two minutes of BOLD recording (normalized by its standard deviation ). After setting a threshold , the in ection points larger than a given value are identied. These points constitute the marked point process in question.

Figure 3: Example of correlation maps obtained from the raw BOLD time series of length n=235 (right panel) and from the derived point process (left panels) for different numbers of points (n=4,7,14,26). The left images represent, as \heat maps", the co-activation of the seed (located at MNI coordinates x=4, y=-60, z=18) with respect to each voxel. Note that a few points already suce to identify well-dened clusters that are 1-4 standard deviations away from chance coctivations. The right panel corresponds to the Pearson correlation computed from the entire length of the raw BOLD time series. Red/blue colors label positive/negative point co-activations in the case of the left maps and positive/negative correlations in the case of the right map. 

Now we ask how much of the raw signal is left out if these in ection points are used to extrapolate a piece-wise linear time-series. To answer this we analyze BOLD time series from the brain of a subject during an experiment in which fMRI data are collected at rest [3]. We proceed to compute the linear correlation between the two time series, the raw and the piece-wise linear one. In panels D and E the results are shown for diferent values of threshold (in units of ) as well as for the correlation of the time series, estimated by the value of the rst autocorrelation coecient. Panel D shows that as the BOLD signal autocorrelation increases, the similarity between the piece-wise linear and the raw signals increases, evaluated in two ways: by the root mean squared error (rmse) and by the linear correlation hri between the two time series. As expected, raising the threshold above zero produces an increasing loss of information about the signal, which is re ected in a monotonic increase in the rmse and a decrease in the hri values (see Panel E).

Figure 4: Why it works: The trace in Panel A is an example of a two-minute recording of a BOLD brain signal during rest (denoted by the dotted box in Fig. 2). The point process is dened by the timing of the peaks and valleys larger than a given threshold, (two are indicated here by arrows). The points, in this case, are only six (dots depicted in the bottom trace) out of the 120 samples of the original time series. The ability of these six points to preserve information about the original BOLD signal can be estimated by its similarity with a piece-wise linear time series (dashed red lines) constructed by joining the peaks and valleys. Horizontal dotted lines in Panel A denote the threshold for the example. Panels D and E correspond to the computed similarity between the BOLD signals and the piece-wise linear signals (evaluated by the correlation hri and rmse values) for diferent autocorrelation and threshold values. Panels B and C correspond to similar calculations using synthetic time series. For Panels B and D was xed at 1. For Panels C and E was 0.85. 

According to the present hypothesis, the functional dependence shown by the BOLD signals in Panels D and E will be replicated using synthetic signals with similar autocorrelation properties. To this end, we generate articial time series with autocorrelation values identical to those of the BOLD signals, using the MATLAB routine f_alpha_gaussian.mfrom [16] (see also source codes at [¹⁷]). Panels B and C show that the behavior with respect to the threshold and for the synthetic and empirical data are very similar.

The results show that the key to understanding why the approach works lies in the correlation properties of the time series under consideration. In synthesis, it is found that signals with long-range correlations are particularly suitable for this type of compression, where in ection points contain most of the information. The results also apply to other signals from any origin, as long as their autocorrelation features are similar.

IV. Further testing using synthetic ground-truth data.

The results in the previous section emphasize the relevance of the individual signal's autocorrelation as the main property related to the ability of a point process to preserve information about the original signal, and consequently, about functional connectivity between signals. We can test this by manipulating the autocorrelation in any given system for which the \ground-truth" crosscorrelations are known. The data reported by Smith et al. [19] can be used for our purpose

The authors in [19] reviewed and compared the available fMRI analysis methods ranging from simple measures of pair-wise linear correlations to sophisticated multivariate approaches. In the process, they generated diverse and realistic simulated fMRI data sets describing dierent underlying networks. These simulations were based on the dynamic causal modelling (DCM) [20] fMRI forward model (see [19] for full details). We used these simulated BOLD time series (downloaded from the authors' site [22]) to test the point process approach in comparison with standard correlation methods. First, for completeness, we brie y describe the essence of the model used in these simulations. The model: Smith et al. simulations used a neural network model coupled with Buxton's nonlinear balloon model [21] for the vascular dynamics.

Each neural network node has a binary external input (square symbols in the top diagram of Fig. 5). The state of the inputs (i.e., active or inactive) is given by a Poisson process which controls the probability of switching states, and can be seen as external signals or as noise at the neural level. Subsequently, the neural signals propagate across the network according to the DCM neural network model, where node interactions are dened by the A network matrix:

Figure 5: Simulated fMRI network from Ref. [19]. Panel A depicts the topology and Panel B the adjacency matrix of the interactions of the nodes, where negative values (labeled red) correspond to self-interactions. Connections are directed: a node in the upper diagonal of the matrix denotes a directed connection from a lower-numbered node to a higher-numbered one. Panel C shows an example of the BOLD time series simulated on this network. 

where z is the neural time series, z is its rate of change, u are the external inputs and C the weights controlling how the external inputs feed into the network (here just through the identity matrix). The o-diagonal terms in A determine the network connections between nodes (arrows in Fig. 5A), and the diagonal elements are all set to -1 to model within-node temporal decay; in this way, controls both the within-node (neural) temporal inertia and the temporal lag between nodes.

The time series of the activity for each node is then fed to a nonlinear balloon model for vascular dynamics [21] output, which is a function of the changing neural activity. The parameters were adjusted by the authors in order to match BOLD time series seen for typical data of brain resting activity recorded with 3 Tesla fMRI technology. Finally, various sources of variability were added to account for realistic expectations. The BOLD time series and the underlying ground-truth network matrices are accessible on the authors' site [22]. In the following paragraphs, we will explore the ability of the point process to extract the underlying network and to compare it with the commonly-used correlation method. Figure 5 shows the ground-truth network considered here (le sim4.mat downloaded from the site [22]). Panels A and B illustrate the topology and the adjacency matrix of the network (A in Eq. 1). It comprises 10 regular modules interconnected by a few links, a typical small-word graph; a total of 50 nodes interconnected by 61 positive odiagonal interactions (40 of which correspond to nearest neighbors), as well as 50 negative (diagonal) self-interactions. Panel Cshows typical traces of the computed BOLD time series recorded at the 50 nodes, simulating data from a fMRI typical session. The data set contains 50 stochastic realizations representing simulated fMRI records of diferent human subjects.

Figure 6: Panels A and B: Comparison of the partial correlation matrices {one calculated with the original BOLD data (Raw) and the other from its derived point process (PP)- as a function of the threshold used to dene the point process. Panel A corresponds to the root mean squared values of the dierences, and panel B to the correlation coecients. Dots correspond to individual results while averages (lled circles) and error bars correspond to the mean and standard deviation of the 50 simulated BOLD records (time series length T=200 samples) Panel C shows the Receiver Operating Characteristic curve obtained for both methods in detecting the presence of links ( = 0:7 and = 2:2) computed for two time series of maximum length (T. Max. indicated in the legend). Panel D corresponds to the area under the ROC curve, computed for both methods as a function of the increasing maximum length, T. Max., of the time series considered. 

Extracting the correlation graph: To benchmark the relative merits of the point process approach, we rst extracted the point process from the BOLD time series, and then computed the covariance matrices for both the point process and the raw BOLD dataset. Following this, the partial correlation from both matrices was calculated (as in [19]), and their dierences compared for various levels of the threshold used to dene the point process. As seen in Figs. 6A and 6B, there was an optimum threshold (for this dataset = 0:7 ) at which the correlation matrices became more similar. We then proceeded to check how well the method performed in predicting the underlying connectivity graph from the time series. Specically, we checked how well the correlation matrices from both the point process and the raw BOLD dataset described the o-diagonal elements depicted in Fig. 5B, i.e., the synthetic ground-truth network. We used the receiver operating characteristic curve (ROC) [23], which benchmarks specicity and sensitivity as a function of a given parameter. To determine whether a connection is predicted or not between two nodes, we chose a threshold of 2 (corresponding to P = 95%) at the (i; j) partial correlation matrices entry. In general, to obtain the ROC curve a given range of relevant parameters must be explored while the true/false positive/negative predictions are counted. Here, for convenience, we chose to explore a range of time series lengths from relatively very short (set to 20 samples here) up to a variable maximum length, T.Max. For the examples presented, the values of T.Max. ranged from 400 to 2000 samples.

Figure 6C illustrates the results, where the family of curves (triangles for raw data and circles for the point process) corresponds to time series of various maximum lengths (T.Max.). As expected, the shortest T.Max (400 samples) gave the lowest condent results, while the longest T.Max. (2000 samples) resulted in a very good estimation of the true network connections. The area under the curve, plotted in Fig. 6D, is a good estimation of the relative goodness of the prediction, where a value of 1 corresponds to a perfect prediction and a value of 0.5 is equivalent to chance. Note that the main motivation for these numerical simulations is to demonstrate that there is close similarity between the PP and Raw ROC curves in Fig. 6D. This similarity is indicative of the good performance of the point process approach compared with the use of the raw time series.