Bayesian reduced rank regression models generalizable neural fingerprints that differentiate between individuals in magnetoencephalography data

Recent magnetoencephalography (MEG) studies have reported that functional connectivity (FC) and power spectra can be used as neural fingerprints in differentiating individuals. Such studies have mainly used correlations between measurement sessions to distinguish individuals from each other. However, it has remained unclear whether such correlations might reflect a more generalizable principle of individually distinctive brain patterns. Here, we evaluated a machine‐learning based approach, termed latent‐noise Bayesian reduced rank regression (BRRR) as a means of modelling individual differences in the resting‐state MEG data of the Human Connectome Project (HCP), using FC and power spectra as neural features. First, we verified that BRRR could model and reproduce the differences between metrics that correlation‐based fingerprinting yields. We trained BRRR models to distinguish individuals based on data from one measurement and used the models to identify subsequent measurement sessions of those same individuals. The best performing BRRR models, using only 20 spatiospectral components, were able to identify subjects across measurement sessions with over 90% accuracy, approaching the highest correlation‐based accuracies. Using cross‐validation, we then determined whether that BRRR model could generalize to unseen subjects, successfully classifying the measurement sessions of novel individuals with over 80% accuracy. The results demonstrate that individual neurofunctional differences can be reliably extracted from MEG data with a low‐dimensional predictive model and that the model is able to classify novel subjects.


| INTRODUCTION
Every brain is shaped by a unique combination of biological factors and environmental influences.Therefore, we may expect unique patterns in the signals generated by each individual brain.Indeed, several studies have reported that it is possible to differentiate individuals based on their functional neuroimaging data (e.g., da Silva Castanheira et al., 2021;Finn et al., 2015).So far, such studies have mainly used summary statistics, such as functional connectivity (FC) or power spectra, to compute correlations between different measurement sessions of the same participants.Even though such correlations can be used as neural fingerprints to identify individuals within a specific subject group, it remains uncertain which factors contribute to the identifiability, and whether there are more fundamental principles of "uniqueness" that are generalizable beyond a specific subject group.Accessing the neuroimaging features that contribute to successful neural fingerprinting could importantly contribute to understanding interindividual neurofunctional dis/similarity in both health and disease.
In the present study, we utilize latent-noise Bayesian reduced rank regression (BRRR) for building a computational model to assess the differentiation between individuals based on different features of neural signals (FC, power spectra).Latent-noise BRRR assumes that the noise and the covariates affect the target variables (here FC and power spectra) through the same latent space, and it is, therefore, well suited for modelling weak effects in high-dimensional brain imaging data where the chosen covariates (e.g., age, a gene) explain only a small part of the variance in the target variables (Gillberg et al., 2016).The reduced rank restriction of the model allows for determining a low-dimensional representation of the target variables based on the chosen covariates.Thus, latent-noise BRRR can be used to simultaneously reduce the dimensionality of data and learn latent components that maximally differentiate between subject groups defined by the covariates.Here, we apply BRRR in modelling individual differences, so each subject is defined as their own group.
BRRR offers several benefits that are not available with model-free approaches, such as Pearson correlations.Since BRRR is a model, it can be trained on one set of subjects and applied to another set of subjects.This could help to assess the consistency of neurofunctional interindividual differences across datasets.Even within a single dataset, BRRR can be used to interrogate the differences underlying within-subject correlations of certain FC metrics, by examining the relationship between latent components and different metrics and frequency bands.BRRR further allows for multiple variables to be incorporated into the model, for example, one may combine FC data with demographic, cognitive and/or clinical data.
In a previous study, BRRR was used to extract latent components from magnetoencephalography (MEG) power spectra based on sibling relations: the extracted components successfully differentiated families from each other (Leppäaho et al., 2019).In that study, the authors chose power spectra as the measure of interest, given its straightforward computation and wide usage.The results of Leppäaho et al. (2019) demonstrated that BRRR works with sensor-level power spectra.In this study, we extend the analysis to source-level data and to FC metrics.Use of different FC metrics has gained increasing interest both in basic research and in clinical settings, for example, in predicting cognitive decline (e.g., Pusil et al., 2019;Vecchio et al., 2018).Here, our aim was to determine whether BRRR can model previously reported differences between FC metrics and frequency bands, what amount of model complexity (i.e., number of latent components) is required to capture the individual uniqueness of the metrics, whether the model is stable and whether it can model novel individuals not included in the training set.This study will serve as a solid foundation for future work in modelling differences between individuals.
FC metrics measure synchronization between pairs of signals.Since there are various ways to quantify signal properties, numerous FC metrics have been proposed (Bastos & Schoffelen, 2016;Colclough et al., 2016;O'Neill et al., 2015).In the present work, we focused on static FC metrics, as they are more commonly used than dynamic FC metrics and thus allow comparison of the results to existing literature.Generally, FC metrics quantify connectivity based on coherence between the phases of the signals (phase coupling) or based on the correlation between the amplitudes of the signals (amplitude coupling).Phase and amplitude coupling have been suggested to be the two primary coupling modes between neuronal oscillations, and they may reflect partly distinct neural mechanisms (Engel et al., 2013;Siems & Siegel, 2020).
In MEG, activity of a neuronal source is detected by several close-by sensors, and each sensor detects signals from several neuronal sources.Thus, there are correlations between sensor signals, which are not caused by real neural interactions.This effect called field spread can be mitigated, although not completely eliminated, by source reconstruction (Schoffelen & Gross, 2009).Because field spread is instantaneous, approaches for eliminating the artificial interactions by discarding interactions with zero phase lag have been proposed (Baselice et al., 2019;Brookes et al., 2012;Colclough et al., 2015;Hipp et al., 2012;Nolte et al., 2004;Stam et al., 2007;Vinck et al., 2011).However, this is done at the expense of losing real zero-lag interactions.To explore its effect, we here compared methods with and without field spread correction.We trained both sensor and source space models; the sensor space results can be found in the supplemental material.
In the present work, latent-noise BRRR was applied to MEG resting-state data from Human Connectome Project (HCP).Our first step was to verify that BRRR could model the differences between metrics that have been previously reported for correlation-based fingerprinting.Thus, we reproduced the correlation-based fingerprinting results for reference.Then, various BRRR models were trained to differentiate between individuals based on power spectra or FC from one measurement session.The latent components of the different models were used to identify subsequent measurement sessions of the individuals, and the identification accuracy of a given model was compared with other models.The latent spaces of the different models were compared by computing, for each model, the distance matrix among all participants.Our second objective was to determine whether the BRRR model could further classify the measurement sessions of previously unseen subjects, based on the metrics that were most successful in our first analysis.

| HCP data
The resting-state MEG data of the Human Connectome Project (HCP) (Larson-Prior et al., 2013;Van Essen et al., 2013) were collected using a whole head Magnes 3600 (4D Neuroimaging, San Diego, CA, USA) system with 248 magnetometer channels and 23 reference channels.The MEG resting-state data included 89 subjects aged between 22 and 35 years.Most of the subjects had either a monozygotic or a dizygotic twin in the dataset (33 sibling pairs in total).The resting-state data of each subject consisted of three consecutive resting-state recordings that were around 6 min in length each.The participants had their eyes open during the recordings.The MEG data had been passed through the HCP preprocessing pipeline with removal of poor-quality channels and segments, bandpass filtering (1.3-150 Hz), notch filtering to remove power line noise (at 59-61 Hz and 119-121 Hz), and removal of cardiac and eye-movement artefacts based on independent component analysis (ICA).

| MEG data processing
The code used for the analysis is available on the GitHub page of the Imaging Language Group: https://github.com/AaltoImagingLanguage/Haakana2023.

| Source reconstruction
The FieldTrip toolbox (Oostenveld et al., 2011) was used to perform source reconstruction.The HCP dataset contains three-dimensional source models of all the subjects, based on the individual structural MRIs, defined on a regular grid of different resolutions in normalized MNI space.An 8 mm grid was chosen, and source reconstruction was done using linearly-constrained minimumvariance (LCMV) beamforming (Van Veen et al., 1997).The source timeseries were clustered using principal component analysis (PCA) into 116 regions of interest (ROIs) defined by the automated anatomical labelling (AAL) atlas (Tzourio-Mazoyer et al., 2002).

| Power spectra
Power spectral density (PSD) was computed for each region of interest (ROI) using Welch's method.This was done using the default settings of Matlab R2022a: the signals are divided into maximum of eight segments with 50% overlap, and a Hamming window is applied to the segments.The PSD was split into 21 frequency bands within which a mean value was computed.The frequency bands were defined similarly to Leppäaho et al. (2019), such that the first band was 1-3 Hz, and the subsequent bands were linearly widened up to a band of 81.8-87.8Hz.

| Functional connectivity
The source-space data was filtered into delta (1 -4 Hz), theta (4 -8 Hz), alpha (8 À 13 Hz), beta (13 À 30 Hz), and gamma (30 À 50 Hz) bands.Similar frequency bands have been used in prior studies of FC fingerprinting (da Silva Castanheira et al., 2021;Sareen et al., 2021), allowing for comparison of the results.FC was computed between ROIs in the five frequency bands.Several common static FC metrics were used (Table 1).Hilbert transform was used to obtain the analytic signals, from which the amplitude envelope and the instantaneous phase were extracted to compute the amplitude metrics and phase metrics, respectively.

| Latent-noise BRRR
Latent-noise BRRR is defined as where Y NÂP contains the P-dimensional MEG data with N observations (the experimental subjects), X NÂM contains covariates for M predictors, Ψ MÂK contains lowdimensional regression coefficients for the K-dimensional latent space, Ω NÂK contains unknown latent factors (i.e., noise in the latent space), Γ KÂP mediates the effects of both covariates and noise on the target variables (i.e., it is a projection of the latent space to the observational space) and E NÂP is independent noise in the observational space.Latent-noise BRRR can be used to learn a K-dimensional representation of the observed data Y that is maximally informative with regard to the selected covariates X, for example, age groups of the subjects or familial relations between the subjects.Latent-noise BRRR is designed especially for modelling weak effects, accounting for only a small part of the variance, in highdimensional data (Gillberg et al., 2016).In the independent-noise formulation of BRRR, the signal and the noise are assumed to be independent, which can result in weak relationships being hidden by the noise.Therefore, in latent-noise BRRR the signal of interest X and the noise Ω are mediated through the same latent space, allowing for the signal model to borrow strength from the noise model (Gillberg et al., 2016).
To solve the inherent rotational unidentifiability of the reduced rank regression, shrinkage priors are set for the coefficients Ψ and Ω as well as the projection matrix Γ, enforcing the components with the strongest effects to the first columns/rows (Gillberg et al., 2016).We followed Leppäaho et al. (2019) in defining the priors.The rows k of Γ and the columns k of Ψ and Ω were given the priors where τ k is defined as and where N is the normal distribution and G is the gamma distribution of shape and rate.This procedure enforces higher shrinkage for higher k since the expected values for the gamma distribution are positive.
T A B L E 1 FC metrics used in the study.The metrics differ with respect to the coupling type and whether the metric is corrected for signal leakage.

Connectivity metric Abbreviation
Coupling

| Model training
BRRR models were trained using MEG power spectra and FC metrics (see Table 1).One model was trained with power spectra, while with each FC metric, five different models were trained, one for each frequency band (see Figure 1a  for the locations of the included and excluded sensors).
Results for the models trained with sensor space data can be found in the Supplementary material.The steps presented in Leppäaho et al. (2019) were followed in training the models.The regression coefficients Ψ and the projection matrix Γ were initialized so that the identifier matrix X explains the maximum amount of variance in the target variables Y and that each latent component explains less variance than the previous component.The parameters were inferred using Gibbs sampling with 500 iterations, of which the first 250 iterations were discarded as the burn-in period.

| Model accuracy
Each subject had MEG data from three consecutive resting-state sessions.In the following, we will use the naming convention of the HCP dataset and refer to these resting-state sessions as "3-Restin", "4-Restin" and "5-Restin", respectively.The models were trained with the first resting-state session (3-Restin), while the two subsequent measurement sessions (4-Restin, 5-Restin) were used as test sets in computing the accuracy of the model.Results for models trained with the two other sessions can be found in the Supplementary material.
Accuracy was estimated as follows.First, the inverse of the projection matrix Γ (computed from the first resting-state session of the subjects) was used to project all the resting-state data to the latent space.In the latent space, cosine distances were calculated between the components of the training sessions and the components of the test sessions (Figure 1c).Cosine distance was chosen because it is rotationally invariant.If the training session closest to a test session was of the same experimental subject, classification was considered successful for that test session.Since each subject had two test sessions, two predictions were made for each subject, so the number of correct classifications per subject was either 0, 1 or 2. The total number of correct classifications was divided by the total number of predictions (2N), resulting in the accuracy score of each trained model.
The identification accuracies computed from the latent components were compared with correlation-based identification accuracies.Correlation-based individual identification from MEG power spectra and FC have been previously reported by da Silva Castanheira et al. ( 2021) and Sareen et al. (2021).In correlation-based identification, Pearson correlation coefficients are computed between the power spectra or FC matrices of separate measurement sessions across all the subjects.Predictions are based on the highest correlation values between the measurement sessions.The final identification accuracy is the number of correct identifications divided by the number of predictions.In the present work, the reported correlation-based identification accuracies are averages across all the resting-state session pairs (i.e., 3-Restin vs. 4-Restin, 3-Restin vs. 5-Restin, 4-Restin vs. 5-Restin).
In evaluating the BRRR modelling approach, we used the correlation-based fingerprinting accuracy as a "gold standard", with the aim to reach, in the best case, approximately equally good accuracy.Our aim was to determine whether BRRR can model the previously reported differences between the different metrics and frequency bands and what amount of model complexity (i.e., number of latent components) is required to capture the individual uniqueness of the metrics.

| Model comparison
Models were compared as follows.As was the case in computing the accuracy, the inverse of the projection matrix Γ was used to project all the resting-state sessions to the latent space.In the latent space, cosine distances were computed between all the measurement sessions (three per subject) across all the subjects (N ¼ 89), resulting in a distance matrix with the size 267 Â 267 (3N Â 3N).The distance matrix was considered to represent the structure of the latent space of a given model.The distance matrices of different models were compared by computing their correlation with the Mantel test (Mantel, 1967) (Figure 1d).In the present work, the Mantel test was used to assess the similarity of the latent space of the different models.We compared distance matrices, instead of directly comparing the latent components, since rotational variance can make individual components appear different in different models even when, as a whole, they span similar (but rotated) latent space.The use of distance matrices also enables comparison of models with different numbers of latent components.

| Cross-validation
Leave-one-out cross-validation was performed with models trained with the best-performing metrics to get an indication of their ability to classify unseen subjects.Since the results of Leppäaho et al. (2019) showed that spectral power relates siblings to each other, we assumed that if the sibling of a given test subject was included in the training data, the predictive performance of the model might be inflated for that test subject, as the model would have already partly seen the test subject by proxy of their sibling.Therefore, with each training and test data split, we also excluded the sibling of the test subject from the training data.
The inverse of the projection matrix Γ was used to project all the subjects and resting-state data to the latent space, but predictions were performed only for the unseen subject.In the latent space, cosine distances were computed between the components of the first restingstate session (3-Restin) of the test subject and the subsequent resting-state sessions (4-Restin, 5-Restin) of all the subjects.If the closest 4-Restin and/or 5-Restin was from the test subject, classification of that session of the test subject was considered successful.Thus, the number of correct classifications for each test subject was either 0, 1 or 2, and accordingly the accuracy of each crossvalidation model was either 0, 0.5 or 1.The reported cross-validation accuracies for each metric and frequency band are the average accuracies of these models.

| Correlation-based fingerprinting
Pearson correlation coefficients were computed between the power spectra or FC matrices of separate measurement sessions to estimate their similarity.The success of individual identification based on the correlations depended on the choice of frequency band and FC metric.With non-leakage-corrected metrics, both withinsubject (self) and between-subject (other) correlations were higher than with leakage-corrected metrics (Figure 2a).Despite the high values of between-subject correlations of the non-leakage-corrected metrics, the within-subject correlation distributions were quite separable from the between-subject correlations, allowing for high individual identification accuracy between the measurement sessions (Figure 2b).

| BRRR model-based fingerprinting
The BRRR models were trained with the first restingstate session (3-Restin), while the two subsequent measurement sessions (4-Restin, 5-Restin) were used as test sets in computing the accuracy of the model.Results for models trained with the two other sessions can be found in the Supplementary material.

| Effect of the number of latent components
The number of the BRRR latent components had a notable effect on the success of the model: a major increase in accuracy was seen between two and six components, and subsequent increments in the number of components had a diminishing benefit in increasing the accuracy (Figure 3a).Beyond six components, the relative distances between the subjects did not change much with additional components (see Figure 3b for PLV and Figures S2-6 for several other metrics).The correlations between distances in the latent spaces with six and 10 components were r ≥ 0:9 (p < 0:001) with most of the training data.In the following analysis, we use 20 latent components in all comparisons, as after that point improvements were minimal.

| Effect of the metric and frequency band
The different models trained with FC data varied substantially in their identification accuracy (Figure 4a, Table 2).Models trained with non-leakage-corrected phase metrics (Coh, PLV) produced the highest accuracies (Coh 96.1% and PLV 95.5% for data at gamma band).Other FC metrics did not reach > 90% accuracy, but gamma band of PLM and AEC came close (both 88.2%).The leakage-corrected FC metrics performed generally quite poorly, with the exception of PLM.Models trained with power spectra performed better than the models trained with leakage-corrected FC metrics (excluding PLM).With 20 latent components extracted from power spectra, it was possible to attain individual identification accuracy of 79.8%.Studying connectivity metrics computed over the whole range of frequencies that were used for the power spectra (1 À 90 Hz) we obtained a maximum classification accuracy of 73.6% for the leakagecorrected methods (PLM), and a maximum classification accuracy of 93.8% for the non-leakage corrected methods (PLV).
The FC metrics that are mathematically related produced distances in the latent spaces that were highly similar: correlations of r > 0:9 (p < 0:001) were observed, on the one hand, between the distance matrices of models trained with non-leakage-corrected phase metrics (Coh, PLV) and, on the other hand, between models trained with leakage-corrected phase metrics (wPLI and iCoh; PLI and iPLV) (Figure 5a).All the correlations between FC metrics and power spectra were below 0.4 (Figure S7).
The correlations among FC metrics were affected by the frequency band (Figure 5a).In general, alpha band produced the highest correlation values between the leakage-corrected phase metrics.PLM did not generally correlate strongly with other leakage-corrected phase metrics, but those correlations were highest in the alpha band (r ¼ 0:52 À 0:55, p < 0:001).Notably, in the gamma band, PLM diverged from the other leakage-corrected phase metrics by showing fairly high correlation with the non-leakage-corrected phase metrics (r ¼ 0:66, p < 0:001) and with AEC (r ¼ 0:54, p < 0:001).AEC also correlated with Coh and PLV (r ¼ 0:77, p < 0:001).The correlations between PLM, AEC and the non-leakage-corrected phase metrics (Coh, PLV) were much lower in the other frequency bands.AEC also correlated with both of its leakage-corrected versions, but the correlations were higher with pairwise orthogonalization (cAECp) (r > 0:9, p < 0:001; non-gamma bands).
The models trained with different frequency bands of non-leakage-corrected phase metrics (Coh, PLV) yielded very similar distances in the latent spaces (PLV shown in Figure 5b).With both Coh and PLV, the correlations among non-alpha frequency bands were r > 0:7 (p < 0:001).The highest correlations were observed between theta and delta bands (Coh: r ¼ 0:80; PLV: r ¼ 0:81) and between beta and gamma bands (Coh: r ¼ 0:83; PLV: r ¼ 0:84).When using leakage-corrected phase metrics, there was virtually no similarity between the models of different frequency bands (all correlations were below 0.1; see wPLI in Figure 5b).PLM was an exception to this finding, but even the PLM correlations remained relatively low, the highest value being between beta and gamma bands (r ¼ 0:30, p < 0:001).Amplitude metrics showed a somewhat different pattern: while correlations between different frequency bands of AEC were overall lower than those of PLV and Coh, applying leakage-correction to AEC increased some correlations but decreased others (AEC and cAECs shown in Figure 5b).

| Effect of the measurement session and source reconstruction
Models trained with 4-Restin and 5-Restin (Figures S8-9) produced similar accuracies as models trained with 3-Restin (presented above).The distance matrices of the models trained with 20 latent components and different measurement sessions were also similar, r ≥ 0:9 (p < 0:001), with most of the training data .Notably with the leakage-corrected phase metrics the similarity depended on the frequency bands (r > 0:9 only in alpha band; Figure S15).
Overall, similar differences among the FC metrics were observed with sensor-space and source-space data, but the accuracies of the sensor-space models were lower than those of the source-space models (Figure S16).However, power spectra yielded slightly higher accuracies with sensor-space than source-space data.The degree to which source reconstruction affected the distances in the latent space depended on the metric (Figure S17).Non-leakage-corrected phase metrics (Coh, PLV) had the least similar distances in the latent spaces between sensor-space and source-space models (r < 0:3 in all frequency bands).In contrast, the amplitude metrics and PLM yielded very similar distances in the latent spaces (r > 0:8, p < 0:001), but this depended on the frequency band.
computed for models trained with metrics and with different numbers of latent components K. (b) Mantel test correlations between distance matrices of models trained with PLV and different numbers of latent components K.

| Cross-validation
Leave-one-out cross-validation was performed on the best performing FC metrics (PLM, AEC, PLV) and the power spectra (Figure 4b) using 20 latent components.Coh was excluded since its results in the previous analyses were nearly identical to PLV.The cross-validation performance followed the trend of the models trained with the full dataset, but the accuracies were expectedly lower.The highest predictive accuracies of the FC models were again in the gamma band (PLM 67.4%, AEC 75.8%,PLV 86.5%).PLV was the best performing FC metric, producing > 80% accuracy in the non-theta frequency bands.The cross-validation model of power spectra performed  almost as well as the model trained with the full dataset, producing an accuracy of 77.5% in the leave-one-out cross-validation.

| DISCUSSION
Previously, BRRR was successfully used to build a low-dimensional neurofunctional model for identifying individuals from sensor-level resting-state MEG power spectra (Leppäaho et al., 2019).The present work aimed to determine whether the same general methodological framework could be applied to model individuality of FC between brain areas, as well as to evaluate whether the BRRR model generalizes to unseen individuals.A further key aim for the potential future use of the modelling of individual differences was to examine whether the observed distance patterns among individuals remained largely unchanged between the well-functioning models, or whether the different neural metrics would reveal distinct dis/similarity patterns among individuals.The performance of BRRR models mostly reflected the differences that are observed in correlation-based fingerprinting.Namely, the performance of a model depended strongly on the metric and frequency band that were used: Non-leakage-corrected FC metrics performed better than leakage-corrected FC metrics, and phase metrics performed better than amplitude metrics.The best performing BRRR models, with just 20 latent components, were able to identify individuals between measurement sessions with over 90% accuracy, approaching the highest correlation-based accuracies.The accuracy was reduced when classifying the measurement sessions of unseen subjects, but an accuracy exceeding 80% was still achieved with some of the models.
Based on our results, BRRR seems to behave well as a model of individual differences.It is simple, stable, and generalizable.Firstly, it can produce a low-dimensional representation of the various metrics, and the resulting low-dimensional latent components can be used to identify individuals between measurement sessions.Secondly, models trained with different measurement sessions produce very similar latent spaces, which remain quite consistent even when the model complexity is increased by adding more components.Thirdly, the model is predictive of novel subjects.These results provide a solid foundation for applying BRRR in future investigations of neurofunctional fingerprints.

| Individual neural fingerprints can be represented by just 20 components
As expected, the accuracy of all the models increased with the number of latent components, but the distances in the latent spaces were quite stable already with 6 components.These observations align with how the latentnoise BRRR is formulated, that is, the amount of additional variance explained decreases with each added component and the number of components needed is relatively low (Gillberg et al., 2016).We chose to use 20 components for our analyses, since the improvements gained with more than 20 components were minimal with most metrics and we wanted to keep the models as simple as possible.With certain FC metrics and just 20 components, it was possible to attain an identification accuracy that was comparable with the correlation-based accuracy.The FC matrices that we used in training the BRRR models (and in the correlation-based identification) contained 6670 unique values, so the reduction to just 20 values was notable.
Overall, sensor-space and source-space data produced similar results between the metrics, but the accuracies of the sensor-space models tended to be slightly lower.Since the difference between the two types of data were not very large, both approaches are viable for training BRRR models, each with its advantages and disadvantages.
Training models with sensor-space data has the obvious advantage that no anatomical data of the subjects is needed.In contrast, using source-space data makes it possible to apply the model across datasets, even when different MEG devices have been used to collect the data, as long as the same atlas is used in source reconstruction.
4.2 | Some FC metrics and frequency bands are more individually distinctive than others The performance of a BRRR model depended strongly on the FC metric and frequency band that was used in training the model.Mathematically related metrics produced similar results, both in terms of identification accuracy and the distances in the latent space.In contrast, differences were observed between amplitude-based and phase-based metrics as well as between leakage-corrected and non-leakage-corrected metrics.Namely, phase metrics performed better than amplitude metrics, and non-leakage-corrected FC metrics (AEC, Coh, PLV) performed better than their corresponding leakage-corrected versions (cAEC, iCoh, iPLV).The non-leakage-corrected phase metrics (Coh, PLV) had good performance at all frequency bands, but there was a mild drop in accuracy when moving towards the lower frequency bands.They also produced quite similar distances in the latent spaces in different frequency bands, especially between neighbouring frequency bands (delta and theta, beta and gamma).Both the non-leakage-corrected phase metrics and AEC had their best performance in gamma band.
Most of the leakage-corrected phase metrics (except PLM) had their best performance in alpha band, whereas PLM had its best performance in gamma and beta bands.There was almost no similarity between the distances in the latent spaces of non-leakage-corrected phase metrics (Coh, PLV) and their leakage-corrected counterparts (iPLV, iCoh).In addition, the distances in the latent spaces of the leakage-corrected phase metrics were quite unique in each frequency band, which aligns with the results of Marzetti et al. (2013), who reported consistent non-zero phase interactions that were frequency specific.In contrast, the leakage-corrected versions of AEC did yield similar distances in the latent spaces as nonleakage-corrected AEC, although the correlations were relatively low in the gamma band.Overall, it appears that zero-lag interactions are important for the distances in the latent space, but their importance varies between amplitude and phase metrics and across frequency bands.
PLM behaved differently from the other leakagecorrected phase metrics (iCoh, iPLV, PLI, wPLI).While the other leakage-corrected phase metrics yielded very similar results in terms of their identification accuracy, PLM had an identification accuracy that was comparable with the non-leakage-corrected metrics (AEC, Coh, PLV).The good performance of PLM was not entirely unexpected, since Sareen et al. (2021) have reported high correlation-based fingerprinting success rate for PLM.The authors of PLM (Baselice et al., 2019) highlight two factors that differentiate PLM from other leakagecorrected phase metrics: First, it can measure the dependency between signals even when the phase difference between them is not constant over time.Secondly, it is more resilient to noise.PLM's exceptionally good performance, compared with the other leakage-corrected phase metrics, is likely related to these two factors.

| Modelling results agree with correlation-based results but present also new phenomena
The differences between the BRRR models reflected the differences seen in the correlation-based identification: The same metrics that performed well in correlationbased identification also yielded BRRR models that performed well.Our correlation-based results were mostly in agreement with previous reports (Colclough et al., 2016;Sareen et al., 2021).The distributions of the correlation values showed that non-leakage-corrected FC metrics produced connectivity patterns that were highly similar both within a particular subject and with other subjects, when compared with leakage-corrected FC metrics.However, even though the between-subject correlations were high, the within-subject correlations were even higher, so there was enough subject-specific uniqueness in the connectivity values to allow for individuals to be differentiated from one another.
There were also some notable differences between the BRRR models and correlation-based fingerprinting.Firstly, the accuracy of the BRRR models trained with leakage-corrected metrics was markedly lower than that of the corresponding correlation-based accuracies.This suggests that the individually distinctive connectivity patterns quantified by leakage-corrected metrics are too complex to be condensed into just 20 latent values.Secondly, the BRRR modelling resulted in a strong difference between the higher and lower frequency bands of AEC and PLM, unlike the correlation-based results.This would indicate that the individually distinctive connectivity patterns quantified by these metrics can be captured by 20 values at high frequencies but are too complex at lower frequencies to be represented by just 20 values.
4.4 | Individual uniqueness may be linked to short-range zero-lag interactions Our results highlight non-leakage-corrected FC metrics and gamma band in their ability to differentiate between subjects.Non-leakage-corrected FC metrics are known to suffer from the effects of field spread even at the source level, leading to artificial zero-lag interactions (for review, see e.g., (Bastos & Schoffelen, 2016;O'Neill et al., 2015)).Non-leakage-corrected phase metrics appear to be dominated by local connections, which are absent in leakage-corrected phase metrics, while leakage-corrected phase metrics can detect long distance interactions (Sekihara et al., 2011;Stam et al., 2007).Long-distance connectivity is also more prominent in leakage-corrected amplitude coupling in comparison with non-leakage-corrected one (Colclough et al., 2015;Hipp et al., 2012).The local and long-range amplitude coupling patterns appear to be also frequency dependent: long-range interactions are more prominent in lower frequency bands in comparison with gamma band, which is dominated by local interactions (de Pasquale et al., 2010;Liu et al., 2010).Together, the previous work and our present results indicate that local connections are more informative in differentiating between subjects than longrange connections or, at least, that local connectivity patterns are more effectively represented by 20 latent values.To what degree the latent components of the nonleakage-corrected FC metrics reflect differences in true connectivity rather than differences in artificial connections due to field spread remains still unclear.However, the good accuracy obtained with PLM implies that learning latent components that are successful in individual identification is not necessarily dependent on zero-lag interactions.

| Limitations and outlook
The previous correlation-based fingerprinting results on this same HCP MEG data set could be convincingly reproduced with the present BRRR model-based approach.Furthermore, the results of the leave-one-out cross-validation suggest that, with certain FC metrics and power spectra, BRRR is applicable to unseen subjects.Nonetheless, to confirm the generalizability of the BRRR approach, future work will need to train and test the model across completely different datasets.
We used the distance matrices among the subjects as a means of quantifying the latent spaces of the different models.The distances in the latent space remained relatively similar with the increasing number of BRRR components, which speaks to the stability of the BRRR method.The intersubject distance patterns were also remarkably similar between mathematically related FC methods, as well as between gamma-range PLM, AEC, Coh and PLV that yielded the best identification accuracies.However, the intersubject distance patterns differed notably between the FC metrics and power spectra, both showing high identification accuracies.It remains to be explored whether and how these different intersubject relationships emerging from different neural metrics link to individual variation of neurofunctional task effects and behaviour.
Another intriguing question for future studies relates to the interpretation of the latent components.What are the neurofunctional factors that they represent?The differences between the models trained with leakagecorrected and non-leakage-corrected FC metrics suggest that short-range connectivity patterns may be well represented by a low dimensional model, but the relationship between the latent components and specific neural activation patterns requires further studies.
The main limitation of our study relates to the choice of parcellation scheme, which is likely to affect the BRRR modelling.Already in the correlation-based identification results there were minor differences in our results in comparison with those of Colclough et al. (2016) and Sareen et al. (2021), who employed different atlases.Relatedly, different source reconstruction approaches may also influence the FC data and thus the predictive performance of BRRR.

| CONCLUSION
In this study, we examined resting-state MEG power spectra and FC metrics in building a stable lowdimensional model of individual differences.Our results show that BRRR can be used to differentiate individuals with over 90% accuracy, with as few as 20 components, and even predict unseen subjects with over 80% accuracy.What remains to be further investigated in the future is the generalizability of the models: Can the models trained with the HCP dataset generalize to other datasets and differentiate individuals who are from different demographics and whose data have been collected at different MEG sites?Another major question is the interpretation of the latent components: how are they related to neurofunctional factors?Influence of other variables, such as the number of subjects included in training, use of task data and different classification schemes, also need to be inquired into to gain a more comprehensive understanding of the limits and capabilities of BRRR modelling.Our present results serve as a solid starting point for these future explorations.
-b); training a single model per FC metric would have been computationally overly demanding.The number of observations N (i.e., the experimental subjects) was 89 with all models, but the dimension P of the target variable matrix Y depended on the feature.It was 21 Â 116 (frequency bands Â ROIs) with power spectra and 116 Â 116 (ROIs Â ROIs) with FC.The identifier matrix X identified each subject as their own class.When the models were trained with sensor space data, only every second channel (A1, A3, A5, …) was used to reduce the computation time, so the dimensions were 124 Â 124 with FC and 21 Â 124 with power spectra (see FigureS1

F
I G U R E 1 (a)-(b) Overview of the analysis pipeline for using power spectra (a) and FC metrics (b) as features.BRRR models were trained with each feature to learn latent components that maximally differentiate between individuals.In the case of FC, separate models were trained with data at each frequency band.(c) The latent components of one model visualized in 2-dimensional space using t-distributed stochastic neighbor embedding (t-SNE) (van derMaaten & Hinton, 2008).Each colour denotes the measurement sessions of one subject.Cosine distances were computed between the latent components of the training sessions and test sessions across all the subjects.(d) Distance matrices from two different models.The distance matrix of a given BRRR model contains the cosine distances between the latent components of all the subjects N and three measurement sessions per subject.Mantel test was computed between two distance matrices to assess the similarity of their latent spaces.
The combined distributions of between-session correlation values of all the subjects.Red: correlation values between measurements of the same subject (self).Blue: correlation values between measurements of different subjects (other).(b) Accuracy of correlationbased identification for different metrics.

F
I G U R E 4 (a) Identification accuracies computed for models trained with different metrics and frequency bands.(b) Identification accuracies of leave-one-out cross-validation with the best metrics.Twenty latent components were used.T A B L E 2 Identification accuracies (%) for 20 latent components extracted from different FC metrics in different frequency bands.

F
I G U R E 5 (a) Mantel test correlations between distance matrices of models trained with different FC metrics.(b) Mantel test correlations between distance matrices of models trained with different frequency bands of PLV, PLM, wPLI, AEC and cAECs.Twenty latent components were used.