{"paper_id":"077ae16b-0cd8-4f7c-9b6c-a84f68e78333","body_text":"PREPRINT\n1\n2\n3\n4\n5\n6\n7\n8\n9\n10\n11\n12\n13\n14\n15\n16\n17\n18\n19\n20\n21\n22\n23\n24\n25\n26\n27\n28\n29\n30\n31\n32\n33\n34\n35\n36\n37\n38\n39\n40\n41\n42\n43\n44\n45\n46\n47\n48\n49\n50\n51\n52\n53\n54\n55\n56\n57\n58\n59\n60\n61\n62\n63\n64\n65\n66\n67\n68\n69\n70\n71\n72\n73\n74\n75\n76\n77\n78\n79\n80\n81\n82\n83\n84\n85\n86\n87\n88\n89\n90\n91\n92\n93\n94\n95\n96\n97\n98\n99\n100\n101\n102\n103\n104\n105\n106\n107\n108\n109\n110\n111\n112\n113\n114\n115\n116\n117\n118\n119\n120\n121\n122\n123\n124\nNonlinear brain connectivity from neurons to\nnetworks: quantification, sources and localization\nGiulio Tani Raffaellia, Stanislav Jiˇr´ ıˇceka,b,c, and Jaroslav Hlinkaa,b,1\nThis manuscript was compiled on November 17, 2024\nSince the first studies in functional connectivity, Pearson’s correlation has been the primary\ntool to determine relatedness between the activity of different brain locations. Over the years,\nconcern over the information neglected by correlation pushed toward using different measures\naccounting for non-linearity. However, some studies suggest that, at the typical observation\nscale, a linear description of the brain captures a vast majority of the information. Therefore,\nwe measured the fraction of information that would be lost using a linear description and\nwhich regions would be affected the most. We considered fMRI, EEG, iEEG, and single unit\nspikes to assess how the observation scale impacts the amount of non-linearity. We observe\nthat the information loss is reduced for modalities with large temporal or spatial averaging\n(fMRI and EEG) and gains relevance on more fine descriptions of the activity (iEEG and single\nunit spikes). We conclude that for most human applications, Pearson’s correlation coefficient\nadequately describes pairwise interactions in time series from current recording techniques.\nFunctional Connectivity | nonlinearity | Mutual Information | fMRI | Electrophysiology\nM\nany complex real-world systems, such as the human brain, Earth’s climate,\nor social or financial networks, do not easily allow probing their full\nintrinsic repertoire. Fully controlled laboratory experiments are unattainable,\nand the researchers increasingly complement traditional theoretical or experimental\napproaches with advanced data analysis on observational recordings of their natural\nbehaviour. A common starting step in such endeavour is to capture the structure\nof statistical interactions between their constituent units. These interactions are\nthen further analysed, e.g., through graph theory, machine learning, or statistical\napproaches.\nIn the field of neuroscience, the term functional connectivity is used for such\nstatistical dependence between remote neurophysiological events ( 1, 2). It has been\nthe key cornerstone of the paradigm shift stressing the importance of understanding\nthe mechanisms of spontaneous brain activity dynamics (\n3, 4). Over the decades, a\nplethora of functional connectivity indices has been proposed and applied. These\nindices form a densely populated zoo of methods, differentially sensitive to various\nforms of the concerned statistical dependence. This variety ranges from the classical\nPearson’s correlation coefficient, assessing the strength of linear dependence between\nobserved subsystems, to the principled use of mutual information. This is an entropy-\nbased measure sensitive to any form of dependence, making it highly attractive for\nstudying the interaction structures of highly non-linear or even chaotic systems.\nOver the years, a number of papers have reported the effectiveness of MI on different\nmodalities such as EEG (5) or fMRI (6) and significant changes in MI that correlate\nwith disease (7, 8).\nMutual information as a measure of dependence may thus appear as a clear\nmethod of choice due to its lack of theoretical assumptions. However, in practice,\nthere are trade-offs that need to be taken into account. The issue of estimate\naccuracy, the need for longer time series (9), as well as computational demands, may\nfavour more basic methods such as the Pearson’s linear correlation coefficient. The\npragmatic question the researchers (should) face is thus: Is the dependence pattern\nunder study so exotic and non-linear as to warrant using more general methods (or\neven, in principle, mutual information) at the expense of decreased statistical power,\nincreased computational demands, and other costs of more advanced dependence\nquantification methods?\nWhile the presence of non-linearity in brain signal dependence structure is\ntheoretically well established and beyond dispute, much less is known about its\nstrength. A principled approach for its quantification in terms of extra-Gaussian\ninformation has so far only been applied to a single modality (10). The investigation\nwas limited to functional magnetic resonance imaging (fMRI), and at a single\nSignificance Statement\nIn neuroimaging, as in other com-\nplex systems fields, the increasing\ninterest in network inference by sta-\ntistical dependencies (i.e. functional\nconnectivity) invites advanced ways\nto quantify it. Various nonlinear\nmeasures, ultimately the Mutual\nInformation, are used as alterna-\ntives to conventional linear Pear-\nson’s correlation coefficient. We\nsystematically assess the amount\nand reliability of detectable non-\nlinearity of brain functional connec-\ntivity across imaging modalities and\nspatial scales. We demonstrate\nmore pronounced nonlinearity in\nmicro-scale recordings, while rather\nlimited and unreliable in more ac-\ncessible, noninvasive, large-scale\nmodalities: functional magnetic res-\nonance imaging and scalp electro-\nphysiology. This fundamentally sup-\nports the use of robust and easily\ninterpretable linear tools in large-\nscale neuroimaging, and important\ninsights concerning the microscale\nconnectivity nonlinearity, including\nthe link to brain state dynamics.\nAuthor affiliations: aInstitute of Computer Science,\nCzech Academy of Science, Prague, Czech Republic;\nbNational Institute of Mental Health, Klecany, Czech\nRepublic; cFaculty of Electrical Engineering, Czech\nTechnical University in Prague, Czech Republic\nGTR and JH designed the methodology and handled\nthe validation and formal analysis. JH conceptualised\nand supervised the study, secured the funding, cared\nfor the administration and managed resources. GTR\nwrote the software for the analysis and produced the\nvisualisation and the first version of the draft. All\nauthors conducted the investigation, curated the data,\nand reviewed and edited the draft.\nNo competing interests.\n1To whom correspondence should be addressed. E-\nmail: hlinka@cs.cas.cz\nPrepr\nint — November 17, 2024 — 1–9\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted November 18, 2024. ; https://doi.org/10.1101/2024.11.17.623635doi: bioRxiv preprint \n\nPREPRINT\n125\n126\n127\n128\n129\n130\n131\n132\n133\n134\n135\n136\n137\n138\n139\n140\n141\n142\n143\n144\n145\n146\n147\n148\n149\n150\n151\n152\n153\n154\n155\n156\n157\n158\n159\n160\n161\n162\n163\n164\n165\n166\n167\n168\n169\n170\n171\n172\n173\n174\n175\n176\n177\n178\n179\n180\n181\n182\n183\n184\n185\n186\n187\n188\n189\n190\n191\n192\n193\n194\n195\n196\n197\n198\n199\n200\n201\n202\n203\n204\n205\n206\n207\n208\n209\n210\n211\n212\n213\n214\n215\n216\n217\n218\n219\n220\n221\n222\n223\n224\n225\n226\n227\n228\n229\n230\n231\n232\n233\n234\n235\n236\n237\n238\n239\n240\n241\n242\n243\n244\n245\n246\n247\n248\nspatial signal averaging resolution given by the common,\nyet relatively coarse and arguably suboptimal Automatic\nAnatomical Labelling atlas ( 11). The results suggested that\nfor this modality and at this spatial resolution, non-linear\ndependence contribution was likely negligible, supporting\nthe practice of using Pearson’s linear correlation for analysis\nof this type of coarse fMRI data. However, this still left\nthe question of the relevance of non-linear dependences\nunanswered for a dominant proportion of neuroscientists,\nas the results could be only speculatively generalized to other\nneuroimaging modalities and spatial scales.\nTo at least partially fill this gap, we here apply the\nprincipled methodology for non-linearity quantification to\nneuroimaging data ranging across modalities from fMRI\nthrough scalp electroencephalography, intracranial electroen-\ncephalography, to single neuron (spike) recordings, adding\nanalysis across orders of magnitudes of spatial averaging, and\nadditional study of the sources, reliability, and localization\nof the observable non-linearity. Moreover, we provide an\nimplementation in open code, allowing to quantify non-\nlinearity in other datasets in neuroscience and beyond, which\nshould provide a powerful tool also to researchers in other\ndisciplines, following early successes of similar methodology\nin climatology (12) and financial networks (13).\nResults\nProbably the simplest, and in some sense canonical, form of\ndependence between quantitative variables is that of purely\nlinear relation between two variables in a jointly normal\n(Gaussian) distribution, the strength of which is captured\nby the well-known Pearson’s (linear) correlation coefficient.\nIndeed, linear correlation is sufficient to determine mutual\ninformation in a bivariate Gaussian distribution through the\nsimple equation\nIGauss = −1\n2 log\n(\n1 −r2)\n. [1]\nIt is exactly the deviation from the jointly Gaussian distri-\nbution that makes the use of Pearson’s linear correlation\ncoefficient potentially suboptimal, calling for the use of\nmore advanced methods such as the universally sensitive\nmutual information. The strength of the deviation from\nGaussian dependence pattern is thus the object of our interest.\nMoreover, we are not concerned with non-Gaussianity of the\nmarginal distributions individually (as these can be easily\ntreated by the use of monotonic nonlinear rescaling such as\nlog-transform, or application Spearman’s correlation instead\nof Pearson’s), but particularly with such deviations from\nGaussianity which concern the very relation between the\nvariables, called copula, which entails the full characterization\nof the dependence structure, invariant with respect to any\nsuch bijective monotonic rescaling of marginals.\nNote that to avoid too technical language, we shall\nuse the terms “linearity’ and “Gaussianity” (of distribu-\ntion/dependence/interaction/pattern/process) interchange-\nably throughout this manuscript, although indeed in a strict\nsense the linearity is a bit wider term in particular contexts\n(one can e.g. imagine linear functional dependence as the best\nfit between two variables with non-Gaussian distributions\nor errors, or linearly coupled process with nonlinear driving\nnoise).\nWe shall leverage that known statistical physics results\nthat the bivariate Gaussian has the maximum entropy, and\nthus the minimal information (under Gaussian marginals),\namong the distributions with the same correlation, and thus\nevery joint distribution that doesn’t have a Gaussian copula\n(we shall call these distributions “non-linear” ), has higher\nMI than expected from correlation by the formula Eq. (1).\nThis allows us to define the distribution non-Gaussianity by\nIextra(X,Y ) =I(X,Y ) −IGauss(r(X,Y )), or its normalized\nvariants.\nWe identify two primary sources of non-linearity: intrinsic\nnon-linearity and non-stationarities. Intrinsic non-linearity\nhappens when the recorded samples are genuinely identically\ndistributed, but they do not have a Gaussian joint probability\ndistribution. For example, this can include a relationship\nbetween absolute values or out-of-phase synchronization\nphenomena.\nWe refer to non-stationarity when the samples are not\nidentically distributed, and the source distribution depends\non time. The samples will, in general, not be distributed\naccording to a multivariate Gaussian, even if the source\ndistributions are all Gaussians. Most estimators require\nthe samples to be i.i.d. or at least from a stationary\ndistribution. The non-stationarity, even of resting-state (rs)\ndata, is well known and studied on its own as a potential\nsource of insight into the brain functioning ( 14, 15). The\neffect’s magnitude depends on the estimator and the non-\nstationarity’s properties. Examples of non-stationarities are\nthe switching between states—each associated with a different\ncorrelation, isolated bursts of activity, and a continuous drift\nof mean, variance, or correlation.\nAlong with these neural sources for the observed non-\nlinearity, others can reside in the acquisition and pre-\nprocessing of the signal. For instance, non-monotonous\ntransformations may easily result in an increase in observed\nnon-linearity, and in particular, observation of nonlinearity\nbetween originally linearly related processes (it is illustrative\nto imagine the effect of taking, as a domain-specific prepro-\ncessing step, absolute value or square of Gaussian signals\nbefore probing their functional connectivity; similar but less\nstraightforward effect is obtained by working on, e.g., Hilbert-\ntransformed data or (band-limited) signal power time series,\ntransformation much more common in neuroscience).\nStrength of non-linearity.We begin by measuring the fraction\nof information in all connections not explained by a bivariate\nGaussian copula. We call this measure Relative Non-Linearity\n(RNL).\nWe start from rs-fMRI and look at RNL over a range of\nregion numbers and sizes to probe the effect of different\ndegrees of spatial averaging. We observe a consistent\npresence of non-linearity for the different region sizes from\nthe Craddock atlas (\n16). The fraction of MI not explained\nby the correlation (Fig. 1A) sits around 4% for all region\nsizes except for very few and large regions. At the same time,\nthe non-linearity observed in the shadow (phase randomised)\ndataset remains below 2%, providing a measure of the bias\nin the absence of non-linearity. The difference in MI between\nempirical and shadow datasets is always significant, except for\nthe atlas with ten regions (when applying strict Bonferroni\ncorrection for multiple comparisons).\n2 — Tani Raffaelli et al.\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted November 18, 2024. ; https://doi.org/10.1101/2024.11.17.623635doi: bioRxiv preprint \n\nPREPRINT\n249\n250\n251\n252\n253\n254\n255\n256\n257\n258\n259\n260\n261\n262\n263\n264\n265\n266\n267\n268\n269\n270\n271\n272\n273\n274\n275\n276\n277\n278\n279\n280\n281\n282\n283\n284\n285\n286\n287\n288\n289\n290\n291\n292\n293\n294\n295\n296\n297\n298\n299\n300\n301\n302\n303\n304\n305\n306\n307\n308\n309\n310\n311\n312\n313\n314\n315\n316\n317\n318\n319\n320\n321\n322\n323\n324\n325\n326\n327\n328\n329\n330\n331\n332\n333\n334\n335\n336\n337\n338\n339\n340\n341\n342\n343\n344\n345\n346\n347\n348\n349\n350\n351\n352\n353\n354\n355\n356\n357\n358\n359\n360\n361\n362\n363\n364\n365\n366\n367\n368\n369\n370\n371\n372\n10 30 50 70 100150200230270300350400450500550600650700750800850900950\n# regions\n0.1\n0.0\n0.1\n0.2RNL\nA\n0.0\n0.1\n0.2RNL\nB\n0.0\n0.1\n0.2RNL\nC\nBand\n0.0\n0.1\n0.2RNL\nD\nEmpiric\nShadow\nNaïve\nShadow\nSignificative\ndifference\n0%\n2%\n5%\nFig. 1. Distribution over subjects of the Relative amount of Non-Linearity (RNL). A) fMRI data using Craddock parcellation and a varying number of regions, B) EEG scalp\nvoltage, C) iEEG voltage, D) EEG band-limited power. The shadow dataset is a linear surrogate of the empirical one, see text for detail.\nTani Raffaelli et al. Preprint — November 17, 2024 — 3\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted November 18, 2024. ; https://doi.org/10.1101/2024.11.17.623635doi: bioRxiv preprint \n\nPREPRINT\n373\n374\n375\n376\n377\n378\n379\n380\n381\n382\n383\n384\n385\n386\n387\n388\n389\n390\n391\n392\n393\n394\n395\n396\n397\n398\n399\n400\n401\n402\n403\n404\n405\n406\n407\n408\n409\n410\n411\n412\n413\n414\n415\n416\n417\n418\n419\n420\n421\n422\n423\n424\n425\n426\n427\n428\n429\n430\n431\n432\n433\n434\n435\n436\n437\n438\n439\n440\n441\n442\n443\n444\n445\n446\n447\n448\n449\n450\n451\n452\n453\n454\n455\n456\n457\n458\n459\n460\n461\n462\n463\n464\n465\n466\n467\n468\n469\n470\n471\n472\n473\n474\n475\n476\n477\n478\n479\n480\n481\n482\n483\n484\n485\n486\n487\n488\n489\n490\n491\n492\n493\n494\n495\n496\nTime sample [ms]\n# units per site\n125 250 5001000200040008000\n1\n2\n4\n8\n16\n32\nA 30 min. Empiric\n125 250 5001000200040008000\nB 30 min. Shadow\n# units per site\n125 250 500\n1\n2\n4\n8\n16\n32\nC 15 min.\n125 250 500\nD10 min.\nTime sample [ms]\n125 250 500\nE 6 min.\n125 250 500\nF 5 min.\n125 250 500\nG 3 min.\nRNL\n0.0\n0.1\n0.2\n0.3\n0.4\nFig. 2. Average RNL across mice varying the time average window size and the\nnumber of units. A-B) RNL computed over the full duration of the resting state\nexperimental block, C-G) average of the empirical RNL computed over epochs of\ndecreasing length.\nIn EEG time series, we observe a significantly higher pres-\nence of non-linearity (one-sided t-test, Bonferroni corrected)\nthan in the shadow dataset for all frequency bands. In the α\nband, the amount is similar to what was observed in fMRI\n(Fig. 1B), while in the others, it tends to be even lower,\nin particular under 2%. The different number of samples\navailable in each band explains the dependence between the\nfrequency band and RNL in the Shadow dataset. See the\nMethods section and the SI for details.\nIn iEEG, non-linearity is more prominent than in the\nprevious two datasets. For all bands, the empirical RNL\ntends to be above 5%. At the same time, in the shadow\ndataset, it is confined below 1% (Fig. 1C). Part of the non-\nlinearity is explained by non-stationarities due to interictal\nepileptic activity (see SI). Moreover, the iEEG was acquired\nduring natural activity (neither standard resting state nor\nsustained task), which might affect the stationarity of the\nsequences and their RNL.\nLastly, we looked at the RNL in single-unit spike rates\nin mice. In this case, we report the RNL averaged across\nsessions. The relative contribution of non-linearity is above\n23% (Fig. 2A) in all cases, over the full range of temporal\n(between 125 ms and 8 s) and spatial/group size (between\n1 and 32 units) averaging. We observe the highest values\nof RNL for the smallest group sizes and faster time scales.\nConversely, averaging over larger groups always reduces the\nRNL, and averaging over more than 2 s has little effect\nwith the RNL between 32% and 37% depending on group\nsize. We found no correlation between the RNL and the\nnumber of regions in a given mice (ranging between 9 and\n16). Thus, we present the results as the average across all\nmice. At the same time, the RNL never surpasses 5% in the\nshadow dataset (Fig. 2B). The magnitude and amplitude of\nthe fluctuations increases for shorter time series as observed\nfor lower frequency bands in EEG and iEEG.\nLocalisation. The amount of non-linearity is reduced for the\nmore accessible modalities. In fMRI and EEG, high temporal\nor spatial averaging masks most of the non-linearity. However,\nit might still be relevant if localised in specific regions. The\nauthors in ( 17) suggest localisation of non-linearity in the\noccipital region. We evaluated non-linearity’s localisation,\nlooking at regions that participate in consistently non-linear\nconnections across subjects.\nFor fMRI data, we used the AAL90 atlas as it provides\nexplicit anatomical labelling for the regions where the non-\nlinearity might be localised. This dataset has a significant\ncorrelation (p =.007) between the localisation of non-linearity\nin the empirical and shadow datasets (see SI). Furthermore,\ncorrelation and localisation increase with reduced denoising\nsteps in preprocessing. This suggests that most region-\nspecific non-linearity is due to artefacts and is removed during\npreprocessing.\nEEG data (Fig. 3) offer a different picture. Here, non-\nlinear relationships are quite limited in the band θ(only 38%\nchannel pairs) while present in more than 95% of channel\npairs in other bands. In the shadow dataset, less than\n0.8% of the relationships are consistently non-linear across\nsubjects. At the same time, each region participates in\nconnections with a non-random amount of non-linearity. For\nthe slower bandsδ,θ, andα, we observe that the degree of the\noccipital regions is higher (and in frontal regions lower) than\nexpected from a random graph, suggesting a localisation of\nnon-linearity. Conversely, for the high-frequency bandsβand\nγ, the results show a predominance of non-linearity in frontal\nand temporal electrodes. This anterior-posterior gradient\nmay be partially attributed to the spatial distribution of\nmost prevalent sources in normal awake EEG ( 18), however\nmore research is warranted.\nReliability. However small, accounting for non-linearity may\nstill be beneficial if the additional information is stable over\nrepeated measures. As the last test, we looked into non-\nlinearity’s reliability across sessions. We compared the subject\nranks for each connection based on Total Mutual Information\n(TMI, i.e., accounting for non-linearity) and those derived\nfrom correlation.\nAs shown in Figure 1A, the presence of non-linearity in\nfMRI data is significant for all atlases with enough regions.\nHowever, the amount is limited and unreliable under our\nmeasure. The reliability for the Pearson’s correlation is, on\naverage, more than twice that of TMI (0.277 against 0.134).\nMoreover, if we estimate MI from correlation, losing the sign\nof the relationship, we drop to similar levels of correlation\nacross sessions (0.153). In this case, however, Pearson’s\ncorrelation predicts TMI in a second session even slightly\nbetter (0.008 or 5% increase in correlation) than TMI predicts\nitself across sessions, showing the practical advantage of using\nlinear functional connectivity measure at this scale.\nRunning the same analysis on EEG and iEEG data offers\na similar picture (see SI), which again excludes a clear\nadvantage of TMI.\n4 — Tani Raffaelli et al.\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted November 18, 2024. ; https://doi.org/10.1101/2024.11.17.623635doi: bioRxiv preprint \n\nPREPRINT\n497\n498\n499\n500\n501\n502\n503\n504\n505\n506\n507\n508\n509\n510\n511\n512\n513\n514\n515\n516\n517\n518\n519\n520\n521\n522\n523\n524\n525\n526\n527\n528\n529\n530\n531\n532\n533\n534\n535\n536\n537\n538\n539\n540\n541\n542\n543\n544\n545\n546\n547\n548\n549\n550\n551\n552\n553\n554\n555\n556\n557\n558\n559\n560\n561\n562\n563\n564\n565\n566\n567\n568\n569\n570\n571\n572\n573\n574\n575\n576\n577\n578\n579\n580\n581\n582\n583\n584\n585\n586\n587\n588\n589\n590\n591\n592\n593\n594\n595\n596\n597\n598\n599\n600\n601\n602\n603\n604\n605\n606\n607\n608\n609\n610\n611\n612\n613\n614\n615\n616\n617\n618\n619\n620\nEmpiric\nShadow minimum\ndegree\nsignificantly\nbelow random\ngraph\nsignificantly\nabove random\ngraph\nmaximum\ndegree\nFig. 3. Degree of the regions in a network weighted by the z-score of TMI compared to surrogates. Light grey regions have degree zero. Dark grey regions have been excluded\nfrom the analysis as often missing or corrupted. The representation is a Voronoi tessellation of the stereographic projection on the xy-plane of standard electrode positions. The\nnasion is facing up.\nSources. To better understand the relevance of observed non-\nlinearity, we will now survey some sources of non-linearity\nthat may be considered spurious.\nEEG offers an example of non-linearity derived from the\nchoice of the observable. Let us consider band-limited power\nand compute the shadow dataset na¨ ıvely from the sequence\nof power values. This would be an acceptable choice as, for\nexample, in MEG, the computation of FC from the signal\nenvelope is customary (19, 20). We observe a large fraction of\nnon-linear information in all bands (Fig. 1D, yellow boxplots).\nHowever, when extracting the power from a surrogate of\nthe original EEG time series, we notice that most of the\nnon-linearity arises from power computation (Fig. 1D, blue\nboxplots).\nIn most bands, the non-linearity in the shadow dataset is\nstill lower than that of the empirical one. This suggests that\nthe transformation from voltage to power is not responsible\nfor the entirety of the observed nonlinearity in bandpower\ndependence. However, the large amounts of spurious non-\nlinearity mask the significance of band α(Fig. 1B).\nAs mentioned in the introduction and also discussed in\ndetail for climate systems elsewhere ( 12), non-stationarities\ncan be powerful sources of apparent non-linearity. An obvious\npotential source of non-linearity in the iEEG data is epileptic\nactivity. We observe how (Fig. 4), measuring the RNL on\na sliding window across the recording of a seizure, up to\nalmost 90% of the total information appears to be due to\nnon-linearity.\nIn particular, we observe how the RNL is sensitive to\nthe signal’s amplitude variation. The first peak for bands\nθto γhappens when the sliding window crosses between\nthe initial regime of low power and the first high amplitude\noscillations. The RNL decreases when the window contains\nonly large oscillations and increases again when the first group\nof electrodes reverts to lower power. The last increase in RNL\nhappens when the windows cross out of the seizure.\nA second example of the effect of non-stationarities comes\nfrom spiking data of individual neurons. Over the 30 minutes\nof recording, even without stimuli, the mice brain’s activity\n0 25 50 75100125150175200225250275\nTime [s]\npower\n(a.u.)Band Electrode \n 0.1\n0.2\n0.3\n0.4\n0.5\n0.6\n0.7\n0.8\nRNL\nFig. 4. Sample seizure showcasing the effect of non-stationarities on the measure\nof non-linearity. Top: traces from a subset of electrodes. Bottom: mowing window\nvalues of RNL and band-limited power. The green solid lines mark the beginning\nand the end of the seizure. The blue dashed lines mark the beginning of the first\nwindows containing samples past the green lines.\nTani Raffaelli et al. Preprint — November 17, 2024 — 5\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted November 18, 2024. ; https://doi.org/10.1101/2024.11.17.623635doi: bioRxiv preprint \n\nPREPRINT\n621\n622\n623\n624\n625\n626\n627\n628\n629\n630\n631\n632\n633\n634\n635\n636\n637\n638\n639\n640\n641\n642\n643\n644\n645\n646\n647\n648\n649\n650\n651\n652\n653\n654\n655\n656\n657\n658\n659\n660\n661\n662\n663\n664\n665\n666\n667\n668\n669\n670\n671\n672\n673\n674\n675\n676\n677\n678\n679\n680\n681\n682\n683\n684\n685\n686\n687\n688\n689\n690\n691\n692\n693\n694\n695\n696\n697\n698\n699\n700\n701\n702\n703\n704\n705\n706\n707\n708\n709\n710\n711\n712\n713\n714\n715\n716\n717\n718\n719\n720\n721\n722\n723\n724\n725\n726\n727\n728\n729\n730\n731\n732\n733\n734\n735\n736\n737\n738\n739\n740\n741\n742\n743\n744\nhad the opportunity to drift or switch through different states\n(most notably, immobile and running periods). If we compute\nthe RNL over shorter windows with fewer state changes\nand then average, we observe a clear reduction of the non-\nlinearity. However, even in the last case of 3-minute epochs,\nthe average RNL stays much above what is observed with\nother modalities.\nDiscussion\nThe fundamental non-linearity of neuronal activity leads to\na relevant fraction of non-linear information in single-unit\nspikes. However, moving to more accessible and non-invasive\nhuman data, these figures are vastly reduced. This might be\nan effect of the scale at which the activity is recorded. Indeed,\nfMRI and EEG imply large spatial or temporal averages.\nThe more significant fraction of non-linear information in\nthe iEEG dataset might support this interpretation. Part\nof this non-linearity may be explained by non-stationarities\nin the form of IEDs affecting the higher frequency bands\nand drift or state switching due to the natural activity of\nthe subjects during the recording. However, this modality\nremains more promising for the detection of significant non-\nlinearity in the human brain, and future work will have to\ntest this with improved control over non-stationarities.\nWhile the presence of some non-linearity is significant\nat all scales, its variability at the individual level makes it\nhard to assess its per-subject localisation, requiring group-\nlevel analysis. Indeed, the non-linearity is reduced and\nunreliable across sessions, which could be explained by a\nstrong contribution of noise in the difference between the\nGaussian and total MI.\nNon-linearity is present at all levels of our analysis, and the\noblivious use of Pearson’s correlation might be inappropriate.\nHowever, we advocated caution when designing studies that\nplan to apply MI. Without accurate control of spurious\nsources of non-linearity, the improvement from using MI\ninstead of correlation might be hardly reproducible, if any.\nLastly, many FC studies seek relatedness in the frequency\ndomain. This practice is widespread in EEG, where there is\nan established relationship between frequency and function.\nIn this context, there’s a similar debate around the need\nto go beyond Gaussian dependencies. Future work will be\nable to quantitatively (beyond pure statistical testing of the\nbinary question of presence or absence) assess the amount\nof non-linearity and the role of averaging in the frequency\ndomain.\nMaterials and Methods\nTo estimate the non-linear content in the relationships between\nregions, electrodes, and units, we compared the Total Mutual\nInformation (TMI) between time series to the estimate from\nsurrogates where only the linear relationships are preserved. This\nallows us to evaluate the global amount of non-linearity and which\nare the regions or electrodes where it is more substantial.\nMI estimator\nLooking for the potential benefits of using MI aligns with our\ndefinition of non-linearity as the deviation from a multivariate\nGaussian distribution in data. Indeed, of all distributions with\nGaussian marginals, a multivariate Gaussian is maximally entropic\nfor a given covariance matrix, i.e., has the minimum MI. Any\ndeviation from Gaussianity will yield a higher MI. Assuming\nGaussian marginal does not imply a loss of generality. Spearman\ncorrelation coefficient and MI are independent of the marginal\ndistribution. Moreover, while approximate Gaussian distribution\nis often assumed, every sample distribution can be mapped to\nGaussian marginals with a monotonous transformation. Indeed,\nwe enforced Gaussian marginals via rank-normalisation to ensure\nprecise non-Gaussianity estimates. We use the sample ranks to\nestimate the percentile πi for each sample xi and replace the\noriginal value with the one corresponding to the same percentile\nin the standard normal distribution N (0, 1).\nWe estimate the MI of two variables through equiquantal\nbinning (also known as equiprobable ( 21)): the samples are sorted\non a grid with the same bin number for both variables. The bins\nfor each variable have variable width so that the sum over the\nother variable always gives the same number of data points. The\nMI is then computed from the estimated probabilities of each bin\npij as the difference between the sum of the marginal entropies of\nthe two variables and their joint entropy:\nMI =−\n∑\ni\npi\nx logpi\nx−\n∑\nj\npj\ny logpj\ny +\n∑\nij\npij logpij [2]\nwith pi\nx =\n∑\njpij, pj\ny =\n∑\nipij, and pl\na≃pk\na\n, a = x,y ,∀l,k∈\n[1,N ] were the equality holds for all l and k only if the number of\nbins N is a divisor of the sequence length S. We chose N =⌊\n3√\nS⌋\nin line with the previously recommended pragmatic heuristic ( 22).\nThis estimate is known to be affected by bias. For high values of\nMI, the estimate is bounded above by the logarithm of the number\nof bins. By construction, in a perfect bi-univocal relationship, the\nnumber of bins in one dimension is also the number of non-empty\nbins in the joint distribution, all sharing the same number of points.\nThe higher the MI, the stronger the underestimate. On the other\nhand, due to the finiteness of the sample, the estimate has a positive\nbias (23). Any fluctuation will result in a greater than zero estimate\nof the MI, even for independent variables. The formulae for\nsmall sample sizes—that show good agreement with our estimates\nfor independent variables—are derived for independent samples.\nHowever, the definition of the estimator imposes bounds on the\nsums of rows and columns.\nThese biases have opposite signs and non-trivial tractability, and\nwe addressed them numerically. We evaluated the MI on samples\nfrom random bivariate Gaussian distributions with predetermined\ncorrelations (thus known mutual information values according to\nEq. (1)) and sample size S equal to the series length. Specifically,\nwe calculated MI for 50000 bivariate random samples of size S\nfor each correlation value ranging from 0 to .995 in increments\nof .005. The average of the 50000 MI estimates approximates\nthe expected sample mutual information for each given mutual\ninformation value. This process yields a monotonous function (with\nlinear approximation applied if needed to ensure monotonicity for\ncorrelations near zero) that relates the true mutual information to\nits expected numerical MI estimate. The inverse of this function\nthen allows us to adjust each estimated MI to produce a more\naccurate bias-corrected estimate of the true mutual information.\nWe derive the maps assuming that the number of bins N and\nsamples S and the true MI are the only determinant factors for\nthe bias. Every combination of a number of bins and a number of\nsamples requires a different map. This approximation holds for a\nwide range of observed relations spanning over 98% of our datasets’\nconnections.\nPreliminary analyses with KNN and KDE estimators didn’t\nshow any substantial change at the price of much slower estimation.\nOther kinds of estimators, such as those including embedding in\nhigher dimensional spaces ( 24, 25) were excluded because, even\nif they account for non-independent samples, they require even\nlarger samples for a correct distribution estimate.\nNon-linearity estimate. We assess non-linearity’s contribution to\nthe TMI, comparing it to the information conveyed only by linear\ncorrelations. For each dataset and subject, we calculated the MI\non 99 random realisations of multivariate time series that preserve\nthe linear structure but remove the non-linear components. If\nthe original time series has a Gaussian dependence structure, the\noriginal data MI would be similar to that of the surrogates, up\nto some random error due to the variance of the estimator and\n6 — Tani Raffaelli et al.\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted November 18, 2024. ; https://doi.org/10.1101/2024.11.17.623635doi: bioRxiv preprint \n\nPREPRINT\n745\n746\n747\n748\n749\n750\n751\n752\n753\n754\n755\n756\n757\n758\n759\n760\n761\n762\n763\n764\n765\n766\n767\n768\n769\n770\n771\n772\n773\n774\n775\n776\n777\n778\n779\n780\n781\n782\n783\n784\n785\n786\n787\n788\n789\n790\n791\n792\n793\n794\n795\n796\n797\n798\n799\n800\n801\n802\n803\n804\n805\n806\n807\n808\n809\n810\n811\n812\n813\n814\n815\n816\n817\n818\n819\n820\n821\n822\n823\n824\n825\n826\n827\n828\n829\n830\n831\n832\n833\n834\n835\n836\n837\n838\n839\n840\n841\n842\n843\n844\n845\n846\n847\n848\n849\n850\n851\n852\n853\n854\n855\n856\n857\n858\n859\n860\n861\n862\n863\n864\n865\n866\n867\n868\nsurrogates. Conversely, if the original data had substantially higher\nMI than the surrogates, this would indicate the presence of non-\nlinear dependencies.\nWe created the surrogates using the multivariate Fourier\ntransform (FT) method ( 26, 27), which generates realisations\nof multivariate linear stochastic processes preserving the individual\nspectra and cross-spectra of the original time series. Specifically,\neach surrogate is obtained by adding identical random phases\nto corresponding frequency bins of the series’s FT while keeping\nthe amplitude’s magnitudes unchanged. The series were then\ntransformed back into the time domain using the inverse FT. These\nsurrogates retain the dependency structure that a multivariate\nlinear stochastic process can explain while destroying “non-\nlinearity” .\nComparing the MI estimate of the data and of the “linear”\nsurrogates, rather than directly using the linear correlation of the\ndata, has two advantages. First, correlation and MI estimators\nhave different characteristics in terms of bias and variance, and\nsurrogates thus allow for a direct quantitative comparison between\nnon-linear and linear connectivity. Second, surrogates represent\na suitable null model for direct statistical testing of differences.\nHowever, while these estimates are valuable for hypothesis testing—\nas we did to assess localisation—we use the mean of these 99\nvalues when interested in the relative difference. We refer to\nthis as “Gaussian” MI, which closely approximates the MI of a\nbivariate Gaussian distribution. The Relative Non-Linearity (RNL)\nis defined as the fraction of the TMI that exceeds the Gaussian\nMI. Note that to obtain robust estimates of RNL, throughout the\npaper it is reported at a global level, i.e. before taking the fraction,\nboth TMI and Gaussian MI are first averaged over all pairs of\nsignals.\nAs the correction map depends on the number of points used\nto estimate the MI, using oversampled data would introduce some\nbias. For low MI—i.e., for most of the connections—the correction\nmap reduces the measured value. This necessary correction is\nsmaller when the map is computed for a larger sample (i.e. lower\nbias). Using an oversampled sequence would lead to using a high-\nsample-count map, while the non-independent samples would still\nbehave as if they were in a smaller number. As this applies to\nthe original data and the surrogates, both measures would be\noverestimated by different amounts depending on their MI. This\nwould lead to a negative bias in RNL.\nWe acknowledge that for fMRI data, given the band-pass filter,\nthe current sampling rate of 0.5 Hz may lead to oversampling.\nHowever, reducing the sampling rate to 0.18 Hz would give time\nseries too short to get any reliable estimate of MI. In this case,\nthe relative amount of non-linear information will have a small\nnegative bias. However, as this affects the shadow dataset (see\nbelow) similarly, every result relative to it remains valid, as will\nthe significance tests.\nShadow datasets. Following (10), we compared each result against\na control analysis using linear “shadow” datasets. This approach\nallows us to account for any potential biases in the generation of\nsurrogate distributions, such as those caused by small sample\nsizes. For each session, a shadow dataset was created as a\nmultivariate FT surrogate of the original, marginally normalised\ndataset, thereby preserving only the original data’s linear structure\n(both autocorrelation and cross-covariance).\nWe then applied the same processing to the original data\nand the shadow, including initial normalisation, generation of\nmultivariate surrogates, and computation of MI and RNL. This\napproach allowed us to replicate the entire procedure using a\ndataset with the same correlation and autocorrelation structure\nand purely linear interactions, ensuring that any potential bias\nin our findings due to the algorithm’s numerical properties was\nadequately controlled. We compared the RNL between the original\nand shadow datasets using a paired t-test. All group-level tests\napplied a family-wise corrected significance threshold of p = 0.05.\nFor Band-Limited Power (BLP) data, we considered two options\nfor generating the shadow dataset. In the na¨ ıve version, we extract\nthe power from each band and then surrogate the power values.\nThis procedure destroys the sequence’s non-linearity and those\nthat might have been introduced while computing the power. In\nthe second version, we generate a surrogate voltage sequence and\nthen extract new power values from this.\nAmount of non-linearity\nWe relied on session-wise averages in each dataset to get the global\nfraction of non-linearity, RNL. In pure linear relationships, the\ntotal information measured on data should be indistinguishable\nfrom surrogates. We estimate the relative amount of non-linear\ninformation in the brain—according to the different modalities—as\nthe relative difference between the average of the total MI across\nall sequence pairs and the average of the Gaussian MI across all\npairs and surrogates.\nDespite individual and experimental fluctuation or any system-\natic bias in the estimate, session-wise relative differences for linear\ndata should stay close to the estimate from the shadow dataset. A\nsignificant difference between empirical data and shadow datasets\nis the sign of the presence of non-linearity.\nLocalisation of non-linearity\nWe evaluate the localisation of non-linearity leveraging the statistics\nof the surrogates. If a region (or, equivalently, electrode) pair has\nlinear interaction, repeated measures across sessions or subjects\nshould give results analogous to the surrogates. This can be\nobserved by looking at the z-score of the measure of empirical data\ncompared to the mean and variance of surrogate measures. For\neach pair of regions, we average the TMI and the 99 surrogate MIs\nacross all subjects. We then use the mean and standard deviation of\nthe surrogates to compute a\nz-score for the empirical measure. We\nused Holm-Bonferroni correctedp-values over empirical and shadow\ndatasets together to assess which connections are significantly\nnon-linear across subjects. We use z-scores instead of the relative\nnon-linearity computed on individual pairs, as the latter introduces\na strong correlation between empirical and shadow datasets due\nto the shared variations in correlation.\nLastly, we treated the significant z-scores as weights in an\nundirected graph and checked how the node degree compares with\na random graph with the same link strength distribution. This\nallows us to visualise regions with significantly higher or lower\namounts of non-linear connections compared to chance.\nThis analysis is possible only on the fMRI and EEG data\nas other modalities do not comprehensively sample the whole\nbrain. We repeated the same steps on the fMRI data with raw\npreprocessing to highlight the role of artefact removal in the\nobserved non-linearity.\nReliablity of non-linearity\nFor the non-linearity to be potentially relevant as an individual\ncharacteristic (e.g. as a diagnostic biomarker), the additional\ninformation provided to the researcher must be reliable across\nsessions. We compared how well measures of FC through\ncorrelation or TMI in one session predict FC (estimated in either\nway) in a different session.\nWe compute the strength of the connectivity in three ways:\nPearson’s correlation of the rank-normalised values, TMI and the\nGaussian MI estimate through Pearson’s correlation. The latter is\nobtained according to Eq. (1), with r being the sample Pearson\ncorrelation over the rank-normalised sequence. This estimate of\ninformation is faster to compute than TMI but only accounts for\nlinear relationships as the Spearman correlation.\nWe assume that high or low connectivity between specific\nregions is a marker of scientific or clinical interest. In this case,\nwe desire that subjects ranking high or low in one session do so\nalso in another. For each connection, we compute the Spearman\ncorrelation of the strength of connectivity of the subjects over two\nsessions. We report averages across connections.\nSources of non linearity\nFor EEG, iEEG and single-unit spikes, we show examples of sources\nof non-linearity at the whole brain level. We compute the RNL\non BLP sequences for EEG, contrasting two ways of deriving the\nTani Raffaelli et al. Preprint — November 17, 2024 — 7\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted November 18, 2024. ; https://doi.org/10.1101/2024.11.17.623635doi: bioRxiv preprint \n\nPREPRINT\n869\n870\n871\n872\n873\n874\n875\n876\n877\n878\n879\n880\n881\n882\n883\n884\n885\n886\n887\n888\n889\n890\n891\n892\n893\n894\n895\n896\n897\n898\n899\n900\n901\n902\n903\n904\n905\n906\n907\n908\n909\n910\n911\n912\n913\n914\n915\n916\n917\n918\n919\n920\n921\n922\n923\n924\n925\n926\n927\n928\n929\n930\n931\n932\n933\n934\n935\n936\n937\n938\n939\n940\n941\n942\n943\n944\n945\n946\n947\n948\n949\n950\n951\n952\n953\n954\n955\n956\n957\n958\n959\n960\n961\n962\n963\n964\n965\n966\n967\n968\n969\n970\n971\n972\n973\n974\n975\n976\n977\n978\n979\n980\n981\n982\n983\n984\n985\n986\n987\n988\n989\n990\n991\n992\nshadow dataset. In iEEG data, we show the evolution of RNL\nthroughout a seizure. We used a sliding window of 45 s with 90%\noverlap for each frequency band. The average power is computed\nas the mean over all electrodes of the absolute value of the Hilbert\ntransform averaged over each window.\nFor single-unit spikes, we considered the sequences with higher\nsampling rates. This included firing rates averaged over windows\nof 125, 250 and 500 ms, allowing for enough samples to compute\nMI on shorter windows. We split the whole sequence into epochs\nof 15, 10, 6, 5 and 3 minutes, computed the RNL for each epoch\nand then averaged over all sessions and epochs.\nData\nThis work relies on five openly accessible datasets spanning four\nmodalities: fMRI, EEG, iEEG and single-unit spikes.\nfMRI data. We used two different datasets of fMRI data: the public\ndataset used in ( 28)(ESO245) and the MPI-Leipzig Mind-Brain-\nBody dataset (29, 30) (LEMON).\nESO245. This dataset contains 10 minutes of resting-state eyes-\nclosed functional magnetic resonance from 245 healthy subjects\n(148 right-handed, 132 females, mean age 29.22/standard deviation\n6.99) acquired as healthy controls as part of the ESO project.\nParticipants were informed about the experimental procedures and\nprovided written informed consent. The local Ethics Committee\nof the Institute of Clinical and Experimental Medicine and the\nPsychiatric Center Prague approved the study design. The\nacquisition included T1-weighted and T2-weighted anatomical\nscans not used in this study. The scanner was a 3T MRI\nscanner (Siemens; Magnetom Trio) at the Institute of Clinical\nand Experimental Medicine in Prague, Czech Republic.\nFunctional images were obtained using T2-weighted echo-planar\nimaging (EPI) with BOLD contrast. GE-EPIs (TR/TE = 2,000/30\nms) comprised 35 axial slices—acquired continuously in descending\norder covering the entire cerebrum (48 × 64 voxels, voxel size = 3\n× 3 × 3 mm3) (28).\nPreprocessing For the present study, we used the data already\nprocessed∗. We report the essential aspects of the processing while\nreferring to (28) for the details.\nThe preprocessing followed the default pipeline of the CONN\ntoolbox (McGovern Institute for Brain Research, MIT, USA) with\n12 head motion parameters and five white matter components in\nthe Component-based Correction (CompCor). The authors in ( 28)\nalso detrended the resulting time series and applied band-pass\nfiltering with cut-off frequencies 0.009-0.08 Hz. This pipeline is the\n“stringent” preprocessing compared to the “moderate” preprocessing\n(6 head motion parameters, one white matter component, cut-off\nfrequencies 0.004-0.1 Hz) and “raw” (no CompCor nor filtering).\nThe results in this paper use the “stringent” preprocessing unless\nexplicitly stated otherwise.\nChoice of the atlas We used two sets of atlases available in ( 28).\nThe first includes 23 different parcellations using the Craddock\natlas and a number of ROIs between 10 and 950 to investigate the\neffect of spatial averaging and region size on non-linearity. The\nsecond is the widely used AAL atlas with 90 regions to assess\nnon-linearity localisation.\nThe normalised cut spectral clustering used in Craddock\nparcellation can yield a number of ROIs smaller than desired—i.e.,\nempty clusters—and some ROIs can fall outside the GM mask for\nsome subjects. This effect is more pronounced with decreasing\nregion size and is observed in this dataset for all sizes larger than\n200. The largest number of ROIs is 840 against a seed of 950. The\naverage number of voxels included in each ROI of the Craddock\natlas thus varies from almost 2×103 with ten regions to about 20\nvoxels with 840.\nFor each set number of ROIs in the Craddock atlas, we removed\nthe ROIs that were empty for any subject from all subjects. We\ndiscarded the three subjects with the most empty regions to avoid\ndiscarding too many ROIs. The resulting dataset for region-size\neffect analysis includes 242 subjects with parcellation in 10 to 691\nregions.\n∗Available at doi.org/10.17632/crx7d22pym.4\nLEMON. We included a subset of the MPI-Leipzig Mind-Brain-\nBody dataset (29, 30) to assess non-linearity test-retest reliability†.\nThe subset contains the 14 subjects (1 female, ages 20 to 35,\nreported in 5-year bins, mode 25-30) with at least three resting\nstate measurement sessions, all with the same TE (to avoid possible\neffects due to scanning parameters). We applied the same CONN\ntoolbox default, “stringent” preprocessing pipeline and chose the\nAAL 90 parcellation.\nEEG data. The EEG data for this study derives from 8 minutes of\nresting-state eyes-closed recording from 215 healthy participants\nin the Max Planck Institut Leipzig Mind-Brain-Body Dataset –\nLEMON. We report here the main features of the data referring\nto the dataset presentation papers (29, 30).\nData acquisition. The EEG data was acquired with a sampling\nfrequency of 2500 Hz using 62 channels (ActiCAP, Brain Products\nGmbH, Gilching, Germany) according to the 10-10 system with\none VEOG electrode. The total acquisition lasted for 16 minutes,\ndivided into 60 s blocks with 8 eyes closed blocks interleaved with\n8 eyes open blocks.\nData preprocessing. We downloaded the preprocessed version of\nthe dataset containing 204 subjects ‡. The preprocessing included\nbandpass filtering (1–45 Hz), downsampling to 250 Hz, removal\nof corrupted channels, and eye movement and heartbeat removal\nvia ICA. We excluded from the analysis 7 electrodes frequently\nmissing (T7, T8, Cz, F7, CP6, PO10, Fp2) and considered the 150\nsubjects with all the remaining 54 channels available.\nWe further processed the data in Python using the MNE ( 31)\nand SciPy (32) packages. At this stage, the sequences are split into\nsegments up to 60 s long. We obtained the five usual frequency\nbands (δ= [1, 4] Hz, θ= [4, 8] Hz, α= [8, 12] Hz, β= [12, 30] Hz,\nγ= [30, 44] Hz) from each segment using an IIR Butterworth filter\nof order 4 in two forward and backwards passes to minimise phase\ndistortion. We obtained the scalp voltage by down-sampling to\n1.25 times the Nyquist frequency of the upper limit of each band\nof the filtered sequence to avoid biases (see “MI estimator”).\nWe derived the scalp BLP from the filtered signal by applying\nHilbert transformation and taking block averages of the modulus\nover 125 ms windows. We obtained the shadow dataset for BLP\nsequences computing FT multivariate surrogates before applying\nthe Hilbert transformation.\nFinally, we obtained the three epochs used for RNL and\nreliability estimation by chaining together segments for a total of\n124 s. We obtained the na¨ ıve shadow dataset as a surrogate of the\nBLP sequences right before RNL estimation.\niEEG data. The iEEG dataset is derived from the open dataset\npublished by the Sleep-Wake-Epilepsy-Center (SWEC) of the\nUniversity Department of Neurology at the Inselspital Bern and\nthe Integrated Systems Laboratory of the ETH Zurich ( 33).\nThe original dataset §, contains 2656 hours of anonymised and\ncontinuous intracranial electroencephalography (iEEG) of 18\npatients with pharmaco-resistant epilepsies. All the patients gave\nwritten informed consent that their iEEG data might be used for\nresearch and teaching purposes. Given the extensive size of this\ndataset, for each subject, we selected the central 124 s of the first\nhour that was at least 45 minutes away from any seizure based on\nmetadata and manual inspection.\nData acquisition. The iEEG signals were recorded intracranially by\nstrip, grid, and depth electrodes. The recorded signal was saved\nat a rate of 512 or 1024 Hz after band-passing between 0.5 and\n120 Hz with a double-pass fourth-order Butterworth filter. An\nepileptologist visually inspected all the recordings, marked the\nonset and end of each seizure, and removed corrupted channels.\nData preprocessing. The 18 subjects have between 24 and 128\nelectrodes of undisclosed type in undisclosed locations. Also, it is\n†Available at: https://fcon 1000.projects.nitrc.org/indi/retro/MPI LEMON/downloads/download MRI.\nhtml\n‡Available at: https://fcon 1000.projects.nitrc.org/indi/retro/MPI LEMON/downloads/download EEG.\nhtml\n§Available at: http://ieeg-swez.ethz.ch/\n8 — Tani Raffaelli et al.\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted November 18, 2024. ; https://doi.org/10.1101/2024.11.17.623635doi: bioRxiv preprint \n\nPREPRINT\n993\n994\n995\n996\n997\n998\n999\n1000\n1001\n1002\n1003\n1004\n1005\n1006\n1007\n1008\n1009\n1010\n1011\n1012\n1013\n1014\n1015\n1016\n1017\n1018\n1019\n1020\n1021\n1022\n1023\n1024\n1025\n1026\n1027\n1028\n1029\n1030\n1031\n1032\n1033\n1034\n1035\n1036\n1037\n1038\n1039\n1040\n1041\n1042\n1043\n1044\n1045\n1046\n1047\n1048\n1049\n1050\n1051\n1052\n1053\n1054\n1055\n1056\n1057\n1058\n1059\n1060\n1061\n1062\n1063\n1064\n1065\n1066\n1067\n1068\n1069\n1070\n1071\n1072\n1073\n1074\n1075\n1076\n1077\n1078\n1079\n1080\n1081\n1082\n1083\n1084\n1085\n1086\n1087\n1088\n1089\n1090\n1091\n1092\n1093\n1094\n1095\n1096\n1097\n1098\n1099\n1100\n1101\n1102\n1103\n1104\n1105\n1106\n1107\n1108\n1109\n1110\n1111\n1112\n1113\n1114\n1115\n1116\nnot reported which electrodes are located in the epileptic foci. To\nimprove data homogeneity, we randomly selected 24 electrodes for\neach subject. We band-passed and down-sampled the time series\nas we did for the EEG electrode voltage and selected a segment\nof 124 s. We extracted three additional subsets of the data. In\nthe first, the sequences have the same starting time but cover\ndifferent time lengths (up to 22 minutes and 44 s) depending on\nthe frequency band to allow in each band the same number of\nsamples as for band gamma with 124 s. The other two correspond\nto windows of 124 s extracted 24 and 48 hours after the first one.\nIn case of a seizure or a total length of the recording under 48\nhours, we selected seizure-free hours, trying to keep the distance\nbetween windows as close as possible to 24 hours.\nSingle Unit Spikes data. The last dataset we include is derived\nfrom the Allen Brain Observatory – Neuropixels Visual Coding\ndataset (\n34). The dataset ¶ contains individual units isolated\nin 58 mice from six probes inserted in the visual areas. The\ndataset section relevant to this work includes spiking times from\nindividual units and metadata about unit localisation and quality.\nWe considered the 26 specimens that underwent the “functional\nconnectivity” experiment. This included 30 minutes of spontaneous\nactivity, i.e., resting-state data.\nData preprocessing. With this dataset, we aimed to explore the\neffects of averaging on non-linearity with the finest granularity.\nWe first selected units with reasonable quality measures (expected\nfraction of missed spikes ≤0.08, maximum inter-spike interval\nviolations = 0.2 ( 35)). Then, we created groups of 32 units from\nthe same brain structure. If a structure had enough units to form\nmore than one group, we determined the groups so that the units\nwere spatially separated. When the number of exceeding units is\nnot enough to form an extra group, we discarded those of lower\nquality weighting in order: ISI violations, quality if isolation from\nother units, and the expected fraction of missed spikes. We kept\nonly the 16 sessions that allowed the identification of at least nine\ngroups.\nIn this case, our data will be the average spike count per unit in\ntime intervals of different lengths. As our estimator is designed to\nwork with continuous variables, we added a slight jitter to the data\nto keep all points distinct. For every session, we derived 42 sets of\ntime series varying the time interval from 125 ms to 8 s (doubling\nthe length at each step) and taking the average spike count of 1 to\n32 units (doubling the number at each step), assessing thus the\neffect of either temporal or spatial averaging.\nACKNOWLEDGMENTS. The publication was\nsupported by ERDF-Project Brain dynamics, No.\nCZ.02.01.01/00/22\n008/0004643, the Czech Science Foundation\nprojects No. 21-32608S and No. 21-17211S.\n1. KJ Friston, Functional and effective connectivity in neuroimaging: A synthesis. Hum. Brain\nMapp. 2, 56–78 (1994).\n2. KJ Friston, Functional and Effective Connectivity: A Review. Brain Connect. 1, 13–36\n(2011).\n3. ME Raichle, et al., A default mode of brain function. Proc. Natl. Acad. Sci. United States\nAm. 98, 676–682 (2001).\n4. BB Biswal, et al., Toward discovery science of human brain function. Proc. Natl. Acad. Sci.\nUnited States Am. 107, 4734–4739 (2010).\n5. H Bakhshayesh, SP Fitzgibbon, AS Janani, TS Grummett, KJ Pope, Detecting synchrony in\nEEG: A comparative study of functional connectivity measures. Comput. Biol. Medicine 105,\n1–15 (2019) application to EEG, good grades to MI with EQB.\n6. AS Mahadevan, UA Tooley, MA Bertolero, AP Mackey, DS Bassett, Evaluating the sensitivity\nof functional connectivity measures to motion artifact in resting-state fmri data. NeuroImage\n241, 118408 (2021).\n7. Z Wang, A Alahmadi, D Zhu, T Li, Brain functional connectivity analysis using mutual\ninformation in 2015 IEEE Global Conference on Signal and Information Processing\n(GlobalSIP). (IEEE), pp. 542–546 (2015).\n8. H Chen, Y Song, X Li, A deep learning framework for identifying children with ADHD using\nan EEG-based brain network. Neurocomputing 356, 83–96 (2019).\n9. HE Wang, et al., A systematic framework for functional connectivity measures. Front.\nNeurosci. 8, 111632 (2014).\n10. J Hlinka, M Palu ˇs, M Vejmelka, D Mantini, M Corbetta, Functional connectivity in\nresting-state fmri: Is linear correlation sufficient? NeuroImage 54, 2218–2225 (2011).\n11. N Tzourio-Mazoyer, et al., Automated anatomical labeling of activations in SPM using a\nmacroscopic anatomical parcellation of the MNI MRI single-subject brain. NeuroImage 15,\n273–289 (2002).\n12. J Hlinka, D Hartman, M Vejmelka, D Novotn ´a, M Paluˇs, Non-linear dependence and\nteleconnections in climate data: Sources, relevance, nonstationarity. Clim. Dyn. 42,\n1873–1886 (2014).\n13. D Hartman, J Hlinka, Nonlinearity in stock networks. Chaos 28 (2018).\n14. C Chang, GH Glover, Time–frequency dynamics of resting-state brain connectivity\nmeasured with fmri. NeuroImage 50, 81–98 (2010).\n15. J Cabral, ML Kringelbach, G Deco, Functional connectivity dynamically evolves on multiple\ntime-scales over a static structural connectome: Models and mechanisms. NeuroImage\n160, 84–96 (2017).\n16. RC Craddock, GA James, PE Holtzheimer, XP Hu, HS Mayberg, A whole brain fMRI atlas\ngenerated via spatially constrained spectral clustering. Hum. Brain Mapp. 33, 1914–1928\n(2012).\n17. SM Motlaghian, et al., Nonlinear functional network connectivity in resting functional\nmagnetic resonance imaging data. Hum. Brain Mapp. 43, 4556–4566 (2022).\n¶Available at: https://allensdk.readthedocs.io/en/latest/visual coding neuropixels.html\n18. B Frauscher, et al., Atlas of the normal intracranial electroencephalogram:\nneurophysiological awake activity in different cortical areas. Brain 141, 1130–1144 (2018).\n19. MJ Brookes, et al., Measuring functional connectivity using MEG: Methodology and\ncomparison with fcMRI. NeuroImage 56, 1082–1104 (2011).\n20. JF Hipp, M Siegel, Bold fmri correlation reflects frequency-specific neuronal correlation.\nCurr. Biol. 25, 1368–1374 (2015).\n21. GA Darbellay, I Vajda, Estimation of the information by an adaptive partitioning of the\nobservation space. IEEE T ransactions on Inf. Theory 45, 1315–1321 (1999) equiprobable\nintervals.\n22. M Palu ˇs, Testing for nonlinearity using redundancies: quantitative and qualitative aspects.\nPhys. D 80, 186–205 (1995).\n23. JA Bonachela, H Hinrichsen, MA Mu ˜noz, Entropy estimates of small data sets. J. Phys. A:\nMath. Theor. 41, 202001 (2008).\n24. RQ Quiroga, A Kraskov, T Kreuz, P Grassberger, Performance of different synchronization\nmeasures in real data: A case study on electroencephalographic signals. Phys. Rev. E -\nStat. Physics, Plasmas, Fluids, Relat. Interdiscip. T op. 65, 14 (2002).\n25. Z Jia, Y Lin, Y Liu, Z Jiao, J Wang, Refined nonuniform embedding for coupling detection in\nmultivariate time series. Phys. Rev. E 101, 062113 (2020).\n26. D Prichard, J Theiler, Generating surrogate data for time series with several simultaneously\nmeasured variables. Phys. Rev. Lett. 73, 951 (1994).\n27. M Palu ˇs, Detecting phase synchronization in noisy systems. Phys. Lett. A 235, 341–351\n(1997).\n28. J Kopal, A Pidnebesna, D Tomeˇcek, J Tintˇera, J Hlinka, Typicality of functional connectivity\nrobustly captures motion artifacts in rs-fMRI across datasets, atlases, and preprocessing\npipelines. Hum. Brain Mapp. 41, 5325–5340 (2020).\n29. A Babayan, et al., A mind-brain-body dataset of MRI, EEG, cognition, emotion, and\nperipheral physiology in young and old adults. Sci. Data 2019 6:1 6, 1–21 (2019).\n30. N Mendes, et al., A functional connectome phenotyping dataset including cognitive state\nand personality measures. Sci. Data 2019 6:1 6, 1–19 (2019).\n31. A Gramfort, et al., MEG and EEG data analysis with MNE-Python. Front. Neurosci. 7, 1–13\n(2013).\n32. P Virtanen, et al., SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.\nNat. Methods 17, 261–272 (2020).\n33. A Burrello, L Cavigelli, K Schindler, L Benini, A Rahimi, Laelaps: An Energy-Efficient\nSeizure Detection Algorithm from Long-term Human iEEG Recordings without False Alarms\nin 2019 Design, Automation & T est in Europe Conference & Exhibition (DA TE). (IEEE), pp.\n752–757 (2019).\n34. JH Siegle, et al., Survey of spiking in the mouse visual system reveals functional hierarchy.\nNature 592, 86–92 (2021).\n35. DN Hill, SB Mehta, D Kleinfeld, Quality metrics to accompany spike sorting of extracellular\nsignals. J. Neurosci. 31, 8699–8705 (2011).\nTani Raffaelli et al. Preprint — November 17, 2024 — 9\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted November 18, 2024. ; https://doi.org/10.1101/2024.11.17.623635doi: bioRxiv preprint","source_license":"CC-BY-4.0","license_restricted":false}