Figure 1c illustrates the general framework of a feedforward ANN (Venkatesan and Anitha, 2006) with one hidden layer. We used a nonlinear regression implemented in Python¹ using ANN. We utilised a tangent hyperbolic activation function as this is widely used (Sharma et al., 2017). The number of units in the hidden layer was set to be the average of the number of nodes in the input and output layers. The network is trained using a standard back-propagation learning algorithm (Rojas, 1996). The solver for weight optimisation was set to “LBFGS” (Limited-memory Broyden–Fletcher–Goldfarb–Shanno), an unconstrained nonlinear optimisation algorithm in the family of quasi-Newton methods. In general, LBFGS has better performance than other available algorithms for small datasets and provides efficient weight optimisation¹. Each pattern was standardised before being input into the ANN by subtracting the mean and dividing by the standard deviation of the pattern.

Equation 1 defines a general mapping between X and Y as follows:

where f ( . ) could be a linear or nonlinear function and E is a zero-mean Gaussian matrix. For nTL-MDPC, we use an ANN to have a more general estimation of f ( . ) than a linear regression allows. According to Equation 2, using the train subset, for any ROIs X and Y, and a set of time points t x and t y , we have:

where X train is the input of the ANN, which is used to predict Y train during training. Below is a step-by-step description of how X train and Y train are related.

Equation 3 defines the input of the hth node in the hidden layer:

where x i is the ith input, w ih is the weight between the ith node in the input and the hth node in the hidden layer, N x is the number of nodes in the input layer, and N H is the number of nodes in the hidden layer.

Equation 4 obtains the output of the hth node in the hidden layer is as follows:

where a h and b h are scaling parameters of the tangent hyperbolic function.

Finally, the output of the jth node in the output layer is obtained by Equation 5 as:

where v hj is the weight between hth node in the hidden and jth node in the output layer, and N y is the number of nodes in the output layer.

After estimating the above parameters and function f ( . ) , the predicted patterns in ROI Y can be obtained by Equation 6 using the test subset as follows:

where Y ̂ test is the predicted pattern in ROI Y at time point t y . Then, Equation 7 computes the explained variance (EV) for each vertex k = 1,…, n y :

Finally, to quantify the multidimensional connectivity metric for each pair of ROIs at each pair of time points, we summarise the above EVs using Equation 8 by averaging across all of the vertices in ROI Y:

To more thoroughly assess the performance in each scenario, we varied the critical parameters, namely the number of trials, vertices, and signal-to-noise ratios (SNRs), within realistic ranges. The number of vertices used was either 5 or 15, i.e., the minimum and maximum number of vertices obtained from the implementation of our clustering approach on our EEG/MEG dataset (Rahimi et al., 2023). To report the final connectivity values and assess the stability of connectivity estimates, each simulation condition was repeated 100 times with different random trial subsamples. Trials were randomly selected without replacement to simulate increasing data availability, and results were averaged across iterations. Standard deviations across runs were used to visualise variability, helping assess the robustness of the connectivity estimates.

To confirm that any apparent advantages of TL-MDPC or nTL-MDPC were not due to an overall increase in their likelihood to identify connections, we assessed whether either showed a propensity to detect false positives. Figure 2a shows this first scenario, representing two independent patterns with no relationship, so that there is no f(.), i.e., Y ≠ f ( X ) . To do so, two pattern matrices with random noise were generated independently using normal distributions (mean = 0 and std. = 1) across 30, 50, 100, 150, and 300 trials.

We then compared TL-MDPC and nTL-MDPC on their ability to detect linear multidimensional relationships over 50, 150, and 300 trials. Figure 2b shows patterns with this relationship, with vertices in ROI X uncorrelated to each other but transformed to ROI Y through a matrix T, so that Y = XT + E , where T is of size n x × n y , and E is the error term. To create these, we produced patterns X through a normal distribution (mean = 0, and std. = 1). For the transformation matrix, we generated matrices using a normal distribution (mean = 0, and std. = 1) with different degrees of sparsity (varying from 10% of the matrix size to 100%, with 10% as the step size). Y was computed by multiplying X and T . Different levels of noise were then introduced from a zero-mean normal distribution with a varying standard deviation (std = 10^std_pow, where std_pow ∈ [−2, −1.5, −1, −0.5, 0, 0.5, 1, 1.5, 2]).

Scenario 3 assessed how well the two methods can identify nonlinear multidimensional relationships, and whether the performance of the nonlinear method was superior. Constructing nonlinear relationships is not as straightforward as making the patterns for the two previous scenarios, as there are infinitely many ways to simulate a nonlinear function. Here, we chose a method that is easily tractable yet inspired by our knowledge of neuronal interactions. We mimicked the process which determines neuronal activation and firing patterns, a linear summation of weighted inputs followed by nonlinear scaling, to create activation patterns that could reflect the kinds of nonlinear relationships found in the brain.

As Figure 2c illustrates, the patterns in ROI X were first linearly transformed through T 0 and the resulting outputs passed through a nonlinear function. We utilised both sigmoid and tangent hyperbolic functions as these are considered similar to how neurons summarise information (Sharma et al., 2017). The resulting matrix then passes through another transformation matrix, T 1 , yielding Y. Specifically, we first generated X and T 0 using a normal distribution (mean = 0, and std. = 1). The output of the first layer is gained through multiplication of X with T 0 ( Y 0 = X T 0 ). Second, this (linear) output is fed through a nonlinear sigmoid or tangent hyperbolic function ( Y 1 = f ( Y 0 ) , where y ij 1 = ( 1 1 + e − y ij o ) , or y ij 1 = tanh ( y ij 0 ) ). The resulting patterns are then transformed through T 1 to produce patterns in ROI Y, so that Y = Y 1 T 1 + E , where E is a zero-mean Gaussian noise of size n t × n y . We introduced zero-mean noise from a normal distribution with a variable standard deviation (std = 10^std_pow, where std_pow ∈ [−2, −1.5, −1, −0.5, 0, 0.5, 1, 1.5, 2]). As in scenario 2, both T 0 and T 1 were generated with different degrees of sparsity, and the same number of trials, vertices, and replications was used.

Our EEG/MEG dataset contains recordings from 18 healthy native English speakers (mean age 27.00 ± 5.13, 12 female) with normal or corrected-to-normal vision. The experiment was approved by the Cambridge Psychology Research Ethics Committee, and volunteers were paid for their time and effort.

The experiment included 250 words and 250 pseudowords. It comprised four blocks presented in a random sequence. One of the four blocks used a lexical decision (LD) task and the other three a semantic decision (SD) task. In the LD block, participants were asked to decide if the presented stimulus was referring to a word or a pseudoword. In SD blocks, they were required to decide whether the presented word was referring to a certain category of words, namely “non-citrus fruits,” “something edible with a distinctive odour,” and “food containing milk, flour or egg.” A stimuli of 10% belonged to these target categories and required a button-press response. As in our previous studies, we only analysed brain responses to real, non-target words. Each stimulus was presented for 150 ms, with an average Stimulus Onset Asynchrony of 2,400 ms. A summary of the stimuli is presented in Supplementary Figure S1.

MEG/EEG recordings were collected simultaneously using a Neuromag Vectorview system (Elekta AB, Stockholm, Sweden) and MEG-compatible EEG cap (EasyCap GmbH, Herrsching, Germany) at the MRC Cognition and Brain Sciences Unit, University of Cambridge, United Kingdom. MEG was acquired via a 306-channel system consisting of 204 planar gradiometers and 102 magnetometers. EEG was collected via a 70-electrode system with an extended 10–10% electrode layout. The data sampling rate was 1,000 Hz.

As we intended to directly compare the nTL-MDPC results to those previously published using TL-MDPC, the preprocessing and source estimation steps are identical to Rahimi et al. (2023). We used MEGIN Maxwell-Filter software to apply signal source separation (SSS) with its spatiotemporal extension to remove noise from spatially distant sources and to compensate for small head movements (Taulu and Kajola, 2005). We used the MNE-Python software package (Gramfort et al., 2013, 2014) to perform the preprocessing and source reconstruction. The raw data from each participant was visually checked, and bad EEG channels were selected for interpolation (max = 9 channels per person, min = 0, mean = 2.85). A finite-impulse-response (FIR) bandpass filter between 0.1 and 45 Hz was applied, as well as the FastICA algorithm to remove eye movement artefacts (Hyvarinen, 1999; Hyvärinen and Oja, 2000). Afterwards, epochs were created using the data from 300 ms pre-stimulus to 600 ms post-stimulus. A flowchart of the different steps of pre-processing is shown in Supplementary Figure S2.

L2-Minimum Norm Estimation (MNE) (Hämäläinen and Ilmoniemi, 1994; Hauk, 2004) was used for source estimation. We assembled inverse operators based on a 3-layer Boundary Element Model (BEM) of the head geometry, gained from structural MRI images. For this purpose, we assumed sources to be perpendicular to the cortical surface (“fixed” orientation constraint). We obtained the noise covariance matrices using 300 ms-baseline periods and then selected the best choice from a group of methods included in MNE-Python (‘shrunk’, ‘diagonal_fixed’, ‘empirical’, and ‘factor_analysis’) (Engemann and Gramfort, 2015). MNE-Python’s default SNR = 3.0 was used to regularise the inverse operator for evoked responses. The source signals from each participant were then morphed to the standard Freesurfer brain (fsaverage).

As in Rahimi et al. (2023), we used six regions of interest including left and right anterior temporal lobe (ATL), left inferior frontal gyrus (IFG), left posterior temporal cortex (PTC), left angular gyrus (AG) and left primary visual area (PVA) to study connectivity within the semantic network. These ROIs were constructed using the Human Connectome Project (HCP) parcellation (Glasser et al., 2016). An analysis of leakage among these ROIs can be found in our previous paper (Rahimi et al., 2023). To assess and quantify potential spatial leakage among our six ROIs, we previously computed a leakage matrix based on point spread functions (PSFs) derived from the resolution matrix of our source estimates (Rahimi et al., 2022). This matrix indicated the degree to which signal from one ROI leaked into others, compared to self-leakage. The results showed that all ROI pairs exhibited low to medium leakage, except for slightly higher leakage from lATL to PTC. Importantly, this pattern did not correspond to the strongest observed connectivity, suggesting that our multidimensional connectivity results are unlikely to be driven primarily by leakage artefacts.

We performed nTL-MDPC to predict X from Y and Y from X at every 25 ms from 100 ms pre-stimulus to 500 ms post-stimulus and reported the EVs (averaged across the two directions) as the final metric. EVs for every pair of ROIs and across latencies are presented in TTMs (Rahimi et al., 2023). In any TTM, each row indicates statistical relationships between Y at a specific time point and X at each time point, while columns show statistical relationships between X at a certain time point and Y across all time points, enabling us to display connectivity at different time lags. A flowchart illustrating the connectivity estimation process is presented in Supplementary Figure S3.

To avoid any potential bias due to the different number of trials between SD and LD tasks (Bastos and Schoffelen, 2016), TTMs were computed independently for each of the three SD blocks and then averaged, allowing a better comparison with LD. SD and LD connectivity were then contrasted using non-parametric cluster-based permutation tests, implemented in MNE Python (Maris and Oostenveld, 2007). This approach was chosen because it does not assume normality, making it suitable for high-dimensional EEG/MEG data, and it provides a built-in correction for multiple comparisons (Candia-Rivera and Valenza, 2022) by evaluating the significance of clusters as a whole, rather than at individual time-lag coordinates. We implemented two-sided t-tests thresholded at a t-value corresponding to an alpha-level of 0.05 and 5,000 randomised repetitions. To avoid spurious or small clusters, we only presented clusters whose size exceeded 2% of the TTMs total size, i.e., with more than 12 elements.

The TTM results derived from both the TL-MDPC and nTL-MDPC methods were analysed to determine whether there was a significant difference in the effects captured by each approach. Specifically, we aimed to assess whether the nonlinear nTL-MDPC approach demonstrated greater connectivity—either in terms of magnitude or across additional lags—compared to the TL-MDPC method. This comparison was conducted using cluster-based permutation tests with the same parameters and statistical configurations as in the previous paragraph, including two-sided t-tests thresholded at a t-value corresponding to an alpha-level of 0.05 with 5,000 randomised repetitions. Notably, this comparison was performed at the cluster level rather than the individual cell level, focusing exclusively on clusters with sizes greater than 12 elements (i.e., clusters larger than 2% of the TTM’s total elements) to avoid spurious or small clusters. Through these tests, we evaluated whether the nonlinear adjustments in the nTL-MDPC method revealed any additional or distinct significant effects between the SD and LD conditions that were not evident when using the linear TL-MDPC method.

Figure 3 demonstrates connectivity scores (EV, y-axis) for the case where there is no true relationship between two patterns, with the x-axis showing different numbers of trials and different curves representing different numbers of vertices for both panels. For both MDPC (panel a) and nMDPC (panel b), all values are close to zero with highly overlapping error bars, suggesting that neither method is prone to false positive errors when no relationship exists. The fluctuations in explained variance reflect expected random variability across runs due to trial-to-trial differences, the stochastic nature of the training–testing framework, and the overall small magnitude of EV. These variations do not indicate an inherent sensitivity of MDPC or nMDPC to trial numbers but rather reflect noise across repeated simulations.

Figure 4 demonstrates how well the MDPC and nMDPC methods perform in the case where multidimensional patterns are linearly related. Red curves show the result of the MDPC method, and blue curves represent the nMDPC results. For both methods and in all cases, EV reaches values between 0.8 and 1 for SNRs above 25 db and is close to zero for SNRs smaller than −10 db. For intermediate SNRs, the two methods perform similarly with highly overlapping error bars, with EV rising for increasing numbers of trials. Thus, both linear and nonlinear MDPC similarly capture true linear relationships between patterns to a similar degree, with MDPC producing slightly greater EV values.

Figure 5 shows the connectivity (EV, y-axis) between pairs of multidimensional patterns with a nonlinear relationship, for different SNRs (x-axis), trials (indicated by line type), and numbers of vertices (three panels), generated using sigmoid (left panels) and tangent hyperbolic (right panels) activation functions. EV is around zero for all cases at very low SNRs (<−10 dB). Unlike for linear dependencies, the EV does not approach 1 at high SNRs. Nonlinear MDPC generally outperforms its linear counterpart, with a numerical difference of approximately 0.1 and with little overlap of error bars. Figures 5e,f show that the worst performance for both methods is obtained with the largest number of vertices and the lowest number of trials. Interestingly, in this case, the linear method outperforms the nonlinear one. This is likely due to the known requirement of ANNs for a lot of training data (Hagan and Demuth, 1999). Thus, linear MDPC provides a good approximation of the nonlinear multidimensional relationships, particularly when a large number of parameters needs to be estimated from a small number of trials. Yet, in general, nMDPC captures some more variance. This improved ability of nMDPC to detect additional nonlinear relationships could allow the identification of additional functional connections in brain data, including connectivity between more regions or across additional time lags.

We computed the functional connectivity of all possible pairs of the six semantic ROIs using nTL-MDPC. Figure 7 shows the inter-regional connectivity matrix (ICM) (Rahimi et al., 2023), a summary of the semantic decision versus lexical decision TTMs for all ROI comparisons, with the upper diagonal (green area) showing the nTL-MDPC results, and the lower diagonal (blue area) displaying the TL-MDPC results. Generally, nTL-MDPC identified greater connectivity with deeper semantic processing between the same pairs of regions as TL-MDPC (in lATL-rATL, lATL-PTC, lATL-IFG, rATL-PTC, rATL-IFG, rATL-AG, rATL-PVA, PTC-IFG, PTC-AG, PTC-PVA, and AG-PVA). As found with TL-MDPC, the connectivity between core semantic regions (including lATL-rATL, lATL-PTC, lATL-IFG, rATL-PTC, and PTC-IFG), as well as AG-PVA, appears to “fan out” from the diagonal over time, i.e., connectivity persists over longer time lags at later latencies on both sides of the diagonal, possibly reflecting bidirectional or recurrent information flow (McClelland and Rumelhart, 1989; Rogers and McClelland, 2014; Clarke et al., 2018; Kietzmann et al., 2019). Connectivity was also identified over short time lags for rATL-IFG, and late effects were revealed for rATL-AG and rATL-PVA.

The results of nTL-MDPC were overall similar to those of TL-MDPC, and indeed, no additional ROI pairs were found to be connected. However, nTL-MDPC captured more connectivity in some cases. Specifically, for some pairs, it revealed significant modulations across a greater number of time lags, particularly those over greater time windows, e.g., between lATL-rATL, lATL-PTC, rATL-PTC, and PTC-IFG. However, these differences were subtle, and care should be taken when interpreting cluster size with cluster-based thresholding. Moreover, when we formally contrasted the two methods directly using cluster-based permutation tests, there were no significant differences at the cluster level (note that cell-wise comparisons were not conducted).

As multiple different methods have now been used to assess connectivity within the same dataset, we summarised the networks extracted using these different methods in Figure 8. This demonstrates the results from: (1) coherence in four frequency bands (Rahimi et al., 2022), (2) TL-UDC (Rahimi et al., 2023), (3) TL-MDPC (Rahimi et al., 2023), and (4) nTL-MDPC analyses, for early (50–250 ms) and late (250–500 ms) time windows. TL-UDC is the unidimensional counterpart of TL-MDPC (Rahimi et al., 2023), in which all the time courses of the vertices within each ROI are collapsed into one, allowing a direct assessment of the effects of multidimensionality. As coherence was computed per frequency band, we represent the connections specific to the gamma band in blue (i.e., rATL-PTC, IFG-AG at the first time window, and PTC-AG at the second time window), and connections consistent across three frequency bands (i.e., alpha, beta, and gamma) in yellow. To summarise across time lags for TL-UDC, TL-MDPC, and nTL-MDPC, we summed the significant t-values in each time window as a metric of connection intensity. This reflects increases in connectivity strength and time lags with significant connectivity across longer time lags, indicating stronger or more reliable connectivity across time and participants. We represented these values with the arrow widths (with the weakest connectivity reflected by the thinnest lines, as in rATL-IFG with nTL-MDPC at the earlier time window, and the strongest connectivity represented by the thickest line, as in AG-PVA with nTL-MDPC at the later time window). Note that, unlike coherence-based connectivity metrics—which assess phase synchronisation between signals, typically within narrow frequency bands—TL-MDPC and nTL-MDPC estimate the extent to which the multidimensional pattern of activation in one ROI at time t₁ can predict the pattern in another ROI at time t₂. These methods quantify connectivity as the proportion of explained variance (EV) in the target pattern, based on a trained transformation function (e.g., an artificial neural network in the nonlinear case). As such, they capture statistical dependencies—both linear and nonlinear-across multivariate spatiotemporal activation patterns, rather than oscillatory coupling. This distinction is critical for understanding what aspects of brain dynamics each method operationalises.

In general, the two unidimensional connectivity methods revealed fewer and weaker task modulations of connectivity, whereas both multidimensional methods detected richer connectivity across the whole semantic network across the whole latency range. Both multidimensional methods identify rich interactions between the four core semantic representation and control regions (lATL, rATL, PTC, and IFG) starting at early time points (around the onset of the stimulus) and becoming stronger across the time window. Even with methods designed to identify nonlinear relationships, the connectivity of AG is relatively sparse. It only shows connectivity with PVA and PTC from earlier time points, and with rATL at later latencies. The general pattern of results was similar for TL-MDPC and nTL-MDPC, with both distinguishing a core semantic network comprising lATL, rATL, PTC, and IFG from a posterior visual attention network consisting of AG, PVA, and PTC, with overlap of these two sub-networks in PTC.

While the contrast between our linear and nonlinear TL-MDPC TTMs did not yield any significant effects, visual inspection of Figure 8 reveals some subtle numerical differences between the two approaches. Using nTL-MDPC, we found that the connectivity values were numerically higher for lATL-rATL and PTC-IFG at early time points, and for lATL-PTC, PTC-IFG, and rATL-PVA at the later time points. However, connectivity values were numerically higher for lATL-IFG and PTC-AG when using the linear TL-MDPC method at the late time window. The finding of a greater number of numerically stronger connections for the nonlinear method could suggest some subtle benefits for this method in the overall reliability of the evidence of connectivity between the regions of the semantic network for the two methods. However, these differences are very limited, and they did not result in any significant differences in the statistical comparison between the TTMs.