Revisiting correlation-based functional connectivity and its relationship with structural connectivity

Patterns of brain structural connectivity (SC) and functional connectivity (FC) are known to be related. In SC-FC comparisons, FC has classically been evaluated from correlations between functional time series, and more recently from partial correlations or their unnormalized version encoded in the precision matrix. The latter FC metrics yield more meaningful comparisons to SC because they capture ‘direct’ statistical dependencies, that is, discarding the effects of mediators, but their use has been limited because of estimation issues. With the rise of high-quality and large neuroimaging datasets, we revisit the relevance of different FC metrics in the context of SC-FC comparisons. Using data from 100 unrelated Human Connectome Project subjects, we first explore the amount of functional data required to reliably estimate various FC metrics. We find that precision-based FC yields a better match to SC than correlation-based FC when using 5 minutes of functional data or more. Finally, using a linear model linking SC and FC, we show that the SC-FC match can be used to further interrogate various aspects of brain structure and function such as the timescales of functional dynamics in different resting-state networks or the intensity of anatomical self-connections.

where B encodes anatomical connections (Eq. (3)), y t denotes observed brain function at time t, x t is the driving noise at time t, and we have y 0 = 0. Let us now assume that x t evolves in a stepwise fashion, with a certain unknown step h, i.e., x t = x k for all t 2 [(k 1)h, kh] where the sequence of random values {x k } fulfills the (stepwise) whiteness condition E[x k x > k 0 ] = s 2 x d kk 0 . Let us now focus on a time epoch t = nh + D, where n is an integer and 0 < D < h. Under the aforementioned model, we can apply (S1) to obtain: where m = n k and we have used the fact that B 1 and e Bt commute. Using the whiteness condition, we can now compute the covariance of y t , denoted by S(D), as: where we have repeatedly applied commutativity of the pertinent matrices, the cross-terms are annihilated due to the whiteness condition, and we have considered n sufficiently large in order for the sum of the following geometric matrix progression to be expressed as: Â n 1 m=0 e 2Bmh ⇡ Â • m=0 e 2Bmh = (I e 2Bh ) 1 , where we have also used the fact that matrix B is negative semidefinite. The average covariance of the observed signal can then be obtained by computing: In order to explore how the timescales of the driving noise dynamics, i.e., the value of h, impacts the observed covariance matrix S, let us first consider the case where h is small. In

Network Neuroscience 1
Liegeois, R., Santos, A., Matta, V., Van De Ville, D., & Sayed, A. H. (2020). Supporting information for "Revisiting correlationbased functional connectivity and its relationship with structural connectivity." Network Neuroscience, 4(4), 1235-1251. https://doi.org/10.1162/netn_a_00166 that case, we can use a second-order Taylor approximation to rewrite Eq. (S4) as: When h is large, Eq. (S4) suggests: In other words, it can be seen from Eqs. (S5) and (S6) that when the driving process exhibits fast dynamics, i.e., when h is small, then S µ B 1 , or equivalently B µ S 1 , whereas when the driving process exhibits slow dynamics, we have S µ B 2 , or equivalently B µ S 1 2 . These findings are illustrated using a toy example in Figure S1.
Another way of exploring the timescales of driving and observed dynamics is to reformulate Eq. (2) as:ẏ where a is a positive scalar and higher values of a reflect faster observed dynamics y t . Assuming white Gaussian driving noise x t , it can be shown that B µ S 1 a (e.g., Oku & Aihara, 2018). Hence observing optimal exponents of the SC-FC match closer to 0 in Figure  4 reflects faster timescales of the observed dynamics y t as compared to the driving noise x t .
This interpretation complements and is consistent with Eqs. (S5)-(S6). Indeed, Eq. (S7) suggests that faster processing dynamics of y t (i.e., higher exponent of B), given a fixed driving noise x t lead to an optimal exponent of S that gets closer to 0. On the other hand, Eqs. (S5)-(S6) suggest that slower input dynamics x t (i.e., with longer steps h) for fixed processing dynamics y t lead to an optimal exponent of S closer to 0.5 (coming from 1). So in both cases having an optimal exponent that departs from 1 and goes towards less negative values corresponds to having a decreased ratio between timescales of y t and x t : in one case, the timescale of y t decreases (y t faster), in the other, the timescale of x t increases (x t slower).

Optimal exponent for intermediate values of h
In the case of fast input dynamics, we have seen in Eq. (S5) that B µ S 1 , and hence the optimal exponent b is -1. When the input gets slower (i.e., larger values of h), B µ S 2 and hence the optimal exponent b is -1/2 (Eq. (S6)). To explore the value of the optimal exponent for intermediate values of h, we used a toy example with M = 10 variables, T = 100000 time points, a minimal step size dt = 0.005, and a random (symmetric and stable) structural connectivity matrix B. Figure S1 shows the optimal value of the exponent b for different values of h, averaged over 200 trials.
Network Neuroscience Figure S1. Impact of step size on the value of the optimal exponent.
It can be seen that as expected in the cases of low and high values of h (i.e., h = 1 ⇤ dt and h = 50 ⇤ dt), the optimal exponent is found to be equal to -1 and -1/2, respectively. Moreover, for intermediate values of h, this toy example confirms that the optimal exponent is found to smoothly go from -1 to -1/2 when the value of h increases.

Regularized versions of the precision matrix
Inverting the correlation matrix might be an unstable operation especially when it is computed from few time points (Brier et al., 2015). Therefore, we considered a regularized version of the precision matrix. The regularization used in the results of the main manuscript (Figures 2&3) relies on the Tikhonov regularization. This consists in adding a full-rank regularization term to the correlation matrix before performing the inversion. In other words, the Tikhonov regularized precision matrix is equal to (S + g · I), where S is the empirical correlation matrix, I is the identity matrix and g is a parameter that was optimized by minimizing the distance between regularized subject precision matrices and the unregularized group precision matrix. More precisely, we followed the procedure described in Pervaiz et al. (2020) that consists in (i) for a given value of g compute the N regularized subject precision matrices, (ii) compute the distance between these N matrices and a unique unregularized group precision matrix (after taking the square of all matrix entries), and (iii) choose the value of g that minimizes the sum of all N distances. The parameter g was varied between 0 and 1 with 0.02 increments, and its optimal value for different values of functional time series is shown in Figure S2. As expected, when more time points are used to compute the correlation matrix, the optimal value of g decreases since less regularization is required. These optimal values of g were considered when using less than one hour of functional data to compute regularized precision matrices in Figure 2C. When using more than one hour of data, time series from different subjects had to be concatenated and in that case we considered g = 0.32 to compute regularized precision matrices in Figure 2C.

Computing SC-FC match using variations of the SC matrix
When exploring the link between SC and FC, the match between these two matrices is often computed only over non-zero SC entries (e.g., Honey et al., 2009). We did not restrict our analyses to non-zero entries mainly because it requires to set a threshold under which an SC entry is considered to be zero and we wanted our results to be as general as possible. To explore the impact of considering such a threshold, we repeated in Figure S3B,C the analysis of Figure 2C but considering only non-zero SC entries to evaluate the SC-FC match. Then, the results shown in the main manuscript were obtained using SC matrices with entries normalized by the volumes of starting and ending ROIs. To further support robustness of our results, we also repeated the analysis of Figure 2C using SC matrices with entries that were not normalized by the corresponding ROIs' volumes ( Figure S3D). SC-FC match for four FC metrics as a function of data duration to compute FC metrics. (A) The match is computed using all SC and FC entries, similarly to what is presented in Figure 2C, (B) the match is computed using only SC (and corresponding FC) entries that are higher than 0.05, (C) the match is computed using only SC (and corresponding FC) entries that are higher than 0.2, and (D) the match is computed using SC with entries that were not normalized by starting and ending ROIs' volumes.
It can be seen from Figure S3B,C that our main conclusions hold even when computing the SC-FC match only over non-zero SC entries, using two different thresholds. In particular, regularized precision-based FC still outperforms other FC metrics when few data points are available. Interestingly, there is an overall increase of the correlation-based SC-FC match when only non-zero SC entries are considered: using one hour of data, this match is equal to 0.24 when all entries are used ( Figure S3A), it is equal to 0.37 when considering only SC entries higher than 0.05 ( Figure S3B), and it is equal to 0.38 when considering only SC entries higher than 0.20 ( Figure S3C). Finally, Figure S3D shows that the main results are reproduced when considering unnormalized SC matrices, even if the overall SC-FC match is lower as compared to other cases. This suggests that normalizing SC entries by corresponding ROIs' volumes is a meaningful operation in the context of SC-FC comparisons.

Linking the optimal exponent and behavioral measures
The procedure described in Figure 4 allows to identify the correlation matrix exponent that provides the best match to the SC matrix. As explained in the Supplementary Methods (Eqs. (S1)-(S7)), this metric provides a measure of functional dynamics timescales. We here tested whether this value, computed for each subject, is related to eight HCP behavioral measures including bodily features (age, weight, BMI), arousal measures (sleep quality), cognitive features (fluid intelligence, visual episodic memory, cognitive flexibility), and task performance (relational processing). Figure S4. Correlation between the optimal correlation matrix exponent of each subject, and eight HCP behavioral measures: age, weight, BMI, and sleep quality (PSQI score), fluid intelligence (PMAT24), visual episodic memory (PicSeq Unadj), cognitive flexibility (CardSort Unadj), relational processing (relational task accuracy) and sleep quality (PSQI score).
It can be seen from Figure S4 that none of these measures seems to be significantly correlated with the optimal correlation matrix exponent. Note that we have explored the link with other behavioral measures classically used in such applications (e.g., Li et al., 2019;Liégeois et al., 2019), and did not find links that survived correction for multiple comparisons. Future work might explore whether the optimal exponent in specific networks is a better behavioral marker, or if multivariate approaches (e.g., canonical component analy-sis) might help revealing a link between functional timescales as captured by the optimal exponent and behavior.

Impact of repetition time (TR)
In order to evaluate the impact of the repetition time, we reproduced the analysis of Figure  2C by removing one functional time point every two functional time points (corresponding to an effective TR = 1.44 sec), and two functional time points every three functional time points (correponding to an effective TR = 2.16 sec).

Figure S5.
Impact of TR on the (average) match between SC and FC evaluated from the precision matrix (dark blue), the regularized precision matrix (light blue), the correlation matrix (red) and the autoregressive matrix (green) of unfiltered fMRI time-series. We considered TR = 0.72 sec, TR = 1.44 sec, and TR = 2.16 sec.
The results are shown in Figure S5 and suggest that using the same number of fMRI time points but with a larger repetition time (TR), precision-based FC and correlation-based FC estimates have a better match to SC whereas AR-based FC is less correlated with SC. This is explained by the fact that fMRI time points carry less redundant information when the sampling period TR increases, which is beneficial for the estimation of the correlation and precision matrices, but penalizes the estimation of the AR-based metric because it exploits autocorrelation of time series. Obviously, considering longer repetition times also comes at the price of longer scanning sessions for a given number of time points.

Subject-level optimal exponent using short acquisition times
The results of Figure 4 were computed using data concatenated over all subjects. In practice, such amount of data is not always available and we here provide the corresponding curves when using only the first HCP run (i.e., 1200 time points) for the three first HCP subjects. Since matrix inversion operations might be unstable when using such small amount of data, we used the regularization procedure described in Figure S2 with g = 0.52 which is the optimal regularization parameter when using 1200 time points. It can be seen that when using whole-brain data (left panels in Figure S6), subject-level curves are similar to the group-level curve shown in Figure 4A even if the optimal ex-ponents seem to have smaller absolute values. Within resting-state networks, the overall shape also seems to be similar at the subject-level (right panels in Figure S6) as compared to the group level ( Figure 4A), with common features such as a high SC-FC match in the visual network and a low SC-FC match in the limbic network. Future work might explore distribution of these curves across subjects and test whether they encode behavioral information, as presented in Figure S4 at the whole-brain level. Figure S7 shows different average structural and functional connectomes (upper row), together with the corresponding distributions of entries (lower row).

Impact of averaging effects
One question we might have when observing results of Figure 2C is whether averaging functional and structural connectomes over different subjects might boost the SC-FC match. To test this, we computed this match using four runs of functional data (⇠ 1 hour) in three cases: (i) the four runs come from the same subject and are compared to the subject's SC matrix (as in Figure 2B), (ii) the four runs come from two different subjects chosen at random (first two runs are the first two runs of one subject, and last two runs are the last two runs of the other subject) and are compared to the corresponding averaged SC matrix, and (iii) the four runs come from four different subjects chosen at random and are compared to the corresponding averaged SC matrix. Experiment (i) was performed for all 100 subjects, and experiments (ii) and (iii) were repeated 100 times. SC-FC match using four runs of functional data in three cases: (i) the four runs come from the same subject, (ii) the four runs come from two different subjects chosen at random, and (iii) the four runs come from four different subjects chosen at random. Figure S8 indeed suggest that averaging functional and structural connectomes slightly increases the SC-FC match as compared to the case where data coming from one subject is used, which might partly explain the increase observed in Figure 2C when using data from several subjects. We finally note that our main findings, i.e. the advantage of (regularized) precision-based FC over other FC metrics in the context of SC-FC comparisons, is reproduced in the three configurations presented in Figure S8.