Stimulus-specific information is represented as local activity patterns across the brain

Modern neuroimaging represents three-dimensional brain activity, which varies across brain regions. It remains unknown whether activity of different brain regions has similar spatial organization to reflect similar cognitive processes. We developed a rotational cross-correlation method allowing a straightforward analysis of spatial activity patterns distributed across the brain in stimulation specific contrast images. Results of this method were verified using several statistical approaches on real and simulated random datasets. We found, for example, that the seed patterns in the fusiform face area were robustly correlated to brain regions involved in face-specific representations. These regions differed from the non-specific visual network meaning that activity structure in the brain is locally preserved in stimulus-specific regions. Our findings indicate spatially correlated perceptual representations in cerebral activity and suggest that the 3D coding of the processed information is organized using locally preserved activity patterns across the brain. More generally, our results demonstrate that information is represented and shared in the local spatial configurations of brain activity.


Introduction
The brain is a three-dimensional structure which may be considered as an ensemble of small volumes called voxels, with the activity of each voxel representing a summary level of neuroglial activity ( Shulman et al., 2004 ;Strelnikov, 2010 ). Differences of activity between neuroglial populations due to the combination of excitatory and inhibitory connections form specific spatial patterns, which are part of neural coding ( Yamane et al., 2008 ). For example, spatial pattern formation is observed in the primary visual cortex where some cortical ✩ All authors reviewed the manuscript and discussed the results.
cases it is possible to distinguish complex sounds on the basis of the auditory activity patterns ( Paquette et al., 2018 ), and one can also test local differences between stimulus-classes using the searchlight MVPA method ( Etzel et al., 2013 ). In this article, we hypothesize that spatial differences of activity between the neighboring neuroglial populations, which form spatial patterns, constitute an important aspect of information coding not only at the sensory stages but are organized in a network of similar activity patterns across the brain at different hierarchical levels. The similarity of spatial patterns across the brain may be due to several reasons: spatial coding may be transmitted between brain areas; non-adjacent brain areas may receive similar spatial information from the same lower-level structures; or similar spatial patterns may represent repeating computational strategies employed by the brain at different scales/levels of computational hierarchy.
To test this hypothesis, instead of searching for differences between two conditions as in MVPA studies, one can explore the images which already reflect differences between stimulus and baseline, e.g. as provided by contrast images in SPM (Statistical Parametric Mapping) analysis. In these images with established differences between stimulus and baseline, one can explore whether stimulus-specific activity patterns are similar between brain regions. The resemblance of activity organization between regions would suggest similar spatial relationships between neighboring voxels which reflect specific information processing (possibly, some sort of spatial coding of information).
More specifically, such estimation of similar activity patterns across the brain would be to correlate a certain multivoxel spatial pattern in a region of interest with other spatial patterns in the brain; in other words, to apply a cross-correlation analysis using all possible 3D rotations of the seed pattern to adapt to various local orientations of the cortex. The cross-correlation technique provides a unique correlation value per voxel, which summarizes the resemblance of spatial activity in the vicinity of this voxel with the seed pattern. In addition to the increasingly popular temporal connectivity analysis, our approach could be considered as the first analysis of spatial connectivity in functional neuroimaging, which searches for the existence of similar spatially encoded information (in the form of 3D time-averaged activity patterns) between brain areas.
In this article, we present an implementation of such a crosscorrelation analysis and use it to analyze datasets which span two sensory modalities. Firstly, we simulated random datasets to confirm the robustness of our analysis. Then, using a dataset of contrast images for face and non-face stimuli, we showed that face-specific representations are shared so as to form a network of similar spatial patterns across the brain. Furthermore, we used a second dataset from a speech processing task. Here, our analyses again showed that there is a network of similar spatial patterns across the brain dedicated to word processing.

Materials and methods
In this section, we start by the description of the proposed analysis pipeline. Afterwards, we describe each used dataset and how the analysis pipeline was applied to it. We make available all the C ++ , MATLAB scripts and datasets used for this article so that any further details can be checked in the scripts.

Analysis pipeline
The analysis pipeline consisted of the spatial cross-correlation analysis, processing of the results of the spatial cross-correlation analysis (space localization, shift correction, z-transform) and statistical analysis of these processed results of cross-correlation. We provide a package with the cross program and Matlab scripts for the users with detailed instructions and examples how to conduct the analysis ( https://osf.io/ze6bp/ ).

Spatial cross-correlation analysis
We built a toolbox in C ++ using ITK (Insight Toolkit ( https://itk.org/) ) allowing cross-correlation analysis using the Fourier transform of the 3D patterns. The main executable file in our toolbox (cross.exe) takes as input the name of the brain image file and the coordinates of the seed pattern (tested with Windows 10). The seed 3D pattern taken at a given location in the brain is rotated in steps of 10 degrees about the z, y and x axes, and for each rotation step the pattern is cross-correlated with the entire brain activity ( Fig. 1 A  and B). The maximal value of correlation obtained during the seed rotations at the given location is saved as a voxel value in the NIFTI file alongside with NIFTI files for the x, y and z angles of rotation corresponding to this correlation. The step of 10 degrees was chosen as the smallest reasonable one not to make the analysis excessively long. The size of the cubic seed can be varied in the command line of the toolbox; in this analysis we took a cubic seed with a side-length of 10 voxels for the cross analysis because its rotation by 10 degrees displaces the border of the cube by less than 1 voxel, thus no voxels are skipped during rotations. Some combinations of the 180°rotations led to the reflections of the seed. Only the seed is rotated while the whole brain always remains in its initial position (defined by the MNI template in our analysis). Given 18 steps of rotations around each axis to obtain all possible non-redundant rotations, the total number of configurations was 18 × 18 × 18 = 5832. There is no physiological reason to require that exactly 1000 voxels (10 × 10 × 10) of the pattern should be correlated with voxels in other areas; in our analysis, we preserved all correlations when the centre of rotation of the pattern (which corresponds to the centre of the seed pattern) lies within the brain. Due to the rotation, the seeds (which are cubes) span a spherical volume.
Since both the rotations and the cross-correlation analysis are computationally intensive, it is important to accelerate the analysis using Fourier transform. The analysis of one subject took about 43 min on the PC optimized for performance (available in Windows options). The whole analysis on 16 subjects took about 4 hours. Importantly, we did not smooth the contrast images, which were used for cross-correlation. Smoothing was used only as part of preprocessing of the images to improve spatial correspondence and signal to noise ratio ( Wang et al., 2005 ) before making contrasts. In our analysis pipeline, we took the seed at the site of the significant activity peak in the classical group level analysis reported in the used datasets, thus we followed the same pre-processing steps as the published results.
When the activity peak is at the periphery of the brain, it is not useful to take the center of the pattern at the peak of the activity because in this case a large part of the pattern would be located outside the brain. When peaks are peripheral, it is better to center the pattern deeper into the brain with respect to the peak of activity so that the periphery of the pattern includes the peak (see Fig. 1

C and D).
We verified whether the considered pattern of activity with a dimension ( "radius of rotation ") of 5 voxels contains the peak of the activity. Given that the center of the pattern is the point q , the distance between that peak p and the point q should be inferior to the chosen radius for the pattern. In order to calculate that distance, we used the Euclidean distance d(p,q) between the peak of activity p and the point q in the chosen coordinates (1): The obtained distance was 4.12 voxels for the 1 st fMRI dataset (WH dataset) and 3.74 voxels for the 2 nd one (ST dataset), so inferior to the distance of 5 voxels.
At the most statistically significant peak in the group analysis, which should correspond to the seed location, we verified in each subject if the correlations at these coordinates are not different from 1. Furthermore, we verified the angles of rotations at the most significant locations at the group level and found that these maximum cor- The seed pattern is turned then compared with target patterns. (B) Activity in the seed pattern is crosscorrelated with the entire brain activity. (C) An example of a peak of activity where the pattern is centered deeper in the brain with respect to the peak. (D) An example of a pattern that has a large part outside of the brain (indicated by the word 'Out' with arrow) -the situation to be avoided by moving the pattern deeper. In both C and D the square represents the pattern and the dot inside the location of the peak of activity. MRI-Cron software was used to display the figures ( https://www.nitrc.org/projects/mricron ). relations are not defined by the geometry of the cortex in the target region but fluctuate with a rather high inter-subject variability. For example, the rotation with a maximum correlation in the orbitofrontal cortex corresponds approximately to 60-90°of rotation around the X axis (SEM = 13.8). Our analysis produces images containing the information on the rotational angles corresponding to the maximal correlation.

Processing the results of spatial cross-correlation
We put the obtained result (crosscorr_final.nii) (see supplementary information) into the same MNI (Montreal Neurological Institute) space as the original image using the 3D matrix in the initial SPM contrast image (3 mm isotropic voxels). In addition, given the presence of a small shift between the cross-correlation image and the initial contrast, we used a translation of the obtained result in MATLAB (v 2017a). The obtained file was renamed (crosscorr_final_co_trans.nii). In order to verify if the latter has correctly been realigned to the initial contrast, we superimposed both images in SPM and checked that the peak in the cross-correlation image corresponds to the chosen coordinates in the initial contrast. Furthermore, we z-transformed the data after subtracting the mean to approximate the normal distribution. The resulting were then renamed (Ztransf_Crosscorr_final_co_transl.nii). Thereafter, we created randomized brain images for every subject on the basis of the ztransformed images in order to use them as a reference for the statistical analysis. Indeed, in order to be considered statistically significant, a given data should be significantly different from a random distribution of the same type of data. The randomization method used was random permutation of the coordinates within the images of z-values without repeating elements. The null hypothesis is that correlations (represented by z values after the z-transform) are randomly distributed in the brain for each subject. According to this null hypothesis, there should be no significant clusters at the group level. When significant clusters in the group analysis are found, it means the rejection of the null hypothesis within the cluster.

Statistical analysis of the results of spatial cross-correlation
SPM12 toolbox ( https://www.fil.ion.ucl.ac.uk/spm/software/spm12/ ) was used to analyze the data; statistical analysis was applied to the maximal values of normalized cross-correlations per voxel after ztransform. We used the two-sample t-test to compare the obtained z-transformed images and the randomized images (5% FWE (Family-Wise Error) inference with P = 0.001 cluster-forming threshold, extent threshold = 0). These p-values are corrected for the multiple dependent comparisons and are based on the probability of obtaining clusters with k, or more, voxels above a threshold on the basis of Random field theory with Euler characteristics ( Nichols and Hayasaka, 2003 ) (here, k is the number of voxels passing the threshold).
Correlation coefficients at the significant peaks were assessed using bootstrap with bias corrected accelerated confidence intervals (resampling 10000 times with a bias corrected and accelerated (BCa) percentile algorithm for confidence intervals ( Carpenter and Bithell, 2000 ). A difference is judged to be statistically significant if confidence intervals do not overlap.
In order to ensure that there were no spurious correlations, we also used the two-sample t-test on the z-transformed images (see details above) to compare the visual and auditory seed regions from the same initial contrast image.
In addition, for both WH and ST datasets, given the distribution of the resulting data, which may not always respect the normal distribution, we also replicated the same analysis with a nonparametric approach using the SnPM toolbox ( https://warwick.ac. uk/fac/sci/statistics/staff/academic-research/nichols/software/snpm ) (two sample t-test, number of permutations = 256, variance smoothing = [8, 8,8]). In this case, we did not apply the z-transformation to the images; we only put them into the same MNI space as the original contrast and performed the above-mentioned translation for shift correction. This analysis is included in the supporting information part (supplementary information).

fMRI datasets and their analyses
Firstly, we created simulated random images of brain activity as control datasets to verify whether or not statistically significant spatial patterns could be found in random brain activity and whether the shape of the brain could influence the results.
In order to analyze functional pattern similarities in the brain during a given task, we chose two fMRI datasets, each of them coming from studies that explored different questions and sensorial modalities. In the first fMRI dataset ( Wakeman and Henson, 2015 ), authors used visual stimuli (facial and non-facial ones) on human subjects and explored the cortical response to face and non-face stimuli presentation. The second fMRI auditory dataset comes from a work where the authors ( Strelnikov et al., 2011 ) used complex auditory stimuli (words) to study cortical responses. We analyzed these fMRI datasets to explore the possibility that both visual and auditory information are represented and shared in spatial 3D patterns of brain activity.
See further details on each dataset below.

Simulated random datasets (RANDOM datasets)
In order to perform a statistical control analysis for our crosscorrelation method, we generated 100 random datasets, which simulated brain NIFTI images having the same dimensions (sizes: x = 53, y = 63 and z = 46 voxels), shape and 3D smoothing as the 1 st fMRI dataset ( WH, 2015 , see below). These images preserved the shape of the brain in the MNI space; randomization consisted in the filling of brain shapes with random values from the uniform distribution in the range of the real values from the datasets, afterwards smoothing was applied. After the smoothing, the inclusive spatial mask was applied once again to restore the brain shape. We analyzed cross-correlation with occipital and temporal seeds as in the 1 st fMRI dataset in these random simulated images to test the possible effects of brain shape and smoothing on the results of cross-correlation (occipital coordinate: x = 16, y = 13 and z = 15 voxels ( ⇔ x = 33, y = − 76 and z = − 8 mm); temporal coordinate: x = 10, y = 32 and z = 21 voxels ( ⇔x = 51, y = − 19, z = 10 mm)). This direct effect of smoothing was statistically tested on the group of 16 images. In addition, we took the difference between random smoothed images to simulate 16 random "stimulation vs. baseline " contrasts; this simulation is close to our approach on the real contrasted images.
Afterwards, we performed an analysis on 100 randomly formed different groups of images using the occipital seed (occipital coordinate: x = 16, y = 13 and z = 15 voxels ( ⇔ x = 33, y = − 76 and z = − 8 mm). We took the difference between random smoothed images to simulate 16 random "stimulation vs. baseline " contrasts per group. To exclude autocorrelations, double the volume of the seed was masked out. After the same analysis as in the real datasets (see Analysis pipeline), we assessed the fraction of false positive voxels at p < 0.05, uncorrected (t > 1.7) and at p < 0.05, FWE corrected (t > 6.27).
To estimate the fraction of true positives, in 5 random groups, we inserted the seed pattern at different coordinates (3 inserted patterns per group, one rotated using ITK Euler3DTransform and two non-rotated). As we expected no false negatives with this approach, we masked one non-rotated pattern per group with different levels of noise. 100% noise was defined as random values sampled from a uniform distribution in the range of the real values from the pattern, which were averaged with pattern values. In the case of 50% noise, random values were sampled from 50% of the range of the real values and similarly for other noise levels.

fMRI dataset 1 ( WH 2015 )
The freely available data were chosen from a work of Wakeman and Henson ( Wakeman and Henson, 2015 ) ( ftp://ftp.mrccbu.cam.ac.uk/personal/rik.henson/wakemandg_hensonrn/ ). The MATLAB scripts attached to the original dataset were used for preprocessing and statistical analysis (using the MATLAB SPM toolbox), resulting in contrast images in the MNI space with 3mm isotropic voxels smoothed by a 3D 8mm isotropic Gaussian kernel before creating the contrasts. In this study, the subjects (n = 16) were presented gray scale images of familiar and unfamiliar faces, and faces scrambled by a 2-D Fourier transform. Original and scrambled faces were cropped using a mask built on the basis of a combination of familiar and unfamiliar faces.
This dataset is composed of 2 contrasts: con_0005.nii and con_0006.nii. They are renamed in the present manuscript as non-face specific and face-specific contrasts, respectively. In each of them, the stimulus condition has been compared to the baseline in the study of Wakemann and Henson, 2015 . The face-specific contrast corresponds to sessions where the visual face stimuli were used (familiar and unfamiliar faces). The non-face specific contrast corresponds to sessions where the control visual stimuli were used (faces scrambled by a 2-D Fourier transform).
Given the nature of the stimuli used by Wakeman and Henson, we chose two regions of interests for the seeds in our analysis: an occipital and a temporal one. In the occipital region, we considered the coordinates corresponding to the peak of the brain activity in the contrasts. In the temporal seed, as a control, we choose the coordinates of the primary auditory area which should not show activity during visual stimuli presentation. The occipital seed was at the peak of the activity ( Fig. 1 C and D) in the two contrasts (x = 33, y = -76, z = -8 mm) in the Fusiform Face Area (FFA), ( ⇔ x = 16, y = 13, z = 15 voxels in matrix coordinates found from SPM Image view). The control temporal seed had coordinates x = 51, y = -19, z = 10 mm ( ⇔ x = 10, y = 32, z = 21 voxels in matrix coordinates, corresponding to the middle of the primary auditory area ( Morosan et al., 2001 )). The corresponding coordinates' areas were verified in the SPM Anatomy (Version 2.2b) and xjView ( http://www.alivelearn.net/xjview ) toolboxes.
To ensure the robustness of the cross-correlation analysis, we employed a seed-reversal approach. In this approach, the peaks of a given correlation-analysis were used as a new seed (in the contralateral FFA). If the technique is robust, this seed-reversal should lead us back to the initial seed.

fMRI dataset 2 (ST dataset)
The second fMRI dataset (indicated in the current work as fMRI auditory) was for auditory processing and originated from our work on speech processing done on 15 subjects ( Strelnikov et al., 2011 ). Subjects lay in the scanner with eyes closed and listened to dissyllabic words. One out of 13 stimuli was randomly repeated, and the subjects were instructed to press the button when they heard a repetition. The dataset comprises the contrasts of words versus silent baseline. In the present work we considered only 14 subjects. The 5 th subject has been discarded from the analysis because we suspected an anomaly regarding the shape of the brain of this subject even after spatial normalization to a standard MNI space. To verify whether the exclusion of this subject influences the general statistical results, we used one sample t-test for this dataset with and without this subject and looked for statistical differences. With this analysis, we did not find salient differences, especially regarding the peaks of activity.
Given the auditory nature of stimuli used by ( Strelnikov et al., 2011 ), we considered the temporal region as a region of interest. We chose a seed (a cube of voxels) with coordinates near the peak of activity in the L STG (Left Superior Temporal Gyrus, in the vicinity of the primary auditory area): x = 10, y = 35 and z = 19 voxels ( ⇔ x = -51, y = -9 and z = 3 mm). We used the contrasts Words_(n).nii (n denotes the number corresponding to each subject). In addition, we performed a supplemen- tary cross-correlation analysis putting the seed at the contralateral to the previous analysis site (R STG, x = 42, y = 33 and z = 15 voxels ⇔ x = 45, y = -15 and z = -9 mm).

Simulated random datasets (RANDOM datasets)
To check the role of brain shape and smoothing, we performed the cross-correlation analysis with random smoothed values within the brain mask, which preserved the shape and external borders of the brain. In the smoothed random images, the only significant effect was observed for the coordinates of the initial occipital and temporal seeds. Similarly, in the random groups with smoothed contrasted random images the only significant effect was observed for the coordinates of the initial seeds. Thus, no significant effects outside the seed regions were found in simulated random datasets. At the level of p < 0.05, uncorrected, the false positive rate in 100 random groups was 4.4 ± 0.4 % (mean ± SD) and at the level of p < 0.05, FWE corrected (used in this study), the false positive rate was 0.024 ± 0.0023 % (mean ± SD). Without masking with noise, all of the inserted non-rotated and rotated patterns were detected by the statistical analysis at p < 0.05, FWE corrected. False negatives at this threshold were observed only for more than 100% noise level (see Figure S1 in Supplementary information).

Cross-correlation analysis in the occipital area
In order to test the specificity of our method we took a seed pattern in the fusiform face area of the face-specific contrast, coordinates: x = 33, y = − 76, z = − 8 mm, equivalent to x = 16, y = 13, z = 15 voxels (here and in other places, we also indicated the voxel coordinates of the seed because they are used in the cross correlation toolbox). As expected, the most significant peak in the group analysis was found at the location which corresponds to the coordinates of the seed pattern ( Table 1 and Fig. 4 ) (a 1-2 voxel displacement is sometimes observed at the group level analysis, though exactly the same seed coordinates are found at the individual level).
The rest of the results section will omit the predictable correlations of 1 found at the seed location. The other significant peaks were found in the contralateral FFA, the right cerebellum, and the left orbitofrontal region ( Fig. 2 A and B).
To verify the consistency of the analysis, we used the contralateral peak from the previous analysis as the seed pattern. If the technique is robust this reverse-seed analysis should be able to detect the initial seed. In fact, the most significant correlation peaks of this reverse-seed analysis were indeed found at the original seed and its vicinity ( Fig. 2 C;  Fig.3 ; Table 2 ). We further verified that this reverse-seed correlation holds for each subject (mean = 0.62, bootstrapped CI = [0.55, 0.66]). Thus, a pattern in the opposite hemisphere can be detected whatever the side of the seed pattern used for the analysis. Next, to verify, as a control, if the method provides other results for non-face-specific data, we analyzed the non-face-specific contrast with the seed pattern at the same coordinates as for the above described facespecific contrast (x = 33, y = − 76 and z = − 8 mm ⇔ x = 16, y = 13, z = 15 voxels)(supplementary information, Table S1).
In addition to the seed coordinates, we found significant peaks in the occipital region (ex: x = 24, y = -88, z = -5 mm; p < 0.001, FWE correction). These peaks did not coincide with those found in the facespecific contrast.
Thus, significant correlations with a seed pattern can be found across the hemispheres both for face-specific and non-face-specific contrasts but at different locations. Furthermore, values of the cross-correlations per subject were above 0.5 and the mean correlation values for non-seed areas were in the range 0.6-0.8 ( Fig. 3 and Fig. 4 ).
Likewise, we noticed that at the seed coordinates the values of the cross-correlation are 1 or 0.999 for all the subjects. This indicates the high accuracy of our method.
Non-parametric analysis with SnPM toolbox gave similar results indicating the robustness of our analysis (see supplementary information, Tables S5 to S9).
Considering the rotations, the seed pattern at the coordinates of the maximum correlation had high inter-subject variability of rotation angles, though the average was close to the tangential to the cortex surface (for example, a seed pattern in the occipital region was rotated on average by about 60-90 degrees, around the x axis for a correlation peak in the orbitofrontal region).

Cross-correlation analysis with auditory activity
To verify the specificity of our method with respect to sensory modalities, we used the seed patterns defined in the auditory cortex for crosscorrelation with both face-specific and non-specific contrasts (at the coordinates x = 10, y = 32, z = 21 voxels ⇔ x = 51, y = − 19, z = 10 mm). Here, the seed pattern is outside the stimulus specific areas. The most significant correlation peaks were found at the seed pattern and its vicinity of the auditory cortex ( Fig. 2 D; supplementary information, Table  S2). This suggests that our technique is stimulus specific.
Conversely, we also compared the images resulting from a seed pattern in the fusiform face area with images resulting from the seed pattern in the auditory region (which served in this case as a non-specific statistical baseline). The results of this analysis (supplementary information, Table S3) were similar to those obtained using randomized images as the statistical baseline (see Methods for details).

Comparison between face-specific and non face-specific contrasts
As a supplementary control analysis, we compared the z-transformed Face-specific and Non-face-specific correlations using paired t-test and found significant differences (Table S4) in the primary visual area, cuneus, orbitofrontal, superior frontal and angular gyrus. Given the same coordinates of the seed (x = 33, y = -76, z = − 8 mm), the significant differences in this test further indicate that anatomical structures do not influence the analysis. (C) Significant clusters in the FFA resulting from cross-correlation in the face-specific contrast with the seed on the contralateral FFA. Ipsi and contralateral significant correlations are shown on the 3D template map. (D) Significant clusters in the auditory area resulting from cross-correlation in the face-specific contrast. In all images, clusters are shown at the uncorrected threshold for illustrative purposes. Arrows indicate the seed. MRI-Cron software was used to display the figures ( https://www.nitrc.org/projects/mricron ).

fMRI Dataset 2 (ST dataset)
In order to verify the replicability of our approach, we conducted a cross-correlation analysis in an independent dataset putting the seed coordinates near the peak of activity in the left STG.
As expected, the most significant peak in the group analysis was found at the location, which corresponds to the coordinates of the seed pattern. In addition, significant correlations were found in areas involved in word processing, for example: L and R STG, temporo-polar areas and supramarginal gyri (Tables S10 and S12), indicating a high specificity of our approach.
Furthermore, similar to the first fMRI dataset (WH dataset), all individual values of the cross-correlation exactly at the seed coordinates were equal to one. The values of the cross-correlations per subject for significant non-seed areas had confidence intervals of 0.48-0.59 according to the bootstrap analysis at p < 0.05 ( Fig. 5 ).
Moreover, we performed an additional cross-correlation analysis taking the seed in the contralateral hemisphere to the previous analysis position (in the R STG). In this reverse-seed analysis, the most significant correlation peaks were at the original seed (R STG), its vicinity and in the L STG (Tables S11 and S13) indicating, as in WH dataset, that a pattern in the opposite hemisphere can be detected whatever the side of the seed pattern used for the analysis.

Discussion
In the present work, we propose a new approach to analyze fMRI data using spatial correlations. As a proof-of-concept, we apply our approach to publicly available fMRI datasets (see Methods), and show that our methodology allows us to find specific spatial patterns of the activity following the presentation of visual (facial, in WH dataset) and auditory stimuli (ST dataset). The validity of our statistical analysis is confirmed by the fact that it always detects the seed pattern at the right locations and with the highest significance both in real and simulated random datasets.

Face processing viewed through long-range spatial correlations
In the first dataset, we found that spatial patterns are shared with the contralateral FFA with high correlation in both group and individual analyses, indicating the conservation of spatial information regarding faces in both left and right FFA. In addition, by performing reverse-seed correlations (setting the seed in the contralateral area) we found the same results in the FFA -confirming the robustness of our approach. In our results, we found that the correlations are high (up to 0.8), suggesting that a large amount of spatial information in brain activity is shared between these areas.
The FFA is known to respond to face stimuli ( de Souza et al., 2008 ;Ghuman et al., 2014 ;Kanwisher et al., 1997 ;Posamentier and Abdi, 2003 ;Sergent et al., 1992 ), however, the spatial conservation and transmission of face-related information between brain areas involving the FFA has never been reported. The present work provides the evidence of stable spatial patterns of brain activity in response to face stimuli, which preserve information across brain structures. In addition, we also found a high correlation in the orbitofrontal cortex, which is in line with several findings which have demonstrated the existence of face-selective neurons in the orbitofrontal cortex ( Barat et al., 2018 ;Posamentier and Abdi, 2003 ;Troiani et al., 2016 ), and shown that Fig. 3. Values of the cross-correlation for each subject in the fusiform areas (Contrast: face-specific, seed in the left fusiform area). Each subject's value (circles) indicates the result of the cross-correlation at the above-mentioned coordinates. The height of the main bars corresponds to the mean value. Error bars indicate the bootstrap confidence intervals (CI; re-sampling 10000 times with a bias corrected and accelerated (BCa) percentile algorithm for confidence intervals ( Carpenter and Bithell, 2000 )). In parentheses, MNI coordinates of the areas (x,y,z) are indicated. this region is actively involved in face recognition ( Barat et al., 2018 ;Posamentier and Abdi, 2003 ;Trautmann et al., 2009 ;Troiani et al., 2016 ;Willis et al., 2014 ). Barat and collaborators (2018) showed that orbitofrontal face cells encode facial stimuli by, first, discriminating them from non-facial ones, and thereafter categorizing them according to their social and emotional dimensions. Face neurons encode these aspects despite differences in face stimuli concerning, for example, identity or head position ( Barat et al., 2018 ). In the present study, we used the dataset from a work where the presented face-stimuli have various expressions (generally, happy or neutral). The presence of the significant correlation in the orbitofrontal cortex in our findings indicates that the social and Error bars indicate the bootstrap confidence intervals (CI; re-sampling 10000 times with a bias corrected and accelerated (BCa) percentile algorithm for confidence intervals ( Carpenter and Bithell, 2000 ). In parentheses, MNI coordinates of the areas (x,y,z) are indicated. emotional dimensions that are carried by faces are represented as spatial patterns of activity. In addition, given that visual information is received by the occipital cortex before the orbital one, the existence of high correlation between these areas suggests that spatial coding may be transmitted or shared between them, e.g. from the occipital (FFA) to the orbitofrontal cortex in the case of bottom-up processing. Similarly, we also found significant cerebellar correlations for face-stimuli -an observation that is in agreement with studies which have implicated the cerebellum in facial information analysis ( Adamaszek et al., 2014 ;D'Agata et al., 2011 ;Posamentier and Abdi, 2003 ;Uono et al., 2016 ).
Thus, the obtained correlations are highly specific to facial information. They are found in the FFA, the orbitofrontal and cerebellum regions, which are known for being activated after facial stimuli presen- Fig. 4. Values of cross-correlation for each subject in the fusiform, visual and other areas. A. Face-specific contrast. B. Non-face specific contrast. For both figures, the seed was in the right fusiform area. Each subject's value (open circles) indicates the result of the crosscorrelation at the above-mentioned coordinates. The height of the main bars corresponds to the mean value. Error bars (in blue) indicate the bootstrap confidence intervals (re-sampling 10000 times with a bias corrected and accelerated (BCa) percentile algorithm for confidence intervals ( Carpenter and Bithell, 2000 )).
In parentheses, MNI coordinates of the areas (x,y,z) are indicated. tations in several studies as described above. In addition to the above agreement with literature, another argument for the stimulus specificity of our results is the fact that the spatial correlation analysis was performed on the "stimulus vs. baseline" contrasts and all the images were acquired under the same conditions and processed in the same way excluding acquisition and processing artefacts in contrast images. In addition, we also showed that using a seed in the auditory region produces high correlations in the auditory areas -further indication of the specificity of our methodology.

Complex auditory information (words) is distributed as spatial patterns of brain activity
Analyzing the second dataset, we found significant correlations (0.5-0.6) between areas involved in word processing, mostly in both L and R STG, temporo-polar areas and supramarginal gyri. These areas, involved in speech processing pathways, are known to process lexical information ( Davis, 2016 ;Davis and Gaskell, 2009 ). Phonetic forms are coded in the STG and, thereafter, lexical information is integrated with semantic representations in the anterior parts of the temporal lobe ( Davis, 2016 ). In addition, we found significant correlations in the supramarginal gyrus, which is considered an auditory-motor interface that organizes speech information into premotor and auditory representations ( Davis, 2016 ;Davis and Gaskell, 2009 ).
In addition, we found correlations in the insula and the cerebellum. These structures are involved in mapping speech information in order to organize motor responses ( Davis and Gaskell, 2009 ). The cerebellum, for example, is involved in speech processing in several studies, especially in the organization of syllables and in the pre-articulatory verbal information stages ( Ackermann, 2008 ;Moberget and Ivry, 2016 ;Strelnikov et al., 2011 ).
Our findings indicate that information regarding lexical processing and organization in different specialized areas is shared as spatial patterns of brain activity.

No meaningful spatial patterns were found in random images of brain activity (RANDOM datasets)
As a control analysis, we simulated smoothed random brain activity in order to verify whether the analysis of these simulated images would lead to significantly high correlations due to smoothing or to the brain shape. Our analysis in the Random datasets preserved brain shapes, the absence of significant results in the same statistical analysis as for the real datasets indicates the absence of the effects of cortical folding ( Tijms et al., 2012 ). Given the same coordinates of the seed, the significant differences in the contrast "face-specific vs. non facespecific " further indicate that anatomical structures do not influence the analysis. Besides, we performed cross-correlation on "stimulation vs. baseline " contrasts for each subject; as morphological aspects are the same for stimulation and baseline, they are likely to have negligible effect in the contrasts per subject. Effects of smoothing were also absent both in the direct analysis of the smoothed random images and in the analysis of their differences, which simulated "stimulation vs baseline " contrasts. The analyses of several groups of random brain images and simulations with inserted patterns statistically confirmed the robustness of the method.

General considerations of local activity patterns
Our analysis shows that correlated stimulus-specific patterns are found throughout the brain. Moreover, this suggests that there are recurring local patterns of meaningful stimulus-related activity. One interpretation of these local patterns could be that there is a partial conservation of stimulus-related information in localized regions in the cortex, and that the activity of any single unit is related to the activity of the neighboring units. The maintenance of these local patterns would also be en-ergetically favorable as compared to longer range activity patterns, and could be an efficient neuroglial code for local copies of stimulus-related information. Likewise, we suppose that brain information is spatially organized in such a way that each neuroglial population, represented by a given voxel, maintains a given level of activation according to its neighbors to maintain this pattern stable. Maintaining patterns stable is also probably related to the necessity to spend less energy in an optimal way to maintain information-representation in close vicinity instead of using the long-distance arrangements.
Spatial analysis is widely spread in electrophysiological research, e.g., to analyze the structure of cortical activity; spatial differentiation is widely used in electrophysiological studies in the form of current source density (CSD) analysis ( Aizenman et al., 1996 ;Mitzdorf, 1985 ). In the latter, spatial differentiation permits to calculate electric current flows, given the measured electric potentials. Though the BOLD signal is an indirect and complex measure of brain activity ( Raichle and Mintun, 2006 ), its correlation with electrical measurements in the brain has been established by numerous studies (e.g., ( Goloshevsky et al., 2008 ;Heeger et al., 2000 ;Kim et al., 2004 ;Liu et al., 2010 ;Rees et al., 2000 ). Given that electromagnetic and metabolic energies are highly correlated during brain activity, Friston ( Friston, 2005 ) proposed that the most suitable summary form of energy to describe brain mechanisms is free energy. In our previous research, we investigated the possibility to apply mathematical formalism, which is similar to the electrophysiological current source density (CSD) analysis ( Aizenman et al., 1996 ;Mitzdorf, 1985 ), to the BOLD signal contrast maps. In this previous work, we observed stable task-related gradients of activity at the group level ( Strelnikov and Barone, 2012 ), necessitating the existence of the stimulus-related processes which act to maintain the described spatial gradients. These oxygenation related gradients are likely to be driven by electrical gradients, which encode the flow of local information. The results from our spatial-correlation analysis are in favor of this interpretation. We compared fMRI, EEG and MEG spatial differential activity during different tasks to see what amount of fMRI differential activity corresponds to the electromagnetic differential activity ( Strelnikov and Barone, 2014 ). Distributed source reconstruction was used to obtain 3-dimensional models of electric and magnetic activity in EEG and MEG prior to spatial differentiation. Using independent datasets with the same stimuli, we demonstrated that the mean spatial overlap of the fMRI differential activity with EEG and MEG may be about 80%. In addition, about 93% of divergence (spatial sources) in fMRI corresponded to the EEG and MEG divergence. Furthermore, previous studies have shown that glutamate related excitatory synaptic transmission accounts for about 70% of total metabolic energy turnover ( Shulman et al., 2004 ) in the brain, whereas GABAergic processes account for only about 15% of total energy turnover by neurons and glia ( Patel et al., 2005 ) -thus suggesting that the stable local BOLD patterns identified in our analysis may mostly result from local excitatory interactions between neurons.
Since brain activity is dynamic in time, it is important to understand how average BOLD contrasts can lead to local patterns of activity. To account for the stability in time of the observed activity propagation and resulting spatial patterns, one can suggest that the repetitive stimulusrelated patterns of activity propagation in the brain form the average patterns reflected in the BOLD signal contrasted images. Informationrelated activity patterns are closely related to the spatial distribution of stable activity flows in the brain. Activity (energy) flows in the brain are generally defined as coherent spatial and temporal changes in the energy turnover of neuroglial units related to information treatment ( Strelnikov, 2010 ); these flows are the result of the stimulus-driven transformations of metabolic energy that propagate in certain directions along the cellular structures (axons, dendrites, synapses, etc.) in neuroglial networks. When a neural signal reaches a neuroglial population, it increases metabolic energy (i.e., activity) in this neuroglial population causing an abrupt increase of its activity. The above described link between the directions of gradients and energy flows was confirmed by the finding that the sources, from which energy/activity flows spread in the cortex, were in the occipital cortex during face processing and in the superior temporal cortex during auditory word processing ( Strelnikov and Barone, 2012 ). It was recently suggested that patterns of brain activity may be related to integrative cognitive processes ( Kuzma, 2019 ).
Future studies using spatial cross-correlation analysis combined with MVPA could be interesting to explore whether decoding methods are more efficient in the regions where correlated spatial patterns are observed. The functional spatial connectivity approach may complement the data on intrinsic connectivity of homotopic brain areas ( Joliot et al., 2015 ). A fundamentally important question is whether during multisensory interactions such as face-voice processing ( Ho et al., 2015 ) patterns of brain activity carry the similar spatial information across modalities, these studies can further develop the models of such phenomena as McGurk effect ( Olasagasti et al., 2015 ). Moreover, in pathologies such as multiple sclerosis, epilepsy and Alzheimer disease, spatial cross-correlations could be a useful methodology for testing how patterns of structure and activity are distributed throughout the cortex ( Payoux et al., 2015 ) in addition to identifying how these patterns of activity may change in early stages of these diseases, and also during the follow up of the patients, leading to a powerful new diagnosis tool complementing the existing techniques ( Belleville et al., 2014 ;Demonet, 2017 ;Frisoni et al., 2017 ) .

Perspectives
In the present study, we considered only contrasts originating from stimulus-specific conditions, which represent the average difference between stimulation and baseline. Therefore, the nature of analyzed images excludes the temporal component from the analysis. As a perspective, we think that the study in the spatiotemporal domain may provide a more complete picture, which is for the moment beyond the scope of the present work.
Likewise, we think that our work has both fundamental and practical importance for future research. By providing strong evidence that information in the brain is organized in a system of the local 3D patterns of brain activity, our approach paves the way for future studies, especially for higher efficiency of activity-based decoding approaches. In addition, applications of spatial cross-correlation in neurological pathologies, particularly neurodegenerative ones, will be useful to study how spatial patterns of activity vary in brain pathology, leading to new advances in the field.
Appendix. Supplementary data: Data generated and analyzed during the current study as well as the cross program and the MATLAB scripts are provided as a source data file. We also provide a package for the users with detailed instructions and examples how to conduct the analysis. Download link: https://osf.io/ze6bp/?view_only = 84670d16805 c4886b37d3ed69f441760 or https://osf.io/ze6bp/

Declaration of Competing Interest
The authors declare no competing interests.