An F-ratio-based method for estimating the number of active sources in MEG

Introduction Magnetoencephalography (MEG) is a powerful technique for studying the human brain function. However, accurately estimating the number of sources that contribute to the MEG recordings remains a challenging problem due to the low signal-to-noise ratio (SNR), the presence of correlated sources, inaccuracies in head modeling, and variations in individual anatomy. Methods To address these issues, our study introduces a robust method for accurately estimating the number of active sources in the brain based on the F-ratio statistical approach, which allows for a comparison between a full model with a higher number of sources and a reduced model with fewer sources. Using this approach, we developed a formal statistical procedure that sequentially increases the number of sources in the multiple dipole localization problem until all sources are found. Results Our results revealed that the selection of thresholds plays a critical role in determining the method's overall performance, and appropriate thresholds needed to be adjusted for the number of sources and SNR levels, while they remained largely invariant to different inter-source correlations, translational modeling inaccuracies, and different cortical anatomies. By identifying optimal thresholds and validating our F-ratio-based method in simulated, real phantom, and human MEG data, we demonstrated the superiority of our F-ratio-based method over existing state-of-the-art statistical approaches, such as the Akaike Information Criterion (AIC) and Minimum Description Length (MDL). Discussion Overall, when tuned for optimal selection of thresholds, our method offers researchers a precise tool to estimate the true number of active brain sources and accurately model brain function.


. Introduction
Magnetoencephalography (MEG) is a powerful non-invasive neuroimaging technique that offers high temporal resolution for studying human brain function (Hämäläinen et al., 1993;Baillet, 2017).Localization of MEG sources has garnered significant interest in recent years since it can reveal the origins of neural signals and offer valuable insights into the complex workings of the human brain.By identifying the sources of neural activity, researchers can study the underlying mechanisms of cognition, perception, and other brain functions (Ahveninen et al., 2006;Giorgetta et al., 2013;Klepp et al., 2015;Pancholi et al., 2023).Additionally, source localization can aid in diagnosing and studying neurological disorders and identifying abnormal brain activity (Oishi et al., 2002;Westlake et al., 2012;Wilkinson et al., 2020;Xu et al., 2021;Giri, 2022;Giri et al., 2022).
MEG source localization methods typically involve solving an inverse problem, which entails estimating current sources within the brain based on the measured MEG data.This is challenging because the measured signals are influenced by several factors, such as the geometry and conductivity of the head, sensor noise, and the ill-posed nature of the problem.Mathematically, the localization problem can be cast as finding the location and moment of the set of dipoles whose field best matches the M/EEG measurements (Mosher et al., 1992).Localization methods can be broadly categorized into distributed and discrete solutions.
Distributed source imaging approaches aim to estimate a density map of active dipoles across the entire cortex.Commonly used methods include minimum norm estimators (MNE) (Hämäläinen et al., 1993(Hämäläinen et al., , 1994)), dynamic statistical parametric mapping (dSPM) (Dale et al., 2020), and standardized low-resolution electromagnetic tomography (sLORETA) (Pascual-Marqui et al., 2002).However, these methods assume a significantly larger number of unknown sources in a discrete surface or volumetric grid compared to the number of MEG sensors.The ill-posed nature of the problem poses a significant challenge, especially in the presence of multiple active regions in the brain (Darvas et al., 2004).Nonlinear source estimation methods, such as Mixed Norm Estimate (MxNE) (Strohmeier et al., 2016) and time-frequency mixednorm estimates (TF-MxNE) (Gramfort et al., 2013), address this issue by incorporating l1-norm penalty regularizers that favor sparse collections of focal dipolar sources.Other sparse approaches include hierarchical reconstructions of cortical and subcortical sources (Gramfort et al., 2012;Krishnaswamy et al., 2017;Rezaei et al., 2021).While these methods have shown some success, they tend to be computationally demanding and have limited accuracy when dealing with complex multi-dipole configurations.
On the other hand, discrete multiple dipole localization methods avoid the ill-posedness associated with distributed methods by finding a small set of equivalent current dipoles (ECDs) whose field best matches the M/EEG measurements in a least-squares sense (Mosher et al., 1992).Dipole localization methods offer a more classical approach to brain source localization and provide more intuitive interpretations of brain activity by estimating the location, orientation and amplitude of neural sources.The most well-known methods are beamformers (Van Veen et al., 1997;Vrba and Robinson, 2001) and MUltiple SIgnal Classification (MUSIC) (Mosher et al., 1992), and their recursive variants Recursively Applied and Projected MUSIC (RAP-MUSIC) (Mosher and Leahy, 1999), Truncated RAP-MUSIC (Mäkelä et al., 2018), andRAP Beamformer (Ilmoniemi andSarvas, 2019).While recursive variants generally perform better than their non-recursive counterparts, they still face limitations such as reduced effectiveness, reliance on high signal-to-noise ratio (SNR), and potential cancellation of correlated sources.Recent advancements in this field have addressed some of these concerns, including Alternating Projections (AP) (Adler et al., 2022), doublescanning (DS-MUSIC) (Mäkelä et al., 2017;Ilmoniemi and Sarvas, 2019), hemispherical harmonics MUSIC (HSH-MUSIC) (Giri et al., 2018), head harmonics MUSIC (H 2 MUSIC) (Giri et al., 2019), and Flex-MUSIC (Hecker et al., 2023).Estimating the number of independent signal components is a prerequisite for dipole localization methods to accurately estimate dipole sources.However, determining the correct number of active sources contributing to the recorded signals remains a fundamental challenge in MEG data analysis (Wendel et al., 2009), significantly impacting the success of brain source localization.We focus on addressing this specific problem here.
Estimating the number of active sources in MEG data poses significant challenges due to multiple factors.First, MEG signals generally exhibit a low SNR, which makes it difficult to differentiate between simultaneously active sources.Second, the presence of correlated sources adds complexity by potentially causing multiple sources to be mistakenly identified as a single source.Last, errors in translational head modeling and variations in individual anatomy introduce additional noise and variability, hampering accurate estimation of the location and strength of underlying sources.
Early attempts to estimate the number of dipoles relied on subjective thresholds (Bartlett, 1954;Lawley, 1956;Chen et al., 1991).These approaches involved setting a threshold that separated the eigenvalues of the data covariance matrix from the complete set of eigenvalues.Chen et al. (1991) proposed a method that detected the number of sources by imposing an upper bound on the eigenvalue magnitudes of the correlation matrix derived from the array output.In addition to conventional eigenvalue-based techniques, a few methods have also employed eigenvectors for estimating the number of sources (Di and Tian, 1984;Jiang and Ingram, 2004).
To overcome the limitations of subjective thresholds, two main classes of methods have been developed for estimating the number of signal sources using the distribution of the eigenvalues of the data covariance matrix.The first class involves techniques based on principal component analysis (PCA) (Green et al., 1988;Yao et al., 2018), independent component analysis (ICA) (Ikeda and Toyama, 2000), and factor analysis (Malinowski, 1977a,b).The second class consists of information theoretic approaches (Wax and Kailath, 1985;Knösche et al., 1998), such as the Akaike information criterion (AIC) (Akaike, 1974) and the minimum description length (MDL) (Rissanen, 1978;Schwarz, 1978).The work from Wax and Kailath (1985) derived the eigenvalue forms of AIC and MDL methods, which can be directly applied to array signal processing problems.These methods aim to strike a balance between model fit and complexity using principles from information theory.However, these approaches assume source independence, which is not always valid in the brain.As a result, they tend to perform poorly (Zhang et al., 1989;Chen et al., 1991;Salman et al., 2015;Yao et al., 2018), especially in the presence of correlated sources, noise, low SNR, and limited time samples.Hence, there is a need for more robust and accurate methods to estimate the number of active sources in MEG signals.
To address these limitations, we propose a robust method for accurately estimating the number of active sources in the brain using the F-ratio statistical approach.Our method introduces formal decision criteria that sequentially increase the number of sources in the multiple dipole localization problem until all sources are found.Our method is based on the F-ratio test, which is commonly used in statistics to compare the variances of two samples.It is sensitive to differences in the variances of the samples and can be used to determine whether adding a source to the model significantly improves the fit of the model to the data.The F-ratio statistical approach allows for a comparison between a full model with a higher number of sources and a reduced model with fewer sources.
We validated the F-ratio-based method on simulated, real phantom, and human MEG data, and compared its performance to that of other state-of-the-art statistical approaches, such as AIC and MDL.We found that the F-ratio-based method outperformed competing methods in terms of accuracy and reliability.One crucial aspect we investigated was the selection of appropriate thresholds for the F-ratio values.We found that this selection played a critical role in determining the overall performance of our method.Through systematic analyses, we identified optimal thresholds that needed to be adjusted according to the number of sources and SNR levels.Importantly, these thresholds exhibited remarkable consistency across different inter-source correlations, translational modeling inaccuracies, and cortical anatomies.When fine-tuned with the optimal selection of thresholds, our F-ratio-based method emerged as a precise and robust tool for estimating the true number of active sources in MEG data.

. Materials and methods
In this section, we provide a concise overview of the notations used to describe the measurement data, forward matrix, and sources.We also present the problem formulation for estimating multiple ECDs in the brain.Subsequently, we describe the F-ratio statistical procedure, which serves as the foundation to estimate the number of active sources in the brain, and outline the experimental procedures we use to assess performance.

. . Measurement model and notations
Consider an array of M MEG sensors detecting signals from Q ECD sources located at positions {p q } Q q=1 .At time t, the MEG measurement vector y(t) can be described as a superposition of the contributions from Q source signals {s q (t)} Q q=1 and additive noise: The topography l(p q ) of the qth dipole at location p q is defined as l(p q ) = L(p q )o, where L(p q ) ∈ R M×3 is the lead field matrix and o ∈ R 3×1 is the orientation vector.Depending on the problem, the orientation vector o may either be known, referred to as a fixedoriented dipole, or it may be unknown, referred to as freely-oriented.Additionally, the measurements are subject to the presence of additive white Gaussian noise, which is represented by n(t) ∈ R M×1 .
Several source localization methods exist for estimating the Q ECD sources, with each dipole source characterized by its location, orientation, and amplitude.Recently, we introduced a method called Alternating Projections (AP) (Adler et al., 2022), which offers several advantages.AP source localization method is robust to forward model errors, can handle high inter-source correlation values, and is effective even in low SNR scenarios.It is important to note that estimating the true number of active sources Q is a fundamental requirement for all the aforementioned dipole localization methods to accurately estimate the dipole source parameters.

. . F-ratio based method
The F-test is a widely used statistical technique that leverages the F-ratio to assess the presence of a significant difference between the variances of two data sets.In the context of determining the true number of sources, this technique holds particular value.It enables us to test the hypothesis that incorporating an additional source results in a substantial improvement in the variance accounted for by the model.By employing the F-test, we can make informed decisions regarding the optimal number of sources to include in order to provide the most accurate explanation for the observed data.
In probability theory and statistics, the F-statistic, also known as the F-ratio, is defined as the ratio between two independent chi-square distributions, denoted as X 1 ∼ χ 2 DOF 1 and X 2 ∼ χ 2 DOF 2 , where DOF 1 and DOF 2 represent the respective degrees of freedom.Mathematically, the F-ratio is expressed as: This formula provides a means to calculate the F-ratio by dividing the observed value of X 1 normalized by its degrees of freedom, DOF 1 , by the observed value of X 2 normalized by its degrees of freedom, DOF 2 .In this study, we use this idea to compare two hypothesized models based on the variance they explain.The first model, referred to as the "reduced model", explains the data with K number of sources: where y R (t) represent the estimated signals and n R (t) represents noise of a reduced model.On the other hand, the second model, called the "full model, " includes one more source compared to the reduced model, with K + 1 sources: where y F (t) represent the estimated signals and n F (t) represents noise of a full model.The estimation of y R (t) and y F (t) signals is achieved by solving an inverse problem using a dipole fitting method.By comparing the residual variance of these two models, we can assess their performance and determine the most appropriate model for the given data.Since we assume that the noise added to the measured signal is white Gaussian noise, it can be deduced that the sum of square errors between the measured signal y(t) and the estimated signals y R (t) and y F (t), denoted as 2 respectively, follows chi-square distributions.Therefore, the F-ratio test can be written as Supek and Aine (1993): where N is the number of time samples.For the fixed-oriented case, the degrees of freedom (DOF) for the reduced and full models are given by DOF R = MN − (3 + N)K and DOF F = MN − (3 + N)(K + 1), respectively.For the freely-oriented case, the DOF for the reduced and full models are DOF R = MN − (4 + N)K and DOF F = MN−(4+N)(K+1), respectively.These formulas account for a total of MN degrees of freedom, with three degrees of freedom deducted for position, one for orientation (in the freely-oriented case), and N for amplitude, for each dipole.
Note that the formulas for estimating the DOF assume independence among all data points, which is not the case in experimental data.To address this issue, in experimental data we implemented a two-step whitening process involving both temporal and spatial filtering.The first step involved temporal whitening, as depicted in the flowchart of the F-ratio method (Figure 1).The MEG data of each trial was whitened temporally using a six-order linear predictive coding (LPC) technique.This method helped alleviate temporal correlations within the data, reducing their impact on the results.Subsequently, the LPC-filtered trials were averaged.The second step was spatial whitening, which was achieved by applying a whitening filter derived from inverting the noise covariance matrix.This step further mitigated interdependencies among the data points, enhancing the reliability of the analysis.A more detailed discussion of the two-step whitening process is presented in Section 2.4.
It is important to note that the calculated F-ratio values are influenced by the residuals obtained after dipole fitting, and therefore, are dependent on the chosen source localization method used for solving the ECD localization problem in MEG.In our study, we specifically examined the behavior of the F-ratio statistical procedure when employing the AP localization method.The AP method solves the inverse problem iteratively and sequentially by minimizing the least-squares (LS) criterion.For a more detailed and comprehensive discussion of the AP method, we refer readers to Adler et al. (2022).
To estimate the number of sources, a systematic approach involving the comparison of a reduced model and a full model is illustrated in Figure 1.The formal comparison begins by initializing the reduced model with zero sources and the full model with one source to compute the F-ratio value.The decision regarding the number of sources involves comparing the resulting F-ratio value, denoted as F 0→1 , with a threshold value, F th .If the reduced model is rejected (F 0→1 > F th ), indicating evidence of at least one source, the analysis proceeds to compare a reduced model with one source against a full model with two sources, represented as F 1→2 .This sequential process continues, increasing the number of sources, until reaching a step where the reduced model cannot be rejected, providing an estimation of the true number of sources.

. . Performance evaluation with simulations
We evaluated the performance of the F-ratio method in diverse simulated scenarios, considering variations in the number of sources, SNR levels, inter-source correlations, translational modeling errors, and cortical anatomies.
The SNR was defined as the ratio between the Frobenius norm of the signal-magnetic-field spatiotemporal matrix and that of the noise matrix, following the approach described in Sekihara et al. (2001).To quantify inter-source correlation, we employed the Pearson's correlation coefficient.To establish the desired correlation among the sources, we utilized the Cholesky decomposition method.Initially, we generated fundamental cosine signals for each simulated source, with randomized phase and frequencies ranging from 10Hz to 30Hz.Next, we applied the Cholesky decomposition to factorize the symmetric positive definite target correlation matrix into the product of a lower triangular matrix and its conjugate transpose.By multiplying lower triangular matrix with the fundamental cosine signals, we generated a set of correlated dipole waveforms.This procedure ensured that the resulting source time courses closely matched the correlation coefficients specified by the target correlation matrix, thus incorporating the desired inter-source correlations in our simulations.
The sensor array geometry was based on the Megin Triux MEG system, which consists of a 306-channel probe unit with 204 planar gradiometer sensors and 102 magnetometer sensors.The MEG source space geometries were modeled using the cortical manifold extracted from MR data of adult humans, employing Freesurfer (Fischl et al., 2004).In our analysis, we used cortical anatomies from four different adult humans.Simulated sources were restricted to approximately 15,000 grid points distributed over the cortex.The reconstructed sources were estimated on a distinct grid of 50,000 points covering the cortex.To avoid the "inverse crime" problem, where identical parameters are used for data synthesis and inversion in an inverse problem, the simulation and reconstruction grids were non-overlapping with an average distance of 0.7 mm between neighboring points (Colton and Kress, 1998).The forward matrix for both grids was computed using the boundary element method implemented in OpenMEEG (Gramfort et al., 2010) within the BrainStorm software (Tadel et al., 2011).Simulated MEG data was generated by randomly selecting sources from the simulation grids.Gaussian white noise was then added to the MEG sensors to model instrumentation noise at specified SNR levels.In order to evaluate the effect of head model errors, we introduced translations to the reconstruction grid before computing the forward matrix.Last, we employed the AP method to solve the inverse problem in the fixed oriented case.All experiments were conducted with 100 Monte Carlo simulations to ensure statistical robustness.

. . Performance evaluation with a real phantom
We assessed the performance of the F-ratio method using the freely-oriented dipoles model with phantom data provided in the phantom tutorial (Taylor et al., 2016) of the Brainstorm software (Tadel et al., 2011).The phantom experiment was conducted using the Megin Neuromag system, which consists of a 306channel probe unit with 204 planar gradiometer sensors and 102 magnetometer sensors.
The data comprised MEG recordings obtained from the sequential activation of 32 artificial dipoles.To activate the phantom dipoles, an internal signal generator was used along with an external multiplexer box that connected the signal to each individual dipole.Each dipole was activated 20 times with an amplitude of 200 nAm, resulting in a total of 20 trials for each experimental condition.It is important to note that the chosen amplitude of 200 nAm falls within the range typically observed in inter-ictal spikes associated with epilepsy, as observed in raw data (Oishi et al., 2002).
For each dipole and each trial, the MEG data was whitened temporally using a six-order LPC technique.In particular, baseline data of a 200 ms pre-stimulus interval was used to compute the LPC coefficients of sixth order.These coefficients were then averaged across sensors and subsequently applied to the post-stimulus data modeled as an moving average (MA) filter.The purpose of this step was to eliminate temporal dependencies in the post-stimulus data, as observed in the auto-regressive (AR) model of baseline data.Following this, the LPC-filtered post-stimulus measurements were averaged across the 20 trials.In addition to temporal prewhitening, spatial prewhitening was also performed on the average data using a regularized noise covariance matrix in the Brainstorm software (Tadel et al., 2011).The regularization process included adding an identity matrix scaled to 10% of the largest eigenvalue of the noise covariance matrix.
To simulate the concurrent activation of multiple sources, we combined averaged data from different dipoles since only one dipole could be activated at a time.To introduce variability and avoid perfect coherence, we added a random delay ranging from 0 to 50 ms for each dipole.
The reconstruction source space was defined as a sphere centered within the MEG sensor array, with a radius of 64.5 mm.It was sampled using a regular volumetric grid of points with a resolution of 2.5 mm, resulting in a total of 56,762 grid points.The forward matrix was estimated based on a single sphere head model using the BrainStorm software (Tadel et al., 2011).The performance of the F-ratio method was evaluated using the AP method for localizing dipoles in the freely-oriented case.

. . Performance evaluation with human MEG data
The effectiveness of the F-ratio method in practical scenarios was assessed using human MEG data recorded from a single human participant during an auditory task.Prior to participation, the subject provided written informed consent, and the study was approved by the local ethics committee (Institutional Review Board of the Massachusetts Institute of Technology), following the principles of the Declaration of Helsinki.During the experiment, binaural sounds (beeps) were delivered to the subject using tubalinsert earphones.These auditory stimuli are known to elicit specific brain responses that are localized in the bilateral primary auditory cortex.A total of 166 trials were recorded, with an interstimulus interval of 2150 ms between each auditory stimulus.The MEG data were acquired using a MEGIN Triux MEG system, which includes a 306-channel probe unit consisting of 102 magnetometers and 204 planar gradiometers.
The forward matrix was estimated using BrainStorm based on an overlapping spheres head model.The reconstruction source space samples the brain volume in an adaptive way, with a higher density near the surface where we expect a higher spatial resolution due to the proximity to the sensors.The density decreased gradually toward the center of the brain, resulting in a total of 33,073 grid points.This grid is constructed using a specific adaptive algorithm in Brainstorm: It begins with a brain envelope containing 10,000 vertices as the initial number.Then, each previous layer is successively shrunk and downsampled by a factor of 3.This operation is repeated for 17 layers or until no more vertices are available.
The performance of the F-ratio method was evaluated using the AP method for localizing dipoles in the freely-oriented case.Before localization, the raw data underwent prewhitening in both the temporal and spatial domains following the same procedure as in the real phantom data.

. . Optimal modeling of MEG data requires the adjustment of F-ratio nominal thresholds
investigate the behavior of the F-ratio method under different experimental conditions, we conducted a thorough simulation analysis.Our objective was to assess the accuracy of estimating the true number of sources by varying the threshold values across various experimental scenarios.These scenarios included different numbers of active sources, varying SNR levels, inter-source correlation values, translational modeling errors, and cortical anatomies.
We observed that the accuracy of estimating the true number of sources using a specific F-ratio threshold was strongly influenced by the actual number of sources Q.This relationship is depicted in Figures 2A-C, where the estimation accuracy varied significantly for different values of Q. Notably, as the number of true sources increased, a lower F-ratio threshold was required to achieve higher performance.Similarly, we discovered a strong correlation between the accuracy of estimating the true number of sources and the SNR level.Figures 2D-F demonstrates this dependency for various SNR values.As the SNR level increased, a lower F-ratio threshold became necessary to achieve higher accuracy in estimating the true number of sources.
In contrast, we made the important observation that the optimal F-ratio threshold remained independent of the intersource correlation level.This finding is illustrated in Figures 2G-I, where we tested different inter-source correlation values (0.1, 0.5, and 0.9).The accuracy in estimating the number of sources peaked at the same threshold value for all correlation levels.This robustness indicates that the optimal F-ratio thresholds were not influenced by the inter-source correlation values of the active sources.Consequently, researchers may rely on a consistent threshold value regardless of the degree of correlation among the sources, enhancing the practical applicability and reliability of the F-ratio method.We obtained similarly robust results when testing the accuracy of estimating the true number of sources across different model errors (Figures 2J-L) and cortical anatomies (Figures 2M-O).In the case of model errors, we introduced registration errors by translating the reconstruction source space relative to the source simulation space.Specifically, we applied translations of 1 mm posterior (X-axis), right (Y-axis), and upward (Z-axis).Importantly, despite the presence of these registration errors, the F-ratio method remained highly robust.Similarly, when evaluating the F-ratio thresholds across the cortical anatomies of four different adult humans, we observed consistent and robust results.
In summary, our observations indicated that the performance of F-ratio thresholds varied significantly depending on the true number of sources and the SNR levels of the data.However, we found that F-ratio thresholds remained robust across different inter-source correlation values, translational model errors, and cortical anatomies.These findings highlight the need to adapt and optimize threshold procedures for the F-ratio test based on the specific number of sources and SNR levels in the data.In the next section, we determined these optimal thresholds.

. . Computation and evaluation of adjusted F-ratio thresholds in simulated data
In this section, we aimed to determine adjusted F-ratio thresholds for accurately estimating the number of active sources in MEG data.To accomplish this, we used a specific cortical anatomy as a reference (referred to as Anatomy 1).Optimal threshold values for Anatomy 1 were computed by identifying the F-ratio value that yielded the highest average accuracy in estimating the number of sources, considering the ρ ∈ {0.1, 0.5, 0.9} inter-source correlations and 1mm translational modeling errors in x, y, z.It may be noted that our study specifically examined the robustness of the method against translational modeling errors as a representative example.The resulting adjusted F-ratio thresholds were obtained for various SNR levels and numbers of sources, as depicted in Figure 3.Our findings indicate that higher threshold values are necessary in scenarios characterized by a high SNR and a low number of sources.
To assess the effectiveness of the adjusted threshold values obtained for the reference anatomy, we conducted tests on three additional cortical anatomies (Figure 4).The performance of the adjusted F-ratio thresholds in estimating the number of sources was evaluated at 0 dB SNR and correlation levels ρ ∈ {0.1, 0.5, 0.9}.Remarkably, the results showed that the performance of the adjusted thresholds was comparable to that of the reference cortical anatomy (Anatomy 1).These findings demonstrate the reliability and robustness of the calculated optimal F-ratio thresholds across a wide range of simulation scenarios, including variations in the number of sources, SNR levels, inter-source correlation values, translational modeling errors, and cortical anatomies.We proceeded by conducting a comparative analysis between the proposed F-ratio method with adjusted thresholds and two commonly used methods, namely the information criterion AIC and MDL method, for estimating the number of sources.These methods rely on likelihood functions derived from information theory to assess and choose the optimal model.The comparison results of the F-ratio, AIC, and MDL methods across different SNR conditions and correlation levels are depicted in Figure 5.The results correspond to cortical anatomy 4 with 1 mm translational modeling error in X. Remarkably, the proposed F-ratio method with adjusted thresholds outperformed both the AIC and MDL methods in terms of accuracy and reliability.

. . Performance of the F-ratio method in estimating the number of active dipoles in phantom data
We assessed the performance of the F-ratio method in estimating the number of active dipoles in phantom data (Figure 6A).The locations of the 32 artificial dipoles of the MEGIN phantom are shown in Figure 6B.To simulate the activation of multiple MEG sources simultaneously, we combined the data obtained from individually activated dipoles.To avoid perfect coherence, a random delay ranging from 0 to 50 ms was introduced between the dipole time courses.Figure 6C illustrates an example of MEG sensor data from two active dipoles with a randomly selected temporal delay of 29 ms.The time courses are displayed following temporal and spatial prewhitening, and had an estimated SNR of 5.5 dB.
We conducted 100 Monte Carlo repetitions of phantom data simulations for each scenario involving 0 to 5 active sources.For each repetition, we applied the optimal F-ratio thresholds (presented in Figure 3) based on the estimated SNR of the phantom data and the corresponding number of tested sources.The rationale behind using the same optimal F-ratio value from simulations to both the phantom and experimental data lies in the observation that, while the performance of F-ratio thresholds varied significantly based on the number of sources and SNR levels, these thresholds remained largely invariant to variations in inter-source correlations, translational modeling inaccuracies, and different cortical anatomies.
The performance of the F-ratio method in accurately estimating the true number of active dipoles is shown in Figure 6D.The method successfully identified the correct number of sources up to 2, surpassing both the AIC and MDL methods.In contrast, the latter methods failed entirely, with no correct estimations among the 100 simulated scenarios (results not depicted).The AIC method yielded mean estimates of 41 for no true sources and 42 for 1 to 5 true number of sources, respectively, across 100 Monte Carlo repetitions.Similarly, the MDL method estimated the number of sources as 41 in all cases of 0 to 5 true numbers of sources across the same 100 Monte Carlo repetitions.It was found that AIC and MDL consistently overestimated the number of sources.
It is important to note that although the performance of the F-ratio method in experimental data was not as remarkable as in the simulated data, we attribute this to two factors.First, the specific configuration of the phantom dipoles played a critical role.The 32 phantom dipoles were closely spaced and shared similar orientations.In the 100 Monte Carlo repetitions, sources were randomly chosen with no constraints in their simultaneous activation.Consequently, there was a substantial probability of selecting adjacent sources, and this likelihood increased as the number of sources grew (3, 4, and 5).Additionally, the time courses of the phantom dipoles had random delays ranging from 0 to 50 ms, resulting in instances of minimal delay and strong correlation.The combined effect of proximate source selection and small time course delay exacerbated the challenge of accurately estimating the number of sources, particularly when dealing with a larger number of sources.

. . Performance of the F-ratio method in estimating the number of sources in human auditory data
To evaluate the effectiveness of the F-ratio method in human data, we utilized it to analyze brain responses captured during an auditory task.We employed the AP source localization method to fit dipoles within the time interval of 100-130 ms relative to the onset of the auditory stimuli, which approximately corresponded to the peak MEG response.The SNR within this interval was determined to be 5.78 dB.
Figure 7A illustrates the adjusted F-ratio thresholds at 5.78 dB SNR, along with the corresponding values estimated from the human data.By employing the F-ratio sequential procedure, we observed the rejection of the reduced models with zero or one source, while finding no evidence to reject (and thus accepting) a model with two distinct sources.To further validate these findings, we plotted the sum of squares of residuals for the competing models, as shown in Figure 7B.This plot reinforces our conclusion of the presence of two sources, as there was no significant reduction in the sum of squares of residuals beyond the model order of 2. Additionally, we visualized the dipoles detected using the AP localization method for the cases of 1, 2, and 3 sources (Figures 7C-E).In the case of the two-source model, the dipoles were localized bilaterally and coincided with the well-established regions in the primary auditory cortices.
It is worth mentioning that alternative methods such as the AIC and MDL yielded different estimates for the number of sources.Specifically, the AIC method suggested 33 sources, while the MDL method indicated only 31 sources.These estimates far exceed the FIGURE Performance of adjusted F-ratio thresholds in simulated data.The thresholds were optimized for Anatomy and applied to three di erent anatomies (Anatomies , , and ).The F-ratio method was employed to estimate the number of sources under varying levels of source correlation (A) ρ = ., (B) ρ = ., and (C) ρ = ., while maintaining a signal-to-noise ratio of dB.
expected number of active sources in human auditory responses and are not in line with the existing knowledge in the field.
Lastly, we conducted a further analysis on the human data with a reduced SNR of approximately 0 dB, achieved by averaging a reduced number of 44 trials (originally 166).As depicted in Supplementary Figure 1, this analysis highlights the robustness of the F-ratio test in accurately detecting 2 sources, consistent with the findings from the higher SNR case.It was found that AIC and MDL consistently overestimated the number of sources even in the presence of lower SNR.This suggests that the overestimation nature of AIC and MDL (in phantom and real data) might be related to the complex nature of human experimental data.

. Discussion
We have proposed and validated an F-ratio-based method to detect the number of active sources in MEG data.We initially aimed to tune the method with a universal F-ratio threshold value that, once selected, could be applied across diverse simulation scenarios.However, we found that the performance of the Fratio method was inherently dependent on the SNR and number of sources.Thus, we have concluded that it is not feasible to determine a single threshold value that guarantees optimal performance across all cases.Instead, we proposed a methodology that adjusts F-ratio threshold values based on the estimated SNR and the corresponding number of tested sources.Our results demonstrated the reliability and robustness of the calculated optimal F-ratio thresholds across a wide range of simulation scenarios, including variations in the number of sources, SNR levels, inter-source correlation values, modeling errors, and cortical anatomies.
However, it is crucial to acknowledge that the adjusted threshold values obtained in this study are specific to the MEG system analyzed and may need to be adjusted for other devices or modalities, such as EEG.When applying the F-ratio method in different devices, it would be necessary to determine appropriate threshold values that are specific to each case.Similarly, the effectiveness of the proposed method is inherently linked to the choice of the source localization technique.In this study, we employed the AP method to compute and validate the F-ratio thresholds.However, different source localization methods may yield varying results and require different threshold adjustments.Therefore, it is crucial to determine the optimal threshold values for a particular set of experimental settings and source localization method to ensure accurate estimation of the true number of active brain sources.Despite these considerations, our proposed method provides researchers with a precise tool to estimate the true number of active brain sources and effectively model brain function.By calculating threshold values that are tailored to the specific modality and source localization method, researchers can enhance the accuracy and reliability of their source estimation process.
A study by Supek and Aine (1993) aimed to evaluate the efficacy of three statistical measures, namely percent of variance, reduced chi-square, and F-ratio, in determining the correct number of sources (model order).The authors advocated for the reduced chisquare method as a reliable measure of goodness-of-fit, whereas they were less favorable to the percent of variance and F-ratio because they ignored noise contributions.Although we agree that the percent of variance has limited utility, we contend that the efficacy of the F-ratio method was underestimated.Indeed, the simulation results presented in Supek and Aine (1993) demonstrated that the F-ratio remained stable across different noise levels and successfully identified the true number of sources.However, it is important to note that the study was confined to simulations on a simple spherical head model, lacked assessments using real data, and did not provide clear threshold decision criteria for determining the correct number of sources.In this work, we proposed a methodology to estimate the number of sources by employing F-ratio threshold values based on the estimated SNR and the corresponding number of tested sources.
The Bayesian multi-dipole estimation method, Sequential Semi-Analytic Monte-Carlo Estimation (SESAME) (Sommariva and Sorrentino, 2014;Sorrentino et al., 2014;Luria et al., 2020) is an iterative Monte Carlo algorithm that approximates the posterior distribution for an a-priori unknown number of dipoles.The output of SESAME is a posterior distribution for a variable number of dipoles and their parameters.From this distribution, a cortical probability map is computed, quantifying for each voxel the posterior probability of containing a dipolar source.Additionally, the method provides a point estimate of the dipole location, determined as the peak of the cortical probability map.In future research, it would be valuable to compare the performance of SESAME against our proposed method.
It is essential to acknowledge the fact that the performance of the F-ratio method in phantom data was not as remarkable as in the simulated data, we attribute this to two factors.First, the specific configuration of the phantom dipoles played a critical role.The 32 phantom dipoles were closely spaced and shared similar orientations.Additionally, the time courses of the phantom dipoles had random delays ranging from 0 to 50 ms, resulting in instances of minimal delay and strong correlation.The combined effect of proximate source selection and small time course delay exacerbated the challenge of accurately estimating the number of sources, particularly when dealing with a larger number of sources.
In the literature, it is well known that the localization errors tend to increase when dealing with sources that have small spatial separation.This phenomenon has also observed in the AP method, as we reported in our previous work (Adler et al., 2019).Second, estimating the number of sources beyond two is generally an exceptionally challenging problem in MEG.Consequently, a significant body of work has focused on solving the problem of localizing up to two or three sources within simulation or controlled environments (Mosher et al., 1992;Mosher and Leahy, 1999;Mäkelä et al., 2018;Adler et al., 2019;Ilmoniemi and Sarvas, 2019;Giri et al., 2020).
Finally, it is important to mention that the MEG human experiment revealed the presence of two active sources.Notably, the F-ratio method consistently exhibited excellent performance in simulations and phantom experiments with this specific number of sources.However, it is crucial to acknowledge that as the number of sources increases, the suitability of the F-Ratio method might diminish, necessitating further investigation in future studies.

. Conclusion
We have validated our F-ratio-based method on simulated, real phantom, and human MEG data.In comparison to other state-of-the-art statistical approaches like AIC and MDL, which rely on certain assumptions that often do not hold in real-world situations, our method demonstrated superior performance in terms of accuracy and reliability.One crucial aspect we emphasized is the selection of appropriate thresholds for the F-ratio values, which significantly impacts the overall performance of the method.We identified optimal thresholds and showed that these thresholds needed to be adjusted for the number of sources and SNR levels.Notably, these thresholds exhibited remarkable consistency across different inter-source correlations, head translation modeling errors, and cortical anatomies.Overall, by fine-tuning the selection of thresholds, our F-ratio-based method provides researchers with a precise and robust tool for accurately estimating the true number of active sources in MEG data.Further research is needed to explore and validate the proposed method in different modalities and with various source localization techniques.By refining the threshold determination process and investigating its applicability across different experimental conditions, we can extend the utility of this method to a wider range of neuroimaging studies and enhance our understanding of the underlying mechanisms of brain function.

FIGURE
FIGUREFlowchart of the F-ratio method for estimating the number of sources.

FIGURE
FIGUREAccuracy (%) of the F-ratio method for estimating the number of active MEG sources (on vertical axis) under di erent threshold values (on horizontal axis).Performance evaluation was conducted across various experimental conditions: (A-C) varying number of true sources Q, (D-F) SNR levels from -to dB, (G-I) inter-source correlation values ρ from . to ., (J-L) di erent model errors, and (M-O) di erent cortical anatomies.Each experimental condition was tested using Monte-Carlo repetitions to ensure statistical robustness.To account for head registration errors, we incorporated inaccuracies into the lead field matrix by applying a translation of mm posterior (X-axis), rightward (Y-axis), and upward (Z-axis).

FIGURE
FIGUREOptimal F-ratio threshold values, adjusted for the signal-to-noise ratio (SNR) level and the number of active sources in the MEG data.