Non-reversibility outperforms functional connectivity in characterisation of brain states in MEG data

Characterising brain states during tasks is common practice for many neuroscientiﬁc experiments using electrophysiological modalities such as electroencephalog- raphy (EEG) and magnetoencephalography (MEG). Brain states are often described in terms of oscillatory power and correlated brain activity, i.e


Introduction
Characterisation and identification of task induced brain states is a common and widely applied practice in the field of functional neuroimaging as complex cognition from the human brain presumably emerges from the orchestration and evolution of a repertoire of brain states. Task induced brain states can loosely be defined as any transient macroscopic configuration of the brain characterised by a descriptive statistic of choice. Brain states estimated from functional imaging modalities such as functional MRI (fMRI), electroencephalography (EEG), and magnetoencephalography (MEG) are traditionally described in terms of correlated brain activity, i.e. pairwise functional connectivity Abbreviations: MEG, magnetoencephalography; AEC, amplitude envelope correlation; LAEC, lagged amplitude envelope correlation; GEC, generative effective connectivity. erates far away from thermodynamic equilibrium ( Battle et al., 2016 ). In thermodynamic equilibrium, there is no entropy production as the system has reached its state with maximum entropy. Whenever a system in equilibrium visits a temporal sequence of states, the probability of this sequence of states is equally likely to the probability of visiting the reverse of this sequence of states. An operating point far away from thermodynamic equilibrium is extremely relevant for the brain. In case this non-equilibrium principle is violated, the brain's properties would become stationary in time ( Battle et al., 2016;Lynn et al., 2021;Perl et al., 2021 ). Non-equilibrium dynamics can be captured by assessing the asymmetry in the temporal sequences of states, i.e. non-reversibility, giving rise to the notion of the "arrow of time " and production of entropy. The concept of non-equilibrium dynamics is not only interesting for its own sake, but has potential explanatory power to relate to various other dynamical properties of the brain, such as time-varying connectivity and turbulence, and may form a basis as to why the brain is functioning at different hierarchical levels ( Deco and Kringelbach, 2020;Deco et al., 2021b;Escrichs et al., 2022 ).
Recent fMRI work has demonstrated that the concept of nonreversibility is superior to conventional functional connectivity in the identification of brain states ( Deco et al., 2022a ). It was demonstrated that whenever the brain is more engaged in processing information from the environment during several cognitive tasks, this leads to an increase in non-reversibility and entropy production compared to the restingstate ( Deco et al., 2022b ). A whole-brain computational model was further used to gain insight into mechanisms underlying non-reversibility and entropy production, which showed that an important driving factor for non-zero non-reversibility is asymmetry in the effective connectivity that neuronal populations perceive from one another . Another study using electrocorticography (ECoG) recordings demonstrated that while functional connectivity showed striking resemblance between several conscious states (e.g. during the awake, ketamine and recovery phase in monkeys), the signature of different conscious states could be better differentiated using measures of nonequilibrium dynamics ( Deco et al., 2022a ).
In order to gain more understanding of the relevance of the "arrow of time " in human data, the next step is to analyse non-reversibility in signals that are characterised by non-stationarity. Therefore, in our current work, we translate concepts of non-equilibrium dynamics to MEG data. We assess whether non-reversibility outperforms conventional functional connectivity in identification of task condition contrasted to resting-state. In this context, we study a sensorimotor task, a working memory task and a language task. Since there is evidence that amplitude coupling provides a reliable estimate for functional connectivity ( Colclough et al., 2016 ), we restrict our functional connectivity analysis to amplitude coupling, i.e. the amplitude envelope correlation (AEC), ( Brookes et al., 2011a;Hipp et al., 2012 ). Non-reversibility in MEG data is assessed using the lagged version of the amplitude envelope correlation (LAEC) ( Basti et al., 2019 ). As a second step, we analyse using surrogate data whether the observed non-reversibility in empirical data is a genuine sign of temporal asymmetry. We test our empirical data against surrogate data that possess a symmetric cross-covariance and hence corresponds to the null-hypothesis of reversibility. Lastly, we investigate the contributing factors for the emergence of non-reversibility in MEG data using neural mass modelling, with two potential candidates, asymmetry in the effective connectivity and heterogenous axonal conduction delays.

Diffusion MRI: estimation of structural networks
The pipeline of the structural network construction has been described in Tewarie et al. (2019aTewarie et al. ( , 2022 . We applied the same method with a slight adjustment mentioned below. We included diffusion MRI data from ten healthy controls (who also underwent MEG recordings) of the Human Connectome project ( Larson-Prior et al., 2013 ). Diffusion MRI data were obtained from the Human Connectome Project ( Van Essen et al., 2013 ). Full acquisition protocol details are described in Sotiropoulos et al. (2013) . Briefly, a monopolar Stejskal-Tanner echo planar imaging sequence was used in a 3T Siemens Connectom Skyra to acquire data at (1.25 mm) 3 isotropic resolution. Diffusion-sensitization was applied with three b-values (b = 1000, 2000 and 3000 s/mm2) and along 90 directions per b-shell. Two repeats were obtained with blip-reversed phase encoding. The minimally processed data were used , where susceptibility-induced distortions, eddy currents and subject motion were all corrected simultaneously using a non-parametric framework ( Andersson and Sotiropoulos, 2016 ) based on Gaussian processes ( Andersson and Sotiropoulos, 2015 ). Fibre orientations were estimated using a parametric spherical deconvolution model and were fed into probabilistic tractography in FSL to estimate structural networks ( Behrens et al., 2007;Hernandez-Fernandez et al., 2019 ). In contrast to previous work ( Tewarie et al., 2019a;, streamlines were seeded from 60,000 standard-space vertices in the white matter (5000 streamlines per seed). Connectivity was quantified as the number of streamlines reaching each vertex normalised by the total number of valid streamlines propagated. Using the automated anatomical labelling (AAL) parcellation ( Tzourio-Mazoyer et al., 2002 ), this connectivity was reduced to a 78 × 78 parcellated connectome, by computing for each pair of regions the mean structural connectivity between all pairs of vertices that they were comprised of.

MEG: Data acquistion and pre-processing
Resting-state and task-based MEG data were obtained from the Human Connectome Project ( Van Essen et al., 2013 ) as part of the HCP MEG2 release. Briefly, data were collected on a whole-head Magnes 3600 scanner (4D Neuroimaging, San Diego, CA, USA) from 89 subjects ( Larson-Prior et al., 2013;Van Essen et al., 2013 ); 95 subjects were included in the release, but resting-state recordings that passed the quality control checks (which included tests for excessive SQUID jumps, sensible power spectra and correlations between sensors, and for sufficiently many well-behaved recording channels) were not available from six. All subjects were young (22-35 years of age) and healthy. Resting state measurements were taken in three consecutive sessions for each subject with little or no break in between, for 6 min each. The data have been provided pre-processed ( Larson-Prior et al., 2013 ), after passing through a pipeline to remove any artefactual segments of time from the recordings, identify any recording channels that are faulty, and to regress out artefacts which appear as independent components in an ICA decomposition with clear artefactual temporal signatures (such as eye-blinks or cardiac interference). Task based data was collected in the same way as resting-state data. We included a motor task (54 subjects), language task (77 subjects) and working memory task (76 subjects).

MEG: Task data
A detailed description of the tasks can be found in Larson-Prior et al. (2013) .
1. Motor task. In the motor task participants are instructed to make simple hand or foot movements after a visual cue. This task is divided into task and rest blocks, with more task blocks than rest blocks. During the task block a participant is instructed to either make a movement of the right/left foot or hand. This task activates regions in the sensorimotor network in the alpha (mu), beta and gamma band ( Crone et al., 1998;Pfurtscheller and Da Silva, 1999 ).
2. Language task. In this task participants listen to auditory narratives (30 s duration) or matched-duration simple arithmetic problems. This is followed by a 2-alternative forced choice question to which participants respond by a right hand button press. Previous MEG studies show that the language network and regions adjacent to this language network area activate during this task. Renvall et al. (2012) , Pulvermüller (2010) . Speech modulation is especially encoded in the lower frequency bands (theta and alpha band) ( Ding and Simon, 2013 ). 3. Working memory task . A N-back task was performed during the recording. Tools or faces are presented to participants in an alternating 0-back or 2-back fashion. Participants were instructed to press on a button with their right index or right middle finger for matched or non-matched responses respectively. This task tests the ability of perception and long term memory ( Baddeley, 2003 ). Electrophysiological responses are expected in the theta and alpha band and in prefrontal and parietal cortical areas ( Brookes et al., 2011b;Collette et al., 2006;Jensen et al., 2002;Klimesch, 2006 ).

MEG: Source localisation
A description of the source localisation of this dataset is provided in Tewarie et al. (2019a) . An atlas-based beamforming approach was adopted to project MEG sensor level data into source-space ( Hillebrand et al., 2012 ). The cortex was parcellated into 78 cortical regions according to the AAL atlas (same as for structural network). This was done by registering each subject's anatomical MR image to an MNI template and labelling all cortical voxels according to the 78 cortical regions of interest ( Gong et al., 2009 ). Subsequently, an inverse registration to anatomical subject space was performed and the centroid voxel for every region of interest was extracted to serve as representative voxel for every region ( Hillebrand et al., 2016 ). Pre-computed single-shell source models are provided by the HCP at multiple resolutions ( Nolte, 2003 ), registered into the standard co-ordinate space of the Montreal Neuroimaging Institute. Data were beamformed with depth normalisation onto centroid voxels using normalised lead fields and estimates of the data covariance. Covariance was computed for broadband data (1-45 Hz) with a time window spanning the whole experiment ( Brookes et al., 2008 ). Regularisation was applied to the data covariance matrix using the Tikhonov method with a regularisation parameter equal to 5 % of the maximum eigenvalue of the unregularised covariance matrix. Dipole orientations were determined using singular value decomposition to select the source orientation that maximises the output signal-to-noise ratio ( Sekihara et al., 2004 ). This complete process resulted in 78 electrophysiological timecourses, each representative of a separate AAL region.

Functional connectivity and non-reversibility in MEG data
Functional connectivity was estimated using the amplitude envelope correlation metric (AEC) ( Brookes et al., 2011a;Hipp et al., 2012 ). Source reconstructed data were frequency filtered into five frequency bands: delta (1-4 Hz), theta (4-8 Hz), alpha (8-13 Hz), beta (13-30 Hz) and gamma (30-48 Hz). This was followed by pairwise orthogonalisation to reduce the effect of signal leakage ( Hipp et al., 2012 ). The amplitude envelope ( ) for every timecourse was subsequently extracted from from these leakage-reduced frequency-filtered timecourses by calculating the absolute value of their analytical signals. The AEC was estimated by computing the Pearson correlation between pairwise amplitude envelopes. The AEC was computed for a window spanning the whole experiment and was estimated between all possible pairs of timecourses forming an (AEC) functional connectivity matrix. and values were averaged to obtain a symmetric functional connectivity matrix.
We capture non-reversibility (i.e. the arrow of time) through the degree of asymmetry obtained by comparing pairwise time series of the forward and the artificially generated reversed backward version of the amplitude envelopes ( ) . Let us consider two amplitude envelopes from two separate brain regions ( ) and ( ) . By flipping ( ) , denoted as ( ) = (− ) , we obtain the reversed backward version ( ) . Now we can estimate the time-lagged cross-correlation of the forward and backward evolution of the amplitude envelopes (LAEC) (1) , , (Δ ) = corr ( ( ) , ( + Δ )) . (2) In order to work with positive values, we use the expression of mutual information for Gaussian variables ( Baker, 1970 ) to transform the expression for the LAEC into We used the abbreviation FS, which stands for functional causal dependency to keep our notation consistent with our recent fMRI work ( Deco et al., 2022a;G-Guzmán et al., 2023 ). We computed the LAEC using a window spanning the whole experiment. However, note that for (relatively) short windows amplitude envelope data are usually not Gaussian. Therefore, in that case, a Fisher transformation should be applied to the correlations values obtained from equations 1 and 2 before transforming these values to mutual information (equations 4 and 5). The extent of non-reversibility is obtained by capturing the asymmetry between , , (Δ ) and , , (Δ )) . This is expressed as the quadratic distance between the forward and reversal time-shifted matrices The notation ‖ ‖ 2 is defined as the mean value of the absolute squares of the elements of the matrix P. Note that we can also obtain nonreversibility for each brain region by evaluating Equation 6 for every row of the difference matrix separately. Lastly, Δ that results in the highest is chosen for further analysis.

Insensitivity to field spread
Due to residual mixing, the reconstructed MEG source signals are instantaneous linear mixtures of the true source signals. This phenomenon is known as field spread ( Schoffelen and Gross, 2009 ). To avoid false positive connections, connectivity measures that are insensitive to field spread are desirable. A connectivity measure is insensitive to field spread if the absence of a connection between two signals implies the absence of a connection between the observed (i.e. mixed) signals. We show that reversibility is insensitive to field spread. Thus, non-reversibility between the observed signals cannot be explained by field spread and hence reflects non-reversibility of the true signals.
Let ( ) and ( ) be the analytic signals from brain regions and and suppose that they are reversible, that is, their covariance function is symmetric: , (Δ ) = , (−Δ ) . In the above definition of the covariance function, the brackets denote averaging over time. Let be a matrix that models the mixing of two signals due to field spread. Note that because field spread is instantaneous, the entries of the matrix are real-valued. The observed signals ′ ( ) and ′ ( ) are related to the true signals by [ ′ ( ) .
Using the fact that , (Δ ) = , (−Δ ) , the covariance function between the observed signals can be written as Since the functions , and , are symmetric, it follows that the covariance function between the observed signals is symmetric as well: . This shows that if the true signals are reversible, the observed signals are reversible as well and it implies that observed nonreversibility between two signals cannot be explained by signal leakage, but reflects non-reversibility between the true signals.
The situation is a bit more complicated when working with amplitude envelopes instead of the signals proper. In particular, reversibility of the true amplitude envelopes does not necessarily imply reversibility of the observed amplitude envelopes. However, a sufficient condition for reversibility of the observed amplitude envelopes is that the multivariate process is reversible. The first and fourth entries of ( ⊗ )( ) are the amplitude envelopes of the th and th brain region, respectively, and the second and third entries are the "cross-amplitude envelope " ( ) ( ) * and its complex conjugate, respectively. Reversibility of the process ⊗ means that the covariance function between any pair of its entries is symmetric. Since the observed process is related to the true process via where ⊗ denotes the Kronecker product of with itself, it follows that reversibility of ⊗ implies reversibility of ′ ⊗ ′ . In particular, non-reversibility of the observed amplitude envelopes implies nonreversibility of at least one pair of entries in ⊗ .

Machine learning classification and statistics
We employed a random forest algorithm to classify task data from resting-state data. We used the in-built implementation in MATLAB, the so-called 'TreeBagger' function, which is based on Breiman's random forests ( Breiman, 2001 ). Classification was performed for functional connectivity and non-reversibility separately. We used five features as input to the random forest classifier, which included mean AEC for every frequency band or non-reversibility per frequency band. Input data were divided into a training (80%) and test set (20%). Based on the outof-bag error we set the number of trees to 100 for every classification (one classification for every of the three tasks). We report the area under the curve (AUC) of the receiver operating characteristic curve of the classification obtained from the test set.
Testing for significant difference between connectivity or nonreversibility distributions was assessed using the Wilcoxon rank sum test. Correction for multiple tests was performed using the false discovery rate ( Benjamini and Hochberg, 1995 ).

Construction of null-data
We construct surrogate data in order to test the null-hypothesis of reversibility in MEG data. We follow the method described in Hindriks et al. (2018) . Let us consider two time-series = [ 1 , 2 , … , ] and = [ 1 , 2 , … , ] , which are transformed into the Fourier domain to and . We now adjust the Fourier coefficients by taking their real parts and multiplying these by phases that are uniformly distributed in the interval [0 , 2 ] . Thus, the adjusted Fourier coefficients ̂ and ̂ are given by and where is random on [0 , 2 ] and independent for different . Transforming back to the time-domain yields surrogate signals ̂ and ̂ . Using the same phase for both and ensures that the auto-and cross-correlation functions of and are retained ( Schreiber and Schmitz, 2000 ). Taking the real parts of the Fourier coefficients ensures that the auto-covariance function between ̂ and ̂ is symmetric, corresponding to the null hypothesis of reversibility.

Whole-brain computational models for non-reversibility
We constructed a whole-brain model to reveal causal mechanisms of non-reversibility in MEG data. To this end, we fit the empirical AEC and non-reversibility by creating the generative effective connectivity (GEC) and by introducing axonal conduction delays. Whole-brain models have three main constituents: the structural connectivity, the coupling function and the local model. 1) We use the average structural connectivity matrix across subjects as described in section "Diffusion MRI: estimation of structural connectomes ", connecting 78 cortical brain regions. 2) We use a standard additive coupling function ( Pietras and Daffertshofer, 2019 ).
3) The Wilson-Cowan model is used as local model and mimic for MEG data ( Wilson and Cowan, 1972 ). This mode has been widely used for modelling electrophysiological brain activity ( Daffertshofer et al., 2018;Deco et al., 2008;Izhikevich, 2007 ).
Our local model consists of two distinct neuronal populations, an excitatory and an inhibitory neuronal population. The dynamics of a local excitatory and inhibitory population are characterised in terms of their mean firing rates ( ( ) = excitatory, ( ) = inhibitory), which evolve due to local interactions between the excitatory and inhibitory units within the populations, as a consequence of some unaccounted external input , and due to excitatory influence from connected nodes through additive coupling. The sum of all inputs is converted using a sigmoid function ( ) = (1 + − ) −1 , with threshold . The dynamics of a system of Wilson-Cowan oscillators with excitatory and inhibitory populations with additive coupling is described by Parameters , with ∈ { , } and ∈ { , } , refer to coupling strength between local populations, corresponds to the generative effective connectivity between regions and rather than the structural connectivity. The generative effective connectivity is the effective weighting of the structural connectivity (see next paragraph for an explanation). Parameter (in −1 ) refers to a relaxation time constant which is assumed to be equal between excitatory and inhibitory populations. The incoming firing rates from distant excitatory populations are tuned by the global coupling strength parameter and incoming firing rates are delayed by a Euclidean distance dependent delay 0 . Time series of ( ) were used as mimic for MEG signals. The external input is tuned such that the working point of the model is just before a Hopf-bifurcation in the linear regime = 4 (see Tewarie et al., 2019b for a bifurcation diagram). Implementation of the model and model parameters are exactly the same as in Tewarie et al. (2020) and differential equations were numerically solved using a 4th-order Runge-Kutta scheme with a sufficiently small time step ( 1 × 10 −4 s) ( Lemaréchal et al., 2018 ).
We optimised generative effective connectivity between brain areas by comparing the output of the model with the empirical measures of forward and reversed cross-correlations of the amplitude envelopes as well as the empirical AEC. Using a heuristic gradient algorithm, we proceed to update the generative effective connectivity such that the fit is optimised: ] .

(11)
Here and , correspond to a AEC transformed mutual information measure Equations 9,10 and 11 are solved recursively until the fit converges to a stable value. Note that for optimization, we also used the forward and reversed cross-correlations of the model simulated amplitude envelopes as well as the model simulated AEC. The generative effective connectivity is initialised using the structural connectivity and the update of G is only restricted to existing connections of the structural connectivity matrix. The only exception are homologuous connections between mirrored regions in each hemisphere given the a-priori information that tractography is less sensitive to identify these connections. We set = 0 . 05 and ′ = 0 . 01 .

Functional connectivity or non-reversibility based classification in MEG
We quantified the AEC and non-reversibility derived from the LAEC in task and resting-state MEG data. As an illustration, Fig. 1 A shows for all subjects separately, the non-reversibility averaged across brain regions in different frequency bands as a function of the lag or delay for the motor task. For all subjects, we observe clear peaks for the nonreversibility, with for some subjects a second peak corresponding to a local or global maximum. The data shows a clear frequency dependency of the delay corresponding to the maximum non-reversibility. This delay is relatively long for the delta, theta and alpha band compared to the beta and gamma band, with lags of 200 ms for delta, theta, alpha and lags of 30-40 ms for beta and gamma bands. For subsequent analysis, we selected the non-reversibility for the the lag that corresponded to the first maximum of non-reversibility for a subject. This first maximum is for most subjects also the global maximum. Note from Fig. 1 A that there is limited variability in this lag (or delay) between subjects.
1. Motor task. Figure 1 B shows whole-brain non-reversibility values (first maximum of non-reversibility) for every condition and frequency band along with the whole-brain average AEC for every condition and frequency band. For the motor task we see strong task induced effects for non-reversibility in all frequency bands, with a very strong effect in the gamma band. The direction of change of non-reversibility for the motor task was as expected, with an increase in gamma and beta nonreversibility. For whole-brain functional connectivity, we observed significant effects in fewer frequency bands. For all frequency bands that showed significant differences between resting-state and motor task condition, we observed the same change of direction for non-reversibility and functional connectivity. Also note that a global increase in functional connectivity in the beta band for the motor task was not present, as this is usually restricted to sensorimotor areas. Figure 1 C shows brain regions with significant effects in non-reversibility in the beta band. Though the effect for sensorimotor regions was most pronounced, a clear increase in non-reversibility is also observed in the visual areas and premotor areas. The same was true for both the theta and alpha band (see Figure S1). For the gamma band, significant increase in non-reversibility was predominantly found in bitemporal regions. Lastly, using a random forest classifier and data from all frequency bands, we demonstrate that classification of task condition by functional connectivity was outperformed by non-reversibility based classification ( < 0 . 001 ; Fig. 1 D).
2. Language task. For the language task (shown only for the auditory narratives), we also observed strong effects for whole-brain nonreversibility in different frequency bands (delta, alpha and gamma; Fig. 1 B). Similarly as for the motor task, non-reversibility is more sensitive to detect task induced effects as difference in whole-brain functional connectivity between task and resting-state condition was only found for the alpha and gamma band and not for the delta band. Again the direction of the effect is similar for functional connectivity and nonreversibility. The language task activates temporal and frontal language related brain areas and regional analysis indeed shows significant effects in non-reversibility in these areas ( Fig. 1 C). At the same time, there is deactivation in regions corresponding to the posterior default mode in the alpha band ( Figure S1). Especially in the gamma band, widespread increases in non-reversibility are observed in bilateral frontotemporal regions. Finally, classification for task condition did not show superior classification accuracy for non-reversibility compared to functional connectivity ( > 0 . 05 ; Fig. 1 D). A potential explanation is that language induced effects are restricted to fewer frequency bands for both nonreversibility and functional connectivity compared to other tasks, which could result in similar classification accuracy for non-reversibility and functional connectivity to identify the task.
3. Working memory task. For the working memory task we also observed significant effects for non-reversibility across frequency bands (delta, alpha, beta, gamma; Fig. 1 B), while whole-brain alterations in functional connectivity were restricted to the delta and alpha band, with an increase in functional connectivity in the delta band and a decrease in functional connectivity in the alpha band. Again, the direction of change in non-reversibility is consistent with the direction of change of functional connectivity, and all task induced effects captured by functional connectivity are also identified using non-reversibility. Working memory usually elicits activation of frontal regions in the theta band. Although a significant effect between task and resting-state was not apparent in whole-brain non-reversibility or whole-brain functional connectivity for this frequency band, regional analysis revealed strong presence of increased non-reversibility in frontal theta regions ( Fig. 1 C). For the alpha band, there was a strong decrease in alpha band non-reversibility in occipital areas as expected. For the beta band, we observed higher non-reversibility in the left sensorimotor regions due to right hand button press involved in the task ( Figure S1). Similarly as for the motor task, task condition was better classified using non-reversibility than functional connectivity ( < 0 . 001 ; Fig. 1 D).

Non-reversibility in MEG reconstructed null-data
Our second step was to test the null-hypothesis of reversibility in MEG data using null-data. We observe for all tasks and frequency bands that the null-hypothesis of reversibility could be rejected ( Fig. 2 ). This indicates that the temporal asymmetry in the cross-correlation of the amplitude envelopes is a genuine feature of the data and does not reflect statistical noise. For resting-state data, the null-hypothesis of reversibility could also be rejected for all frequency bands, however, the differences in the distributions between the non-reversibility in the observed data and the surrogate or null data were less pronounced for resting-state data compared to task data.

Whole-brain modelling of non-reversibility in MEG
We lastly investigated contributing factors to non-reversibility in MEG data using neural mass modelling, with two potential candidates, asymmetry in the effective connectivity and heterogeneous axonal conduction delays. Individual Wilson-Cowan oscillators were coupled using the effective connectivity rather than the structural connectivity. We first did not include axonal conduction delays and ran our optimisation of effective connectivity with ′ = 0 as benchmark. Effective connectivity is in this case merely optimised by functional connectivity, and hence no asymmetry is introduced. Figure 3 A shows that although simulated Fig. 1. Classification of functional connectivity vs non-reversibility . Whole-brain non-reversibility is depicted as a function of lag or delay for different frequency bands for the motor task (Panel A). Every line depicts the behaviour of non-reversibility for one subject. Note clear peaks with frequency specific maxima for the delays. Panel B shows the whole-brain non-reversibility for the lag corresponding to the maximum non-reversibility for every subject (dot in the distribution). * refers to < 0 . 01 and * * refers to < 0 . 001 . Panel C shows non-reversibility for brain regions that showed significant difference or contrast between task and resting-state (FDR corrected). Non-significance is depicted by grey regions. Random forest based classification of task condition is depicted in panel D with whole-brain functional connectivity or whole-brain non-reversibility as features or input for the classification.

Fig. 2.
Non-reversibility in MEG reconstructed null-data . Whole-brain non-reversibility for the lag corresponding to the first maximum of non-reversibility for every subject (dot in the blue distribution) is shown. The same is depicted for null-data (in red). A star * refers to significance of < 0 . 001 . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 3. Whole-brain modelling of non-reversibility . We simulated whole brain non-reversibility using Wilson-Cowan oscillators and by making use of generative effective connectivity. Panel A shows the Pearson correlation between the simulated FC (AEC) and empirical FC (red), and the Pearson correlation between the simulated non-reversibility (NR) and empirical non-reversibility (blue) as the function of iterations during the optimisation. Panel B shows the empirical NR matrix and the simulated NR matrices. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) FC can adequately match empirical FC when ignoring asymmetry in the effective connectivity in the simulations (Pearson correlation between simulated and empirical FC reaches 0.85), non-reversibility is poorly reconstructed. This can also be visually inspected by the lack of structure in the simulated non-reversibility matrix for this case and its lack of resemblance with the empirical non-reversibility matrix ( Fig. 3 B).
After introducing asymmetry in the effective connectivity (by setting ′ ≠ 0 ), we see a similar level of fit for empirical FC. The Pearson correlation between simulated and empirical FC again approximates 0.85. In contrast to the previous case, the level of fit for empirical non-reversibility deviates from zero and reaches approximately 0.2 ( Fig. 3 A). This can also be observed from visual inspection of the empirical and simulated NR matrices ( Fig. 3 B). A similar goodness-offit can be observed when heterogeneous axonal conduction delays are included. When introducing both asymmetry and heterogeneous axonal conduction delays, again no improvement on the level of fit for FC could be observed. However, there was a clear effect on the level of fit for empirical non-reversibility that reaches 0.25 (see also the visual similarity between the empirical and corresponding simulated NR matrix). Hence, these results show that both asymmetry in the effective connectivity and heterogenous axonal conduction delays contribute to non-zero nonreversibility.

Discussion
We adopted a recently introduced model-free framework to study breaking of temporal symmetry of MEG amplitude envelope data to identify task condition in comparison to conventional functional connectivity. This framework characterises breaking of detailed balance, a hallmark of any living system and is rooted in thermodynamics. Indeed, non-reversibility derived from the lagged amplitude envelope correlations (LAEC) outperformed conventional connectivity in characterisation of task condition. Non-reversibility revealed rich spatiotemporal structure across different task conditions and resting-state. Well-known task-induced spatial and frequency specific signatures were retrieved for non-reversibility such as activation of sensorimotor cortices during a motor task in the beta band, orbitofrontal cortices during a working memory task in the theta band and language related temporal and frontal areas during a language task in the delta band. Moreover, using null-data we identified that non-reversibility was a genuine characteristic of MEG data and could not be obtained from a symmetric or reversible system. Lastly, using neural mass simulations we demonstrate that asymmetry in effective connectivity and heterogenous axonal conduction delays play a major role in shaping non-reversibility in MEG data.
Key result in our work is the notion that we could reject the nullhypothesis of reversibility in MEG data. This was not only the case for different task conditions, but also for resting-state MEG data. Estimation of non-reversibility in our data and in previous work ( Deco et al., 2022a;Kringelbach et al., 2023 ) comes with small magnitudes. Our surrogate data could provide information about the lower limit for these small non-reversibility values. In addition, non-reversibility showed task induced modulations, which were both frequency specific and in agreement with expected spatial activation maps. Cognitive neuroscience has a long tradition to study bottom-up and top-down processes in the brain ( Sarter et al., 2001 ). This very much relates to the concepts of extrinsic and intrinsic brain dynamics that can be captured using non-reversibility (i.e. the former is the effect of the extrinsic environment on the brain and the latter internally driven brain dynamics) ( Deco et al., 2022a ). The advantage of MEG measurements is that top-down and bottom-up processes usually involve spatiotemporal activation in distinct frequency bands, such as top-down processes in the alpha band and bottom-up processes in the gamma band ( Jensen et al., 2014 ). Our work shows that non-reversibility is more sensitive to bottom-up processes in the gamma band than conventional functional connectivity. Overall, it seemed that the more a task enforces a participant to be engaged with the environ-ment, the higher the non-reversibility, especially for the gamma band. Furthermore, our results showed that classification of task condition was more accurate using non-reversibility than conventional functional connectivity, and hence, non-reversibility could pave the way for a more detailed characterisation of top-down and bottom processes, which remains challenging using conventional functional connectivity.
Similar as in previous fMRI work , asymmetry in effective connectivity is an important causal entity for the emergence of non-reversibility in MEG data. Even though correlations between simulated and empirical non-reversibility were only moderate after introducing asymmetry in effective connectivity, these correlations were absent when asymmetry in effective connectivity was ignored. While it can be hypothesized that temporal asymmetry in amplitude envelopes could partially be induced by heterogeneous axonal conduction delays, addition of this entity to the simulations also resulted in a good the level of fit between empirical and simulated nonreversibility. Also note that previous work using the same functional connectivity metric and conventional structural connectivity informed modelling did not reach this level of fit between simulated and empirical functional connectivity ( Cabral et al., 2014;Tewarie et al., 2019a ). One potential way to further improve the fit for empirical non-reversibility patterns is to include regional heterogeneity in the simulations. Recent work has demonstrated that regional heterogeneity in large scale brain models can greatly improve the fit of empirical functional connectivity and future work could analyse whether the same would hold for nonreversibility ( Deco et al., 2021a ).
A few methodological issues should be acknowledged. First, nonreversibility has only been assessed on amplitude envelope data rather than phase data. Applying non-reversibility to phase locking or coherence methods is straightforward as it does not require the transformation of coupling values to mutual information. However, we leave the implementation implementation of non-reversibility to phase locking and coherence methods for future work. Second, trials for task data were not divided into pre-and post-stimulus periods as task induced effects on non-reversibility were clearly visible even without this separation. Third, we left out time-frequency representation of the data as we consider these to be well-known for the general reader. Fourth, we have compared functional connectivity and non-reversibility in terms of classification accuracy for specifying task condition. However, from a mechanistic or electrophysiological viewpoint it could well be the case that these different concepts contain complementary information. Fifth, we have compared functional connectivity to non-reversibility as a first step and proof of concept. However, future work could compare classification accuracy for task condition based on non-reversibility versus other lag-based measures such as Granger causality or transfer entropy ( Friston et al., 2013;Vicente et al., 2011 ). We would like to stress that the standard implementation of these lag-based measures are unlike non-reversibility insensitive to temporal asymmetry in functional interactions. Lastly, recent work shows that the power envelope of neuronal oscillations is characterised by positive kurtosis and positive cokurtosis (Hindriks et al., 2023). We have therefore used the Wilson-Cowan model rather than the normal form of the supercritical Hopf bifurcation as it is more straightforward to capture this empirical phenomenon of positive cokurtosis when additive coupling is used.
We showed that non-reversibility is insensitive to residual mixing of the source-reconstructed MEG signals. Hence, there is no need to apply pairwise orthogonalisation prior to calculation of non-reversibility. This makes it an ideal measure for assessing interactions, not only for MEG signals, but for electroencephalographic (EEG) and electrocorticographic (ECoG) signals as well, which suffer from the same mixing problem ( Schoffelen and Gross, 2009 ). Non-reversibility can hence be added to the list of mixing-insensitive interaction measures that are relatively insensitive to primary leakage (though not secondary leakage), such as the imaginary coherence ( Nolte et al., 2004 ) and the (weighted) phaselag index ( Stam et al., 2007;Vinck et al., 2011 ). Unlike these measures, which can only be applied to complex signals, non-reversibility can be applied to real signals as well and this allows to study functional connectivity in broadband signals. Thus, non-reversibility is likely to be a useful measure for analysing interactions in different experimental scenarios.
In conclusion, we have adopted the new non-reversibility framework derived from the lagged amplitude envelope correlation (LAEC) to analyse task-induced brain states in MEG data. Non-reversibility is a genuine characteristic of MEG data and outperforms conventional functional connectivity in classification of task conditions. Whole-brain computational modelling demonstrates that non-reversibility emerges when two neuronal populations are exposed to asymmetry in connection strengths. Furthermore, this new framework opens avenues to investigate bottomup and top-down process in cognitive neuroscience.

Data and code availability statement
Data used in the manuscript is available from the Human Connectome Project database. Code can be found at https://github.com/ Prejaas/MEG _ nonreversibility

Declaration of Competing Interest
None of the authors report any conflict of interest.

Data availability
We have used open-access MEG data obtained from the Human Connectome Project.