To compare different TMFC methods, we first used empirical fMRI data from the Human Connectome Project (HCP)39 and the Consortium for Neuropsychiatric Phenomics (CNP)40. In particular, we considered two block design tasks (working memory and social cognition tasks, N = 100) from the HCP dataset39,41, two event-related tasks (stop-signal and task-switching tasks, N = 115) from the CNP dataset40,42 and resting-state data from both datasets. To construct empirical FC matrices, we used a set of functionally defined regions of interest (ROIs) covering the whole brain43. The ROIs were defined as spheres with a radius of 4 or 5 mm. We discarded ROIs for which data were incomplete for at least one subject. As a result, we utilised 239 ROIs for the HCP dataset and 246 ROIs for the CNP dataset.Next, we performed a series of simulations to compare the sensitivity and specificity of TMFC methods in experiments with block and event-related designs, different SNRs, sample sizes, duration of events, mean interstimulus intervals (ISIs) and number of events. For block designs, we considered the CorrDiff approach and PPI methods, and for event-related designs, we considered the PPI methods and BSC methods. We also compared the effectiveness of fast data acquisition for a fixed number of scans and fixed total scan time. Then, we assessed the impact of HRF variability on the sensitivity of TMFC methods. Finally, we performed simulations with symmetric (undirected) and asymmetric (directed) ground-truth synaptic matrices to identify sources of PPI matrix asymmetry.Empirical comparison of TMFC methodsWe found that all TMFC methods, except cPPI, produce similar unthresholded matrices for the block design (working memory task, Fig. 2a). The correlations between CorrDiff matrices and symmetrised sPPI and gPPI matrices ranged from 0.66 to 1. We also found that the correlations between these TMFC matrices and RSFC, BGFC, TSFC and cPPI matrices ranged from −0.13 to −0.28.Fig. 2: Correlations between unthresholded RSFC, BGFC and TSFC matrices and FC matrices obtained by different TMFC methods.To evaluate the similarity between the raw FC matrices, we calculated Pearson’s r correlations between lower diagonal elements. sPPI and gPPI matrices were symmetrised. All PPI terms were calculated with the deconvolution step. a Results for the block design: working memory task. b Results for the event-related design: stop-signal task.For the event-related design (stop-signal task), similar TMFC matrices were produced by the sPPI, gPPI, BSC-LSS and BSC-FRR methods, with correlation coefficients ranging from 0.63 to 0.96 (Fig. 2b). At the same time, correlations between these TMFC matrices and the RSFC, BGFC, TSFC and cPPI matrices ranged from −0.09 to 0.004. Remarkably, the BSC-LSA approach produced a random-like matrix correlated with all other FC matrices with a correlation of 0.01–0.09. The similarities and differences between these FC methods are also confirmed by calculating the overlap between thresholded matrices (Supplementary Figs. S1, S2). Analogous results were obtained for empirical data from other tasks with block and event-related designs (social cognition task and task-switching task), which demonstrates that the observed effects are specific to a design type rather than a particular task (Supplementary Figs. S3, S4).Notably, the matrices generated by the cPPI method are similar to matrices obtained using the RSFC, BGFC, and TSFC approaches, which mainly reflect task-independent spontaneous fluctuations10. Correlations between the cPPI, RSFC, BGFC and TSFC matrices ranged from 0.69 to 0.99 for block and event-related design (Fig. 2a, b). The partial correlation between two PPI terms controlling for physiological (BOLD signal in both regions) and psychological variables (task regressor), as implemented in cPPI, was as high as the simple correlation between these two PPI terms and the correlation between physiological regressors (Fig. 3). Therefore, we concluded that the cPPI approach is unable to separate task-modulated (extrinsic) from task-independent (intrinsic) sources of FC. However, comparison of FC matrices based solely on empirical data does not allow us to say which method better reflects the true TMFC. To answer this question, we applied a biophysically realistic simulation approach that enabled us to control the ground-truth TMFC.Fig. 3: The cPPI approach produces FC matrices similar to the RSFC and TSFC matrices.a Correlations for the block design task (working memory task) and resting-state data from the HCP dataset. b Correlations for the event-related design task (stop-signal task) and resting-state data from the CNP dataset. c Illustration of the cPPI approach. The working memory task time series were taken from the two ROIs as an example. All PPI terms were calculated with the deconvolution step. To evaluate the similarity between matrices, we used Pearson’s r correlation. The color scales were adjusted for each matrix based on the maximum absolute value and were assured to be positive and negative symmetrical.Large-scale neural mass simulationsOur simulation approach was based on the coupled oscillator model for FC proposed by Mateo et al. (2017)44. Using optogenetic manipulations and concurrently measuring local field potential, arteriole diameter and blood oxygenation in wake mice, they showed that correlations between ultra-slow BOLD fluctuations (i.e., FC measured by fMRI) are caused by synchronised ultra-slow fluctuations in arteriole diameter, which in turn are caused by ultra-slow modulation of the envelopes of synchronised gamma-band oscillations (Fig. 4a). Accordingly, we modelled gamma-band activity of 100 interconnected regions using the large-scale Wilson-Cowan model, and also modelled arteriole diameter dilations and blood oxygenation change using the Balloon-Windkessel model (Fig. 4b). For details about the simulation procedures, see Methods. The matrices of gamma-band neuronal synchronisation measured by the phase-locking value were closely matched to synaptic weight matrices (Fig. 4c, d). In accordance with previous empirical observations by Mateo et al. (2017)44, we observed strong correlations between simulated ultra-slow fluctuations of the gamma-band envelope and time-shifted BOLD signal. For the event-related design, the simulated BOLD signal was correlated with the gamma-band envelope at r = 0.79 with a 3.5 s time lag (Fig. 4e–j). Similar results were obtained for the block design, where the correlation was r = 0.81 with a 3.5 s time lag (Supplementary Fig. S5).Fig. 4: The emergence of ultra-slow BOLD signal fluctuations from ultra-slow modulations of fast oscillatory activity.a Coupled oscillator model for functional connectivity proposed by Mateo et al.44. Our simulation approach was based on the empirical evidence44 that correlations between ultra-slow BOLD fluctuations are caused by ultra-slow modulation of the envelopes of synchronised gamma-band oscillations. b Scheme for generating the BOLD signal from gamma-band oscillations simulated by the Wilson-Cowan model. c Ground-truth symmetric synaptic weight (wji) matrices for a single subject. For the long-range connections between Wilson-Cowan units (wji), we used three synaptic weight matrices corresponding to two task conditions (“Cond A” and “Cond B”) and interim “Rest” periods. Each synaptic weight matrix consisted of 100 brain regions, which corresponded to the minimum number of areas in brain atlases125, and four functional modules, given the presence of such modules in the human brain5. The highest synaptic weights were within each module during the “Rest” periods5. Under the task conditions, we increased synaptic weights between functional modules relative to the “Rest” periods8. d Gamma-band neuronal synchronisation estimated by the phase-locking value (ρji) for a single subject. The simulation was performed for the event-related design with one hundred 1 s events per condition and mean ISI = 6 s. e Example of simulated synaptic activity (SA) and gamma-band envelope for one of 100 connected brain regions (ROIs). f The time series of the gamma-band envelope and BOLD signal generated by the Balloon-Windkessel model based on simulated SA. We used the standard parameters of the haemodynamic model26, which have previously been used in whole-brain RSFC simulation studies25,126. g One second of simulated SA. h Power spectral density (PSD) of simulated SA. The main peak at 40 Hz (gamma-band oscillations). i Cross-correlation between the gamma-band envelope and BOLD signal. The maximum correlation r = 0.79 corresponds to a time lag of 3.5 s. j Power spectral density of the BOLD signal.The cPPI method fails to estimate TMFCWe first considered simulations without co-activations to investigate whether different TMFC methods produce FC matrices similar to ground-truth synaptic weight matrices for a sample size N = 100, SNR = 0.4, and TR = 2 s. For differences between simulated conditions (“Cond A-B”), almost all tested methods produced TMFC matrices similar to the ground truth (Fig. 5a–g). In contrast, the cPPI method produced matrices similar to BGFC and TSFC (Fig. 5h), which replicated the above-mentioned findings from the empirical data analysis (Fig. 2). Knowing the ground-truth FC patterns, we revealed that the cPPI method, unlike other TMFC methods, does not show the modulation of FC between compared task conditions (the “Cond A-B” effect), but rather shows the sum of all FC during both conditions, that is, the “Cond A + B” effect (Supplementary Fig. S6). Since task-unrelated FC is present in both conditions, we see both task-unrelated and task-related effects in the cPPI matrix. Other TMFC methods effectively remove task-unrelated effect by subtracting FC in one condition from another. Thus, the cPPI method does not contrast FC in one condition relative to another and does not remove task-unrelated effects (functional connections that is high at rest). Therefore, we excluded the cPPI method from further analysis, as it is unable to assess TMFC.Fig. 5: Illustration of FC matrices produced by different TMFC methods based on simulated data.An example of simulation results for sample size N = 100, SNR = 0.4 and TR = 2 s. a Expected FC matrices based on ground-truth synaptic weight matrices. b–d For the CorrDiff, sPPI and gPPI methods, we considered the block design with ten 20 s blocks per condition. e–g For the BSC-LSA, BSC-LSS and BSC-FRR methods, we considered the event-related design with one hundred 1 s events per condition and mean ISI = 6 s. h The cPPI, BGFC and TSFC matrices were calculated based on the block design simulation. Analogous results were obtained for the event-related design simulation. To evaluate the similarity between matrices, we calculated Pearson’s r correlations between lower diagonal elements. All PPI terms were calculated using the deconvolution step. sPPI and gPPI matrices were symmetrised. The color scales were adjusted for each matrix based on the maximum absolute value and were assured to be positive and negative symmetrical.Despite its simplicity, the CorrDiff approach allows the computation of easily interpretable FC matrices separately for task and rest blocks; these matrices adequately reflect the underlying ground-truth synaptic weight matrices (Fig. 5a, b). The sPPI and gPPI approaches enable estimation of both TMFC and task-independent FC (Fig. 5c, d). Task-independent FC is based on beta coefficients for physiological regressors delivered from selected ROIs and is similar to FC during rest periods. In addition, the gPPI approach can be used to calculate FC separately for each of the task conditions (i.e., Condition > Baseline). However, the gPPI matrices for each task condition should be interpreted with caution. The sign of the PPI estimates between nodes that exhibit high connectivity during rest periods depends on the deconvolution procedure and mean centering of the psychological regressor prior to PPI term calculation. With deconvolution and mean centering, the PPI estimates between these nodes become negative, deviating from the ground truth, see Fig. 5d (for more details, see “Methods”). The BSC approaches enable the calculation of FC matrices for each of the task conditions that are consistent with the ground truth, but do not allow to calculate FC for rest periods (task-independent FC) since rest periods are usually not modelled explicitly (Fig. 5e–g).Co-activations spuriously inflate TMFC estimatesNext, we considered simulations with co-activations to investigate how different TMFC methods address artificial inflation of TMFC estimates due to simultaneous activation of brain regions without task-related modulation of synaptic weights between them (Fig. 6a–c). Sensitivity (true positive rate, TPR) and specificity (true negative rate, TNR) were calculated based on TMFC matrices thresholded at α = 0.001 with false discovery rate (FDR) correction (Supplementary Information 1).Fig. 6: Inflation of TMFC estimates due to co-activations.Simulation results for sample size N = 100, SNR = 0.4 and TR = 2 s. a Expected influence of co-activations on FC estimates (artificial TMFC). To model co-activations, we added simultaneous haemodynamic responses for different functional modules to the simulated BOLD signal without changing the synaptic weights between them. Functional modules №1 and №3 are co-activated in “Cond A”, while modules №2 and №4 are co-activated in “Cond B”. Artificial FC inflation was expected within and between these modules. b Expected FC matrices based on ground-truth synaptic weight matrices (true TMFC). In “Cond A”, synaptic weights were increased between modules №1 and №2 and modules №3 and №4. In “Cond B”, synaptic weights were increased between modules №1 and №4 and modules №2 and №3. c If the TMFC method fails to eliminate co-activations, we will observe FC changes between all ROIs (artificial + true TMFC). d, e Raw and thresholded TMFC matrices obtained using the gPPI method without (w/o) the deconvolution step and with (w/) design matrix upsampling, similar to the gPPI implementation in the FSL or CONN toolbox. TMFC matrices were thresholded at α = 0.001 (two-sided one-sample t test, false discovery rate (FDR) correction). d Without task activation regression, we observed FC changes between almost all ROIs (low specificity). e After finite impulse response (FIR) task regression, artificial TMFC was removed, leaving mostly true TMFC (high specificity). f Sensitivity (true positive rate, TPR) and specificity (true negative rate, TNR) of different TMFC methods for the block design with ten 20 s blocks per condition. g TPR and TNR for the event-related design with one hundred 1 s events per condition and mean ISI = 6 s. All TMFC matrices were calculated with and without FIR task regression, as well as with and without upsampling of the design matrix before convolution, except for CorrDiff and BSC-FRR. We did not perform upsampling for the CorrDiff and BSC-FRR methods. All PPI terms were calculated with and without the deconvolution step. sPPI and gPPI matrices were symmetrised. Boxplots whiskers are drawn within the 1.5 interquartile range (IQR), computed from 1000 random resamplings with replacement. The color scales were adjusted for each matrix based on the maximum absolute value and were assured to be positive and negative symmetrical.As a result, we found that if co-activations are not removed from the fMRI time series before TMFC analysis, they spuriously inflate TMFC estimates in three cases. First, the sPPI method practically does not eliminate co-activations and therefore has near-zero specificity in both block and event-related designs (Fig. 6f–g). Second, the gPPI method demonstrates low specificity when applied to event-related designs without deconvolution (Fig. 6d, g). Third, the specificity of all TMFC methods decreases if the design matrix is not upsampled before convolution with the haemodynamic response function (Fig. 6f, g).Upsampling of the design matrix is used to improve the convolution procedure and is implemented in many popular neuroimaging software packages (SPM, FSL, AFNI, CONN toolbox). However, it may be absent in some in-house TMFC analysis scripts. The most prominent effect of upsampling on specificity can be seen for the gPPI method with deconvolution in event-related designs (Fig. 6g).To better isolate TMFC from co-activation effects, it has been proposed to regress out task activations using finite impulse response (FIR) functions prior to TMFC analysis45. FIR task regression substantially improved the specificity of all TMFC methods in both block and event-related designs (Fig. 6f, g). For instance, the gPPI method without deconvolution had a specificity of 24% before FIR task regression and 99% after (Fig. 6e). The downside of FIR task regression is that it slightly reduces sensitivity, most notably for the sPPI and gPPI methods without deconvolution and the BSC-LSA method.In all subsequent sections, we consider simulations with co-activations and perform TMFC analyses with FIR task regression and design matrix upsampling, unless otherwise stated. We did not perform upsampling for CorrDiff and BSC-FRR, since the CorrDiff method does not rely on general linear models, and FRR implementation in the GLMsingle toolbox requires the task design matrix to have the same temporal resolution as the data to be convolved19. We will not report specificity further since in no case did it fall below 95%.Influence of noise and sample size on the sensitivity of TMFC methodsIn this section, we describe the robustness of TMFC methods to high noise levels and low sample sizes. For the block design, the PPI methods with deconvolution were the most sensitive and robust to noise (Fig. 7a, Supplementary Table S3). The least sensitive method was CorrDiff since it needed a sample size >100 to achieve >80% sensitivity at high SNR = 0.5, while the PPI methods with deconvolution needed a sample size >50 (Fig. 7c).Fig. 7: Sensitivity of all TMFC methods depending on the signal-to-noise ratio (SNR) and sample size.Simulation results for TR = 2 s. a, b Sensitivity of different TMFC methods at sample size N = 100 for (a) the block design with ten 20 s blocks per condition and (b) event-related design with one hundred 1 s events per condition and mean ISI = 6 s. c, d Sensitivity of different TMFC methods depending on the sample size for (c) the block design and (d) event-related design. PPI terms were calculated with (w/) and without (w/o) the deconvolution step. sPPI and gPPI matrices were symmetrised. Boxplots whiskers are drawn within the 1.5 interquartile range (IQR), computed from 1000 random resamplings with replacement.For the event-related design, the BSC-LSS method was the most sensitive (Fig. 7b, Supplementary Table S4). The difference between the BSC-LSS method, the BSC-FRR method and PPI methods with deconvolution was relatively small and more pronounced at high noise levels (SNR < 0.40). At the same time, the BSC-LSA method had the lowest sensitivity due to the multicollinearity problem17. The BSC-LSS and BSC-FRR methods, as well as the PPI methods with deconvolution, needed a sample size N > 50 to achieve >80% sensitivity at SNR = 0.5 (for the event-related design), while the BSC-LSA method needed a sample size N > 100 (Fig. 7d).In contrast to previous simulation studies13,20, we did not detect a noticeable increase in the sensitivity of the gPPI method compared to the sPPI method. Additional Bayesian analysis provided evidence for the absence of differences between these methods for the block and event-related designs (Supplementary Tables S3, S4).It has also been previously suggested that deconvolution may benefit only event-related designs and can be omitted for block designs8. Indeed, some popular neuroimaging packages implement the PPI method without the deconvolution step (e.g., FSL and CONN toolbox). Here, we show that the deconvolution procedure substantially increases the sensitivity of the PPI methods in both block and event-related designs (Fig. 7a, b). Without deconvolution, these methods failed to achieve >80% sensitivity for sample sizes N < 100 and SNR = 0.5 (Fig. 7c, d). Therefore, in the remaining sections, we will consider the sPPI and gPPI methods only with deconvolution, unless otherwise stated.Sensitivity of the TMFC methods for event-related designs with different timing parametersNext, we independently varied the temporal parameters of the event-related design, including event duration, mean ISI, number of events and data acquisition (TR), to assess their impact on the sensitivity of TMFC methods with a sample size N = 100 and medium SNR = 0.40. Shortening the event duration substantially decreased the sensitivity of all TMFC methods (Fig. 8a, Supplementary Table S5). The BSC-LSS method was slightly more sensitive than the BSC-FRR method and the PPI methods, which was more noticeable for short event durations <1 s; in contrast, the BSC-LSA method achieved >80% sensitivity only for long event durations >2 s.Fig. 8: Sensitivity of TMFC methods for event-related designs with different temporal parameters.Simulation results for sample size N = 100 and SNR = 0.4. The default event-related design consisted of one hundred 1 s events per condition and mean ISI = 6 s. a Sensitivity depending on the duration of a single event. b Sensitivity depending on the mean ISI. c Sensitivity depending on the number of events per condition (Ne). d Sensitivity depending on the repetition time (TR) for a fixed scan time (23.6 min). e Sensitivity for the event-related design with short event duration (100 ms), fast data acquisition (≤1000 ms) and fixed number of scans (1238). Fixing the number of scans resulted in variable scan times and numbers of events. The panel on the left side represents raw sensitivity (true positive rate, TPR). The panel on the right side represents normalised sensitivity per single event (TPR/Ne). sPPI and gPPI matrices were symmetrised. Boxplots whiskers are drawn within the 1.5 interquartile range (IQR), computed from 1000 random resamplings with replacement.Shortening the mean ISI slightly decreased the sensitivity of the BSC-LSS, BSC-FRR and PPI methods and substantially reduced the sensitivity of the BSC-LSA method (Fig. 8b, Supplementary Table S6). The BSC-LSA method had >80% sensitivity only for slow event-related designs with mean ISI = 12 s. For the rapid event-related designs (ISI ≤ 4 s), the PPI methods were slightly more sensitive than the BSC-LSS and BSC-FRR approaches. At longer ISI (≥6 s), the BSC-LSS approach was more sensitive than the BSC-FRR and PPI methods.Reducing the number of events also substantially decreased the sensitivity of all TMFC methods (Fig. 8c, Supplementary Table S7). The BSC-LSS approach was slightly more robust to shortening scan time than other methods. The BSC-LSS, BSC-FRR and PPI methods needed at least 80 events per condition to achieve >80% sensitivity.Finally, fast data acquisition (TR < 1 s) yielded the maximum sensitivity for all methods (Fig. 8d, Supplementary Table S8), while reducing the fMRI temporal resolution from the typical TR = 2 s to the frequently used TR = 3 s reduced sensitivity below 80%. Notably, the PPI methods were more sensitive than the BSC-LSS and BSC-FRR methods for TR = 3 s.We also considered TMFC simulations with task designs parameters derived from empirical HCP and CNP tasks. For the working memory task with 8 blocks per condition, block duration = 27.5 s, interleaved by 15 s rest blocks, TR = 0.72 s, total scan time ≈ 10 min (810 dynamics), SNR = 0.4, and sample size N = 100, the sensitivity of the CorrDiff and PPI methods with deconvolution was 61% and 99%, respectively. For the social cognition task with 5 blocks per condition, block duration = 23 s, interleaved by 15 s rest blocks, TR = 0.72 s, total scan time ≈ 7 min (548 dynamics), SNR = 0.4, and sample size N = 100, the sensitivity of the CorrDiff and PPI methods with deconvolution was 23% and 97%, respectively. For the stop-signal task with 96 “Go” and 32 “Stop” events, mean ISI = 1 s (ranged from 0.5 to 4 s), event duration = 1.5 s, TR = 2 s, total scan time ≈ 7 min (184 dynamics), SNR = 0.4, and sample size N = 115, the sensitivity of the BSC-LSS and PPI methods with deconvolution was 3% and 15%, respectively. The sensitivity of BSC-LSS was lower than that of PPI methods due to the very short ISI. When we halved the number of “Stop” events (i.e., only considered “Correct Stop” events), the sensitivity dropped to 0% and 8%, respectively. For the task-switching task with 24 “Switch” and 72 “No Switch” events, mean ISI = 3 s, event duration = 1 s, TR = 2 s, total scan time ≈ 6 min (208 dynamics), SNR = 0.4, and sample size N = 115, the sensitivity of the BSC-LSS and PPI methods with deconvolution was 13% and 15%, respectively.Therefore, block design tasks from the HCP dataset are well suited for TMFC estimation. Although the total duration of these tasks is relatively short, it is compensated by short TR. At the same time, event-related tasks from the CNP dataset may not have sufficient sensitivity to capture whole-brain TMFC. In particular, these tasks are unbalanced and have very few events of interest (32 “Stop” and 24 “Switch” events).Rapid synchronisation can be revealed even with typically slow fMRI data acquisitionOne intriguing question is the principal ability of BOLD fMRI to assess rapid modulations of gamma-band neuronal synchronisation evoked by short events given typically slow data acquisition (TR = 2 s). In the previous section, we showed that the most sensitive TMFC method (BSC-LSS) revealed a fairly small number of true positives (TPR = 15%) for the event-related design with one hundred 100 ms events per condition (Fig. 8a). When we doubled the number of events, the sensitivity was increased to TPR = 57%. Therefore, typically slow fMRI sequences can in principle detect task-related neuronal synchronisation on the order of 100 ms. However, the sensitivity to such rapid synchronisation with slow data acquisition is relatively weak, and many events are required to detect them.Importance of fast fMRI data acquisition for TMFC analysisAnother possible way to increase the sensitivity of TMFC methods to rapid neuronal synchronisation is to employ fast fMRI data acquisition techniques. Increasing the fMRI temporal resolution from TR = 2 s to TR = 500 ms resulted in 99% sensitivity for the event-related design with one hundred 100 ms events per condition when the BSC-LSS method was applied. At a fixed task duration (scan time = 23.6 min), decreasing TR from 2 s to 500 ms increased the number of fMRI data points fourfold (from 616 to 2464 scans). Thus, with the same scanning time, fast data sampling enables one to increase the temporal degree of freedom and thereby significantly increase the statistical power of TMFC methods.Furthermore, fast fMRI data sampling may improve sensitivity per se, that is, by more precise insights into neuronal temporal dynamics rather than simply by providing more data points for a fixed scan time46. To test this assumption, in contrast to our previous simulations with fixed scan time and number of events, we considered a fixed number of data points (scans). As a result, the sensitivity of TMFC methods remained at the same level at different temporal resolutions, TR = 500/700/1000 ms, despite a reduction in total scan time and number of events per condition, Ne = 50/70/100 (Fig. 8e, left panel). This meant an increase in normalised sensitivity per single event (TPR/Ne) at shorter TRs (Fig. 8e, right panel). The BSC-LSS method had the highest sensitivity per event for short event duration and fast fMRI data acquisition. Therefore, despite the sluggishness of haemodynamic processes, a more accurate characterisation of ultra-slow BOLD-signal fluctuations allows for more precise insights into rapid task-related modulation of gamma-band synchronisation.Haemodynamic variability markedly decrease sensitivity of all TMFC methodsIn the previous sections, we considered simulations with fixed parameters of the Balloon-Windkessel haemodynamic model. However, it is well known that the HRF shape varies across brain regions and subjects47. To study this issue, we performed simulations with variable haemodynamic parameters providing a time-to-peak range of 3 to 7 s, consistent with empirical studies48. Here, we compared the best-performing methods (gPPI, BSC-LSS, BSC-FRR), which assume canonical HRF shape in general linear models (GLMs), and the CorrDiff approach, which does not use GLMs. The gPPI and BSC-LSS methods were compared with and without upsampling of the task design matrix prior to convolution with canonical HRF. The upsampling procedure can improve sensitivity for an fMRI signal with a fixed canonical HRF, however its interaction with variable HRF is unknown. Additionally, we compared the gPPI method with and without deconvolution, since this procedure assumes a canonical HRF shape.As a result, we revealed a marked decrease in sensitivity of all TMFC methods (Fig. 9). For the block design with SNR = 0.4, the sensitivity of the gPPI method with task design upsampling and deconvolution dropped from 78% to 5% (Fig. 9a). Disabling upsampling and deconvolution further decreased sensitivity of the gPPI method. As the SNR increases to 0.7, the sensitivity of the gPPI method increased to 71%. Therefore, the gPPI method requires a relatively high SNR to estimate TMFC under conditions of variability in the HRF shape.Fig. 9: Influence of haemodynamic variability on the sensitivity of TMFC methods.Simulation results for sample size N = 100, TR = 2 s and variable parameters of the Balloon-Windkessel haemodynamic model. Sensitivity of different TMFC methods (a) for the block design with ten 20 s blocks per condition and (b) event-related design with one hundred 1 s events per condition and mean ISI = 6 s. The BSC-LSS and gPPI methods were implemented with (w/) and (w/o) upsampling of the design matrix. PPI terms were calculated with (w/) and without (w/o) the deconvolution step. gPPI matrices were symmetrised. Boxplots whiskers are drawn within the 1.5 interquartile range (IQR), computed from 1000 random resamplings with replacement.For the default event-related design with SNR = 0.4, the sensitivity of the BSC-LSS and gPPI methods with upsampling dropped from about 93% to 36% and 14%, respectively (Fig. 9b). At the same time, the BSC-FRR sensitivity decreased from 91% to 44%, making it more robust to HRF variability than the BSC-LSS method. This robustness could potentially be due to fractional ridge regression or the lack of design matrix upsampling in the BSC-FFR method. It turned out that this was due to the second option. In the case of variable HRF simulation, the absence of upsampling of the task design matrix prior to convolution with canonical HRF increased the sensitivity of the BSC-LSS method from 36% to 45% (Fig. 9b). This may be due to the fact that convolution of an upsampled design matrix with a canonical HRF (convolution with high temporal resolution) implies greater model dependence on the canonical HRF shape than convolving the design matrix without upsampling (convolution with low temporal resolution).Thus, the BSC-LSS method without upsampling was the most robust to haemodynamic variability. The BSC-FRR method without upsampling was slightly less sensitive than the BSC-LSS method, but still much robust than the gPPI method. Without upsampling, gPPI sensitivity was increased from 14% to only 18% for SNR = 0.4 (Fig. 9b). The gPPI method was half as sensitive than the BSC-LSS and BSC-FRR methods. Disabling deconvolution, reduced gPPI sensitivity to 3%, even though deconvolution assumes the canonical HRF shape. At higher SNRs, the difference in sensitivity between the gPPI and BSC methods became less noticeable.Genuine and spurious asymmetry of the PPI matricesPreviously, we considered only symmetrised PPI matrices because averaging the upper and lower diagonal elements has become a standard procedure in TMFC analysis since they are considered quite similar8,11. Here, we use empirical and simulated data to determine when PPI matrices become asymmetric, making the averaging procedure problematic. Below, we report results only for the gPPI method since the sPPI and gPPI methods produce nearly identical results.For the block design (working memory task), the correlations between the upper and lower diagonal elements of the group-mean gPPI matrices without and with deconvolution were 0.85 and 0.77, respectively (Supplementary Fig. S7a). Without deconvolution, the correlation coefficients for individual subjects ranged from 0.77 to 0.89 with a mean of 0.83 (Supplementary Fig. S7c). With deconvolution, individual correlation coefficients ranged from 0.56 to 0.75 (mean 0.68). For the event-related design (stop-signal task), the correlations between the upper and lower diagonal elements of the group-mean gPPI matrices without and with deconvolution were 0.90 and 0.40, respectively (Supplementary Fig. S7b). The correlation coefficients for individual subjects ranged from 0.78 to 0.93 with a mean of 0.90 without deconvolution (Supplementary Fig. S7d). In contrast, when deconvolution was applied, individual correlation coefficients ranged from 0.24 to 0.41 (mean 0.34). Similar results were obtained for other tasks with block and event-related designs taken from the HCP and CNP datasets (Supplementary Fig. S8). Therefore, empirical data showed that deconvolution slightly increases the asymmetry of gPPI matrices for block designs and substantially increases the asymmetry for event-related designs.Next, we used simulations with symmetric ground-truth matrices to determine which parameters of the fMRI experiment can artificially increase the asymmetry of the PPI matrices. We found that the main factors for artificial matrix asymmetry are a low SNR (Supplementary Fig. S9a), small sample size (Fig. S9b), short event duration (Fig. S9c) and small number of events per condition (Fig. S9e). In addition, the asymmetry of gPPI matrices increases slightly with larger TRs (Fig. S9f) and is practically independent of the mean ISI duration (Fig. S9d). Therefore, the large asymmetry of the gPPI matrices for the event-related designs from the CNP dataset is most likely related to the short scan time and low SNR.Finally, we considered simulations with asymmetric ground-truth matrices (Supplementary Fig. S10), fixed HRF, and without adding co-activations to test whether the gPPI method could in principle provide information about the true causal directionality. The causal influence that one neuronal system exerts over another at a synaptic or neuronal population level is referred to as effective connectivity (EC). Here, we calculated correlation between the asymmetric ground-truth matrix and the asymmetric gPPI matrix, as well as the ratio between correctly identified sign of connections to the total number of non-zero ground-truth connections (correct sign rate, CSR, see Supplementary Eq. S33). CSR of 50% means the sign was determined by chance. CSR of 100% means that the signs of all connections presented in the ground truth were correctly identified. For the block designs, we also calculated task-modulated EC (TMEC) matrices using the regression dynamic causal modelling (rDCM) method49,50,51 (see “Methods”). The rDCM method, which is a conventional EC method, was used as a reference. As rDCM requires a relatively high SNR49, we used SNR = 5 and twice the total scan duration. If the gPPI method fails to correctly estimate the direction of information flow at a high SNR, then it will also fail at lower SNRs. A systematic comparison of TMEC methods such as Granger causality or structural equation modeling is beyond the scope of the current study.As a result, the gPPI method with deconvolution was able to reflect the actual direction of the information flow for the block design with twenty blocks per condition (Fig. 10). The correlation between the group-mean asymmetric gPPI matrix and ground-truth matrix for the “Cond A – Cond B” difference was 0.93 (CSR = 100%), and correlation between the group-mean rDCM matrix and ground-truth matrix was 0.88 (CSR = 99%). We also ensured that the asymmetry of gPPI regression coefficients was not due to amplitude differences between ROIs during task conditions (see Supplementary Information 9, Figs. S11, S12).Fig. 10: gPPI with deconvolution can reveal the direction of information flow under fixed HRF, high SNR and long scan time.Simulation results for the block design with twenty 20 s blocks per condition, sample size N = 100, TR = 2 s and very high SNR = 5. a Asymmetric ground-truth matrix of task-modulated effective connections. In “Cond A”, synaptic weights were increased from module №1 to №4, from №4 to №3, from №3 to №2, and from №2 to №1. In “Cond B”, synaptic weights were increased in the opposite direction. b Group-mean difference between rDCM matrices calculated for the “Cond A” and “Cond B” blocks. c Group-mean asymmetric matrix generated by gPPI with (w/) deconvolution. Deconvolution allows modelling of psychophysiological interactions at the neuronal level. d Group-mean asymmetric matrix generated by gPPI without (w/o) deconvolution. Without deconvolution psychophysiological interactions, are modelled at the haemodynamic level. To evaluate the similarity between matrices, we used Pearson’s r correlations. The color scales were adjusted for each matrix based on the maximum absolute value and were assured to be positive and negative symmetrical.When the scan duration was halved, the correlation between the gPPI and ground-truth matrices decreased to 0.76 (CSR = 88%), and dropped to 0.12 (CSR = 56%) when SNR was reduced to 0.5. Similar results were obtained for the event-related designs: the correlation between the gPPI and ground-truth matrices was 0.85 (CSR = 99%) for two hundred events per condition, decreased to 0.57 (CSR = 75%) when the scan duration was halved, and dropped to 0.03 (CSR = 51%) when SNR was reduced to 0.5. Without deconvolution, the correlation between the ground-truth and gPPI matrices was always close to zero ( | r | < 0.05, CSR ≈ 50%).Therefore, gPPI without deconvolution failed to estimate the actual direction of the information flow, determined by asymmetric ground-truth synaptic weights, even in the best-case scenario (fixed HRF, high SNR, long scan duration). At the same time, gPPI method with deconvolution was able to reveal the true causal directionality in the best case. However, none of the connections estimated by gPPI survived the FDR-corrected threshold of 0.001 (even though the connection signs were correctly identified). Moreover, when we shortened the scan duration, reduced the SNR, and, most importantly, introduced HRF variability, the ability of the gPPI method to correctly identify the direction of information flow was reduced to almost zero.
