Randomized Multiresolution Scanning in Focal and Fast E/MEG Sensing of Brain Activity with a Variable Depth

We focus on electro-/magnetoencephalography imaging of the neural activity and, in particular, finding a robust estimate for the primary current distribution via the hierarchical Bayesian model (HBM). Our aim is to develop a reasonably fast maximum a posteriori (MAP) estimation technique which would be applicable for both superficial and deep areas without specific a priori knowledge of the number or location of the activity. To enable source distinguishability for any depth, we introduce a randomized multiresolution scanning (RAMUS) approach in which the MAP estimate of the brain activity is varied during the reconstruction process. RAMUS aims to provide a robust and accurate imaging outcome for the whole brain, while maintaining the computational cost on an appropriate level. The inverse gamma (IG) distribution is applied as the primary hyperprior in order to achieve an optimal performance for the deep part of the brain. In this proof-of-the-concept study, we consider the detection of simultaneous thalamic and somatosensory activity via numerically simulated data modeling the 14-20 ms post-stimulus somatosensory evoked potential and field response to electrical wrist stimulation. Both a spherical and realistic model are utilized to analyze the source reconstruction discrepancies. In the numerically examined case, RAMUS was observed to enhance the visibility of deep components and also marginalizing the random effects of the discretization and optimization without a remarkable computation cost. A robust and accurate MAP estimate for the primary current density was obtained in both superficial and deep parts of the brain.


Introduction
This study concentrates on electro-/magnetoencephalography (E/MEG) imaging of the brain activity (He et al. 2018). The present focus is on the hierarchical Bayesian model (HBM) (Calvetti et al. 2009;Lucka et al. 2012) which allows one to find a focal and robust reconstruction by exploring a posterior probability distribution following from a conditionally Gaussian prior model. Our aim is, in particular, to develop a fast maximum a posteriori (MAP) estimation technique which would be applicable for both superficial and deep areas without additional a priori knowledge of the brain activity, such as physiological depth weighting (Calvetti et al. 2015(Calvetti et al. , 2018Homa et al. 2013). While high-density measurements (Seeber et al. 2019) and advanced signal processing strategies (Pizzo et al. 2019) have recently been shown to be essential in distinguishing deep activity, this study focuses on the importance to reduce the random effects of the numerical discretization and optimization errors on the reconstruction process.
We introduce a randomized multiresolution scanning (RAMUS) method in which the MAP estimate of the brain activity is refined gradually in the reconstruction procedure. RAMUS aims at reducing the random effects of the numerical discretization on the final estimate. It processes the well and ill-conditioned parts of the source space separately which has been suggested for ill-posed problems, e.g., in (Pursiainen 2008;Liu et al. 1995;Piana and Bertero 1997a). A multiresolution decomposition provides an approximative split between detectable and undetectable parts for different source depths, as the maximal source localization accuracy varies strongly with respect to the depth (Tarkiainen et al. 2003;Cuffin et al. 2001a, b;Grover 2016;Wang et al. 2009) with only the low resolution fluctuations being visible in the deep part of the brain (Pascual-Marqui et al. 1999;Pascual-Marqui 1999). At each resolution level, a MAP estimate is evaluated via the iterative alternating sequential (IAS) algorithm and the inverse gamma (IG) hyperprior which has been found to be advantageous for detecting deep activity (Calvetti et al. 2009).
The previous results suggest that HBM can find a focal solution deep in the head via the Markov chain Monte Carlo (MCMC) sampling techniques, especially, if the activity can be constrained into a region of interest (ROI) (Calvetti et al. 2009;Lucka et al. 2012). However, processing large data sets involving temporal measurement sequences with an advanced MCMC approach without a priori knowledge of a ROI might be computationally too expensive for the practical use. Therefore, finding a robust and fast approach to distinguish activity reliably is crucial regarding the practical applications. In this proof-of-the-concept study, we consider the detection of simultaneous somatosensory and thalamic activity with numerically simulated data. This setup models the detection of the somatosensory evoked potentials and fields (SEP/F) in response to the electrical stimulation of the median nerve, particularly, thalamic (deep) P14/ N14 and somatosensory (superficial) P20/N20 component peaked at 14 and 20 ms post-stimulus, respectively (Buchner et al. 1988(Buchner et al. , 1995(Buchner et al. , 1994aHaueisen et al. 2007;Attal and Schwartz 2013;Fuchs et al. 1998).
In the numerical experiments, both a spherical and realistic model has been used to analyze the source reconstruction discrepancies with RAMUS. The results suggest that a randomized set of decompositions (Mallat 1989;Clark et al. 1995) is essential to marginalize out the possible modeling errors due to projecting the source space into different resolution levels which, again, is necessary in order to achieve the depth-invariance of the final MAP estimate.

Observation Model
For the EEG source modelling, we employ the finite element method and the current preserving H(div) approach (Pursiainen 2012a;Pursiainen et al. 2016;Miinalainen et al. 2019) in which the primary current distribution of the neural activity is assumed to have a square-integrable divergence P ∈ H(div) = { |∇ ⋅ ∈ L 2 ( )} in the source space denoted by S . The observation model is where ∈ ℝ m is the measurement vector, ∈ ℝ m×3K is the lead field matrix, ∈ ℝ 3K is the unknown primary current distribution with K denoting the total number of the source positions, and ∈ ℝ m is the measurement noise vector which is modelled as Gaussian random variable with zero mean and covariance matrix of the form 2 ∈ ℝ m×m . In this numerical study, the diagonal covariance is used for simplicity as it allows fixing the noise level with a single parameter, i.e., the standard deviation . We refer to ℝ 3K as the source space S for the inverse problem of finding given the data . The number of sources is three times the number of their positions, as each position is assumed to have three sources oriented along the Cartesian coordinate axes.

Hierarchical Bayesian Model
In the HBM framework, the prior of is not fixed but random. It is determined by the realization of the so-called hyperparameter . The hyperparameter follows an a priori assumed distribution, i.e., the hyperprior. Consequently, the prior is a joint density given by p( , ) ∝ p( ) p( | ) of and . The conditional part of the prior p( | ) corresponds to a zero mean Gaussian density with a diagonal covariance matrix predicted by the hyperprior p( ) . The hyperparameter is of the same dimension as with each entry defining the variance of its respective entry in . The density of the hyperprior is long-tailed, implying that is likely to be a sparse vector with only few nonzeros, which is advantageous for finding a focal reconstruction of the brain activity. As a hyperprior, one can use, e.g., the gamma (G) or inverse gamma IG( | , 0 ) density (Calvetti et al. 2009), whose shape and scale are controlled by the parameters and 0 , respectively. IG is a conjugate prior for a Gaussian distribution with an unknown variance (here the conditional prior), meaning that the corresponding posterior (here the actual prior) is also Gaussian. Again, G is a conjugate prior with respect to the reciprocal of the variance (O' Hagan and Forster 2004).
The posterior probability density of , following from the classical Bayes formula (O'Hagan and Forster 2004), is of the form i.e., it is proportional to the product between the prior density p( , ) , and the likelihood function p( | ) ∝ exp(−(2 2 ) −1 ‖ − ‖ 2 ) given by the measurement noise model (Schmidt et al. 1999).
We consider finding the inverse estimate via the iterative alternating sequential (IAS) MAP estimation method  (Nagarajan et al. 2006), respectively, while the first step of the iteration concides with the classical minimum norm estimate (MNE) (Hämäläinen et al. 1993). A recent comparison between IAS and other brain activity reconstruction techniques can be found in (Calvetti et al. 2018). The numerical exploration of the posterior density p( , | ) is subject to the numerical discretization, i.e., the numerical definition of the source space S for and the resulting lead field matrix. We aim to reduce the effect of the discretization via the following two strategies motivating the introduction of the RAMUS approach: (1) The reduction of the source space is essential to improve the ability of a solver to recover focal sources both in deep and superficial locations. Furthermore, since a sparse source space results here in source reconstruction of low spatial resolution, a source space refinement during the reconstruction process of this study is crucial. (2) A randomized set of decompositions enables averaging out (marginalizing) the effect of the discretization error.
The theoretical justification of (1) and (2) are given in the following sections "Coarse-to-Fine Optimization" and "Randomized Scanning", respectively.

Coarse-to-Fine Optimization
The EEG source imaging problem is severely ill-posed (Grech et al. 2008) and it is well-known that most of the solvers suffer from depth bias effects (Pascual-Marqui 1999; Koulouri et al. 2017;Awan et al. 2018). A way to reduce the ill-conditioning in the computations is by introducing coarser (sparse) source space, i.e., regularization by discretization (Hansen 2010; Kirsch 2011), or by approximating the source distribution as a linear combination of spatial basis functions (redundant dictionaries) as proposed in Haufe et al. (2008). With a dimensionality reduction, the linear system to be solved is often over-determined and stable estimates can be obtained. However, this comes at a cost of poor resolution reconstructions due to large discretization errors. The idea of employing a multiresolution approach (Mallat 1989), where a progressive refinement in the source space is performed in order to obtain more accurate estimates, has been proposed for the E/MEG problem for example in Gavit et al. (2001); Malioutov et al. (2005). The source space S can be decomposed via the direct sum where is determined by the noise level. S + and S − represent the sets of the detectable and undetectable source distributions, respectively. If possible, it is advantageous to decompose S into S + and S − as, thereby, one can avoid source localization errors related to the indetectable distributions S − (Pursiainen 2008;Piana and Bertero 1997b;Liu et al. 1995). In E/MEG, a coarse enough source configuration can be distinguished, i.e., it belongs to S + , while a dense one has modes that cancel each other and might be indetectable, i.e., in S − (Fig. 1). The coarsity is specifically Fig. 1 1st from the left In E/MEG, a coarse enough source configuration can be distinguished, i.e., it belongs to S + , while a dense one has modes that cancel each other and might be indetectable, i.e., in S − (Fig. 1). 2nd and 3rd from the left: An example of subdividing the grey matter compartment to subdomains in the case of a coarse (center) and fine (right) resolution. Here, the sparsity factor s, i.e., the ratio between number of subdomains for two consequtive resolution levels, would be four. An example of mapping a subdomain from a coarse to fine resolution is given by {2} → {2, 28, 29, 30} important considering deep activity for which the magnitude of the lead field is comparably low and, therefore, any deep source configuration is likely to belong to S − . For a given lead field matrix , the maximum possible number of detectable sources and, thereby, the maximal dimension of S + is determined by the maximum number of nonzero singular values which coincides with the smaller dimension of , that is, the number of the data entries m.
In the coarse-to-fine reconstruction strategy, the aim is to first limit the source space S to a subspace S + by restricting its resolution, to gradually increase its resolution, and to eventually obtain an approximation for the whole space S . A nested set of restricted subspaces with different resolutions referred here to as a multiresolution decomposition is obtained recursively by selecting a uniformly random set of source positions from a given source space and associating those with the original set of positions through nearest interpolation. The coarsest resolution level is associated with the index = 1 . When moving from the -th resolution level to the ( + 1)-th one, the number of source positions is assumed to grow by a constant sparsity factor s > 1 . An example for a dual resolution decomposition and a mapping of the subdomains between them can be found in Fig. 1.
In the IAS MAP estimation process, once the activity has been found at a coarse reconstruction level, the support of the candidate solution will shrink along with the increasing resolution ( Fig. 2). That is, the size of the details found is subject to the resolution level. Therefore, the final estimate is found as a combination of the estimates obtained for the different levels. In order to distinguish the weakly detectable activity, especially, the deep components, the number of the dimensions in the initial set should be of the same size with m, following from the maximal dimensionality of S + .

Randomized Scanning
Since a sparse source space is likely to induce a bias to the consequent estimates, we propose to use a random set of (initial) sparse source spaces that aims to reduce the propagation of random discretization and optimization errors. The relationship between the global posterior optimizer * and k for the original source space S and its restriction S k , respectively, can be modeled via the equation where k and k represent a discretization and optimization error, respectively. Of these, k depends of the quality of the MAP optimization method and vanishes in the ideal case, while k is fixed. If the degrees of freedom in S 1 , S 2 , … , S D have an independent and identical random distribution, the respective discretization errors 1 , 2 , … , D can be modeled as independently and identically distributed random variables and, by the law of large numbers and the central limit theorem, the discretization error term 1 D ∑ D k=1 k of the mean is an asymptotically Gaussian variable with expectation ̃ and the rate of convergence 1 with respect to the number of source spaces (Liu 2001). Consequently, the random effects of the discretization errors can be marginalized via estimating * in multiple randomly (independently and identically) generated source spaces. The expectation ̃ can be regarded as the remaining systematic discretization error which is specific to the set S 1 , S 2 , … , S D , i.e., the resolution level, and is related, for example, to the relationship between the maximal achievable level of detail and the structure of the actual unknown brain activity.
Since the outcome of the optimization process for each given source space is a priori sensitive to the discretization errors, the estimate for k is found using the one for k−1 as the initial guess. This approach is motivated by the present gradually progressing coarse-to-fine subdivision due to which the subsequent optimizers will be nearly similar. We consider it necessary in order to maintain each estimate in the vicinity of the global optimum and, thereby, the norm of the optimization error k as small as possible. Namely, using a fixed initial guess might mean that, instead the global optimizer, a local one is found for some of the source spaces as depicted in Fig. 3. The global optimum might correspond to a situation in which both a superficial and deep source are detected, while the deep activity might be undetected at a local one.
Technically, updating the initial guess makes the estimate for k dependent on the previous one obtained for k−1 , i.e., the sequence of the estimates is a time-homogeneous Markov chain. We regard the present approach as a surrogate transition rule (Liu 2001) estimating the outcome of an ideal optimization method which would find the global optimum precisely with k = 0 , thereby, resulting in the identity Fig. 2 Once an approximation for a non-zero source has been found at a coarse resolution level (left) the its support will shrink at the finer levels (right) which will hold approximately, if the surrogate rule is accurate enough.

RAMUS
We propose the following algorithm for RAMUS to reduce the random discretization and optimization effects when finding a reconstruction for the unknown parameter with the IAS MAP estimation method.
1. Choose the desired number of the resolution levels L and the sparsity factor (the ratio of source counts) s between each level. The number of the sources at a given resolution level will be K = Ks ( −L) , where = 1, 2, … , L is the index of the resolution level, the larger the value of the index the finer the resolution. 2. For each resolution level = 1, 2, … , L , create a random uniformly distributed set of center points 1 , 2 , … , K . Find source point subsets B 1 , B 2 , … , B K applying the nearest interpolation scheme with respect to the center points. That is, each subset B j consists of those source positions of the total source space S , whose nearest neighbor with respect to 1 , 2 , … , K is j . The average number of source positions associated with B j is approximately given by the sparsity factor s. The resolution of this subdivision grows along the number of the center points. The unknown parameter is assumed to be constant in each subset, and the actual source count is assumed to stay unchanged regardless of the resolution. 3. Repeat the first two steps to generate a desired number D of independent multiresolution decompositions 1 , 4. Start the reconstruction process with the decomposition 1 and a suitably chosen initial guess (0) . 5. For decomposition k , find a reconstruction ( ) with the IAS MAP technique with the initial guess ( −1) for the resolution levels = 1, 2, … , L. 6. After going through all the decompositions, obtain the final estimate for the decomposition (basis) k as the normalized mean where the denominator follows from the need to balance out the effect of the multiplied source count following from the interpolation of a coarse level estimate to a denser resolution level. 7. If k < D , move to the next decomposition, i.e., update k → k + 1 , and repeat the previous step with the initial guess (k−1) for the resolution level = 1.

Obtain the final reconstruction as the mean:
Technically, this process is equivalent to first evaluating the mean (7) for each resolution level and then the normalized mean (6) over the different resolutions, showing that an approximation of the form (3) is, in fact, obtained for each set of independent and identically generated source spaces.
Since the final reconstruction is obtained as a mean over all the reconstruction levels, also the potential systematic discretization errors will be averaged with an equal weighting. This approach is used, as different resolution levels localize

Fig. 3
An estimate for the global posterior optimizer k obtained for the source space S k is found using the estimate for k−1 as the initial guess (" Randomized Scanning" section). We consider this approach necessary in order to maintain the estimates as close to the global optimum as possible. Namely, using a fixed initial guess might mean that the global optimizer is not found for some of the source spaces. Left: The global posterior optimizer is found for the posterior of space 1 (solid contours). Right: For space 2 (dashed contours), it is found (solid grey path), if the final estimate obtained in the case one 1 is used as the initial guess for 2 (grey circled point), while a local optimizer is obtained (dashed grey path) with the original initial guess 1 (black circled point) resulting in an optimization error. The global optimum might, in practice, correspond to a situation in which both a superficial and deep source are detected, while the deep activity might be undetected at the local one 1 3 different details ("Coarse-to-Fine Optimization" section). Consequently, the details found for the most levels are likely to gain the highest intensity in the final reconstruction. A schematic illustration of the resulting data flow has been included in Fig. 4.

Numerical Implementation with Zeffiro Interface
The forward and inverse solvers applied in this study were implemented in the Matlab (The MathWorks Inc.) as a part of the Zeffiro Interface (ZI) code package which is openly available in GitHub 1 . ZI is a tool enabling finite element (FE) based forward and inverse computations in electromagnetic brain applications. The forward approach of ZI together with the basic version of the IAS source reconstruction approach have been validated numerically in (Miinalainen et al. 2019;Pursiainen 2012b). ZI generates a uniform tetrahedral finite element (FE) mesh. Each source distribution is obtained by picking the first K entries of the randomly (uniformly) permuted set of the tetrahedron centers for the brain compartment. Due to the uniform mesh structure, this strategy leads to an evenly distributed set of source points. . ZI allows performing the source reconstruction routines using either a CPU or a GPU (graphics processing unit) type processor. Today, effective GPUs are available in power PCs an workstations but most laptops are still limited to CPU processing. Therefore, to compare the performance difference between GPU and CPU platforms, the computing time for forming a random set of multiresolution decompositions and inverting a given measurement data vector were evaluated for NVIDIA Quadro P6000 workstation GPU and Intel i7 5650U laptop CPU.

Numerical Experiments
In the numerical experiments, we used the realistic population head model 2 (PHM) (Lee et al. 2016), consisting of five layers (white matter, grey matter, cerebrospinal fluid (CSF), skull, and skin) and the three-layer Ary model in which concentric 87, 92 and 100 mm spheres present grey matter, skull and scalp layer. The cerebellum and vetricle layers included in the PHM were modeled as part of the grey matter and CSF, respectively. The conductivity of each layer can be found in Table 1. The PHM and Ary model were discretized with a uniform point lattice with the resolution 0.85 and 1 mm, leading to 24M and 30M tetrahedral elements and 4M and 5M nodes, respectively. In both cases, a single lead field matrix was generated for 10000 randomly distributed synthetic source positions. The lead field matrix entries were evaluated in SI units, i.e., Ohm/m and 1/m 2 for EEG and MEG, respectively. Each point contained three sources oriented along x-, y-and z-directions. Since the grey matter compartment of PHM does not include the thalamus, the source space was extended to cover both the white and grey matter compartment. Note that the lead field matrix and the corresponding source space have to be generated only once after which the space can be decomposed in multiple ways, e.g., different resolutions, as is the case in the proposed RAMUS process.

Simulated Measurements
For the Ary model, a total of 102 sensor points were distributed over the upper hemisphere. Using those, both electrode and radial magnetometer measurements of the electromagnetic field were simulated as shown in Fig. 5. The magnetometer locations were obtained by scaling the radial component of the source locations by a factor of 1.2. The electrodes were modeled using the complete electrode model . The inner and outer radius of the ring were 5 and 10 mm, respectively. The average contact resistance of each electrode was assumed to be 1 kOhm. In Fig. 4 A schematic visualization of the data flow during the reconstruction process for the multiresolution decompositions 1, 2, … , D each one with resolution levels 1, 2, … , L . The final estimate (6) obtained for the decomposition k is used as the level-one initial guess for the decomposition k + 1 . This sequential strategy for selecting the initial guess aims to minimize the effect of the optimization errors as suggested in Fig. 3. Note that with a good enough initial guess the global optimum is always found, meaning that the differences between the optimization results can be associated with the discretization errors which are modeled here as independently and identically distributed random variables the case of PHM, an EEG cap with 72 ring electrodes (10 mm outer and 5 mm inner diameter, 1 kOhm resistance) was attached on the head model. Two current dipoles were placed in shallow and deep parts of the grey matter. The source locations can be found in Table 2. Physiologically these could be interpreted as the somatosensory (superficial) P20/N20 and thalamic (deep) P14/N14 component, i.e., the 20 and 14 ms post-stimulus peaks. Activity for both locations occurs at the same time in the SEP/F response to the median nerve stimulus (Buchner et al. 1988(Buchner et al. , 1994a(Buchner et al. , 1995. When active simultaneously, the deep source was assumed to be slightly stronger in magnitude compared to the superficial one to enable the visibility of the deep part. This situation occurs momentarily in the median nerve stimulation, since the thalamic source obtains its maximum before the somatosensory activity increases in magnitude. As the measurement error term, we used zero mean Gaussian white noise with standard deviation of 3% respect to maximal signal amplitude. To investigate the noiserobustness of the source reconstruction, 5% noise was used in a single test. For the generality of the results, the maximum data entry of each dataset was normalized to one. The accuracy of the source recovery was analyzed in two 60 mm diameter spherical ROIs centered at the source locations (Fig. 5).

IAS MAP Iteration
The previous experience shows that, in order to distinguish deep activity (Calvetti et al. 2009), the hyperparameter values for the hyperprior have to be set as small as possible without risking the numerical stability of the reconstruction process. In the present study, the scale parameter 0 was chosen to be 1E-10 and the shape parameter was given the smallest possible value 1.5. These values were found to work generally well and they are supported also by the earlier studies (Calvetti et al. 2009). Ten iteration steps were performed to obtain a MAP estimate for a single resolution level. A single step was utilized in a single test.

Validation Tests
We analyzed the performance of the RAMUS reconstruction approach both visually and numerically in the tests (A)-(I) using the Ary model. The spherical domain was used in order to optimize the clarity of the results. In addition to these reconstructions, one test (J) was performed using PHM, i.e., the realistic model. The specifications for (A)-(J) can be found in Table 3.
The accuracy obtained in the cases (A)-(I) was analyzed by comparing the average position (center of mass), orientation and magnitude of the reconstructed distribution within the ROI to that of the actual dipole source. These average estimates were obtained with respect to the final reconstructed distribution of 10000 sources in each case (A)-(I) and for both single and multiple resolution reconstructions. In addition, the relative magnitude (between 0 and 1) of the distribution was calculated for each ROI. The source was classified as detected, if the relative maximum exceeded the value 0.1, and otherwise undetected. This threshold criterion was chosen as it represents roughly the limit of a visually detectable source. In (A)-(I), we varied the number of multiresolution decompositions, sparsity factor, hyperprior, the source magnitudes, and the measurement modality (EEG or E/MEG). When combining the lead field matrices for E/ MEG, the MEG lead field matrix and data was scaled so that the Frobenius norm, i.e., the 2-norm of all the entries, was equal to that of the EEG lead field matrix.
The case (I), was studied using three alternative approaches in addition to the basic multiresolution scheme. In the first one of these, the noise level was increased to 5%.

3
The second one involved only the coarse resolution level with otherwise unchanged parameters. In the last one, only single IAS MAP iteration was performed on each reconstruction level, meaning that the estimate obtained coincided with MNE (Calvetti et al. 2009). A total of 100-400 source positions at the coarse level, i.e., a number roughly comparable to that of the data entries ("Coarse-to-Fine Optimization" section), was found to work appropriately in the detection of the deep activity. When the sparsity factor between s = 8 and s = 5 , the source position count was within this interval at the coarsest level of a three-level multiresolution decomposition for the initial set of 10000 source positions. At the coarsest level, each source position was associated to about s 2 ( s L−1 with L = 3 ), i.e., between 64 and 25 finest-level source positions, respectively. The number D of multiresolution decompositions was chosen to be comparable to this number, slightly below or above that, in order to guarantee sufficient averaging over all the possible random basis choices.

Results
The results obtained in the numerical experiments have been included in Tables 4, 5 and Figs. 6,7,8,9,10 and 11. In each case, the deep and superficial component have been analyzed separately. Histograms for the cases (A)-(I) illustrate the accuracy of the reconstructed source with respect to the source position (mm), orientation (deg), amplitude, and the relative maximum of the current density within the ROI (Figs. 6, 7, 8 and 9) with respect to the global maximum.   The last one of these is utilized as a measure for the distinguishability of the source within the ROI. Examples of the reconstructions in the cases (A)-(I) are illustrated in Fig. 10, and the distributions obtained for (J), i.e., the realistic PHM, Fig. 6 The results of the deep source localization for the numerical experiments (A-I) conducted in the spherical domain (Table 3). The distributions of the position (mm), angle (°) and relative logarithmic (log10) amplitude difference to the exact dipole source, computed in the ROI, have been analyzed as histograms. The sample size is 50. Each reconstruction in the sample has been obtained by reconstructing the activity in the whole brain for an independent random realization of the noise vector and associating the total integrated activity in each ROI to the corresponding (deep or superficial) dipole source. Additionally, the histogram of the relative maximum in the ROI is given. The solid vertical line shows the median for each distribution, and the dashed lines mark the 90% confidence interval. In general, the results show that the IG hyperprior is necessary for detecting the deep source. The accuracy and reliability of the results increase along with the number of multiresolution decompositions. Furthermore, using E/MEG instead of EEG increased the accuracy of the deep source localization, while EEG was advantageous with respect to the amplitude of the deep source. The results are not visualized for the cases in which the localization criterion (relative maximum > 0.05) was satisfied by less than 5% of the reconstructions Fig. 7 The results of the superficial source localization for the numerical experiments (A-I) conducted in the spherical domain (Table 3).
In contrast to the case of the deep source, the superficial one is detected accurately in each case where its amplitude differs from zero. The most accurate results were obtained, when the deep source was absent. E/MEG yielded superior result compared to EEG 1 3 are shown in Fig. 11. The additional cases evaluated for (I), are presented in Fig. 8 and 9. The histograms in Figs. 6, 7, 8 and 9 illustrate the numerical accuracy of the RAMUS reconstruction approach. Case (A) suggests that the activity in both superficial and deep areas can be reconstructed in EEG, when applying IG as hyperprior. In (A), the superficial source is found with the median positioning accuracy of 8 mm, angle difference of 4.5° and logarithmic (log10) relative amplitude error of − 0.25, i.e., the amplitude of the reconstructed source is 56% compared to that of the actual one. For the deep source these errors are 15 mm, 12 deg, -0.65 (22% amplitude), respectively. Furthermore, as shown by the relative maximum, the superficial source always maximizes the (global) reconstruction, and the relative maximum within the deep ROI is around 50% of the global one in median.
Based on (B) and (C), it is obvious that the reconstruction accuracy is better, if only one of the two sources is active. Furthermore, increasing the intensity of the superficial source decreases the reconstruction accuracy for deep one which is shown by the case of (D) for which the median position, orientation amplitude, and relative maximum for the deep source are 18 mm, 17°, − 0.85 and 0.25, respectively. That is, the accuracy is lower than in (A). In (E), a sparsity factor of 5 was used instead of 8, meaning that the resolution difference between the subsequent levels was less steep, resulting in a weaker distinguishability of the deep source. The same observation was made in the case (F) in which 20 randomized decompositions instead of 100 were used. The deep activity was absent in (G), where we used the G hyperprior instead of IG, confirming the necessity of IG as the hyperprior. In (H) and (I), the use of the E/MEG lead field was observed to improve the deep localization accuracy around 7 mm and orientation accuracy about 8° with respect to the corresponding cases (A) and (D) of EEG data, while the superficial localization accuracy was practically unchanged for (H) and deviated less than 2 mm and 1 deg for (I). The results for E/MEG were visually more focal than the ones obtained with EEG (Fig. 10).
In the three additional tests performed with the parameter setting (I), the increased 5% noise level led to 5 mm and 3 deg lower positioning and orientation accuracy for the deep source, and a smaller 1 mm and 1 deg deviation for the superficial one. In the case of the coarse-level MAP iteration with 3% noise, 2 mm and 1 deg position and orientation improvement was observed for the deep source. For the superficial one, there was a 2 mm deviation in the position, while the orientation accuracy remained unchanged. The coarse-level estimate was visually less focal compared to the ones obtained with multiple resolution levels. MNE detected only the superficial source for which 1 mm position deflection and 3 deg orientation improvement were obtained compared to the basic case (I).
In (J), simultaneous localization of the simulated radial thalamic and tangential somatosensory component was found to be feasible with the realistic PHM model (Fig. 11). Similar to the spherical case, the deep activity had a lower amplitude than the superficial one. In the somatosensory area, with the physiological normal constraint, i.e., the assumption that the primary current is oriented along the inward surface normal, the activity was localized in the posterior wall of the central sulcus similar to the synthetic P20/ N20 component used in generating the data.
The results concerning the computing times have been included in Table 4. Those show that a superior performance was obtained with GPU processing which provided

Discussion
The present numerical results suggest that via the proposed randomized multiresolution scanning (RAMUS) technique one can obtain a robust and accurate MAP estimate for the primary current density in both superficial and deep parts of the brain. RAMUS was observed to enhance the visibility of deep components and also marginalizing the effect of the discretization without a remarkable computation cost. The noise-robustness of RAMUS was shown for 3% and 5% noise levels. As expected, the effect of the noise was observed to be the most obvious with respect to the deep source.
Utilizing a multiresolution approach was found to be crucial per se for the reconstruction quality, since maximal achievable accuracy for the deep components is significantly lower than for the superficial one. Detecting the deep source necessitated the presence of a coarse resolution level in the MAP estimation process, i.e., a sparsity factor s larger than one. The superior results were obtained with s = 8 . Decreasing the value of s, i.e., increasing the source count, quickly diminished the detectability of the deep component which can be observed based on the results obtained for s = 5 . The distinguishability of the deep source in the final estimate was determined by the number of the source positions at the coarsest level which, in this study, was observed to be around 100-400 roughly matching the sparsity factors between s = 8 and s = 5 . Investigating this interval was motivated by the fact that the maximal number of the detectable sources in the numerical system is determined by the number of the data entries ("Coarse-to-Fine Optimization" section) which is 102 for EEG and 204 for E/MEG, i.e., roughly of the same magnitude. In practice, the optimal size of the coarse system should also take the physiological modeling aspects into account and might be, therefore, also considerably larger than the present choice. For example, if the neural activity is limited to a priori known ROIs a larger number might be well-motivated. A comparison between the  (Table 3). In each image, the actual source and the center of mass of the reconstruction w.r.t. a ROI centered at the actual source position, are marked by the red and purple arrow, respectively. The case (I), was studied using the following three alternative approaches in addition to the basic multiresolution scheme. MNE: only single IAS MAP iteration was performed on each reconstruction level, meaning that the estimate obtained coincided with MNE. Noise 5%: the noise level was increased to 5%. Coarse-level MAP: only the coarse resolution level was applied in the MAP estimation process with otherwise unchanged parameters ▸ single (coarse-level) and multiple resolution results showed that the refinement of the resolution during the reconstruction process improves the focality of the reconstruction and its accuracy in the superficial areas. Nevertheless, the coarse-level reconstruction was marginally superior in the deep part, emphasizing that here the finer resolution levels slightly affected the coarser level outcome, which is here presumably optimal for the weakly distinguishable deep Fig. 11 The reconstruction (I) of the primary current density for the numerically modeled deep (thalamic P14/N14) and superficial (somatosensory P20/ N20) activity obtained using the population head model (PHM). On each row, the left column shows the amplitude and the right one the normal component in the direction of the surface normal. On 3rd and 4th row, the normal activity has been constrained into the outward direction. In each image, the actual source and the center of mass of the reconstruction w.r.t. a ROI centered at the actual source position, are marked by the red and purple arrow, respectively activity. Thus, it is important to adjust the decomposition parameters appropriately. The marginalization of the discretization errors via random scanning was perceived to be vital in order to optimize the robustness of the reconstruction which was observed to grow along the number of the multiresolution decompositions utilized.
When coupled with the iterative alternating sequential (IAS) algorithm, RAMUS constitutes essentially a repetitive MAP optimization process for HBM (O'Hagan and Forster 2004;Calvetti et al. 2009). Marginalizing the result over a given number of random multiresolution decompositions can be associated with computing an equal number of MAP estimates. Since the computational cost of the IAS algorithm is largely determined by the product between the lead field matrix and a candidate solution which is parallelized effectively in both CPU and GPU processors. Here, the latter option was found to achieve the fastest performance with the total computation time for a single reconstruction being 14 seconds which would be feasible in processing a larger dataset. Overall, the computational effort of evaluating the MAP via the RAMUS technique can be regarded as moderate compared to a full MCMC sampling based conditional mean (CM) estimate for the posterior density which has been evaluated in (Calvetti et al. 2009) within a ROI. Namely, achieving a full convergence of MCMC would require thousands of iterations (Liu 2001) and the effort of one iteration step is comparable to a single step of IAS. Thus, MCMC would be a slower option. Even though an optimization method, RAMUS can be also interpreted as a surrogate for CM, as it, on one hand, increases the robustness of the source reconstruction via sampling, but, on the other hand, does not provide as extensive information about the posterior density itself as an actual Bayesian sampler does.
The results obtained suggest that the IG hyperprior (O'Hagan and Forster 2004) is necessary in conjunction with IAS, when it is coupled with RAMUS, as the deep activity was not detected with G. Since here the cases of the G hyperprior and single-step MAP can be associated with the 1-norm regularized MCE and MNE (Uutela et al. 1999;Hämäläinen et al. 1993) ("Hierarchical Bayesian Model" section), respectively, it also seems that RAMUS would not enable correcting a depth bias related to either of these estimates. Previously, in (Calvetti et al. 2009), IG was found to perform well for the deep part, when a region of interest was used. Based on the present results, RAMUS provides the means to utilize the advantage of the IG within the whole brain and with a high imaging resolution, while maintaining the computational cost on an appropriate level.
We emphasize that the present conditionally Gaussian prior, in its current formulation, is depth, resolution and decomposition invariant. That is, additional physiological or operator based weighting or prior conditioning (Homa et al. 2013;Calvetti et al. 2015Calvetti et al. , 2018 is not necessary in order to balance the depth performance of the MAP estimate. Our interpretation for this is that RAMUS can correct the depth localization inaccuracies that are otherwise found with MAP estimates, as it, via the multiresolution approach, decomposes the source space into a set of a visible and invisible fluctuations, explores both sets, and also helps to marginalize the random numerical discretization and optimization errors out of the final estimate. Central here are the visibility of the deep activity at low resolution levels (Pascual-Marqui 1999), the concept of the sensitivity decomposition (Liu et al. 1995) and forming such through projections and multiresolution decompositions which have been investigated in the context of other inverse problems, e.g., in (Piana and Bertero 1997a;Pursiainen 2008).
In addition to the investigated properties, the choice for the scale parameter was also observed to be important in order to guarantee proper function of RAMUS. In each MAP estimation process, the present value 1E−10 was found to work well in detecting activity for both the spherical and realistic head model. The workable range for the scale parameter was observed to be from 1E−10 to 1E−08 similar to the previous findings (Calvetti et al. 2009). Outside this interval the deep activity was not found appropriately or the orientation accuracy of the estimates was lost. In the latter case, the estimate was locked into the direction of a Cartesian coordinate axis, meaning that, due to overly strong focality condition, only single component in the estimate differed from zero in the end of the iteration. Locking was also observed for MAP optimization sequences considerably longer than 10 iteration steps. With a sufficiently large scale parameter there was no locking, but the reconstructions obtained were also less focal.
The results of this article concern only the present numerical framework in which a deep and superficial source were detected simultaneously. Future work will include testing and analyzing the performance of the RAMUS approach with real experimental SEP/F data with the goal to distinguish cortical and sub-cortical activity, e.g., the P14/N14 (deep) and P20/N20 (superficial) components occuring in the stimulation of the median nerve. A comparison with other inverse methods capable of deep localization, such as LORETA and Beamforming (Pascual-Marqui et al. 2002Jonmohamadi et al. 2014), will also be important. Further method development topics will include a deeper investigation on the inverse effects of the hyperprior and decomposition parameters as well as finding alternative strategies to update the initial guess for the IAS MAP estimation technique. In the latter case, for instance, an averaged initial guess obtained with respect to several multiresolution decompositions might provide a potential alternative for the current approach which relies on a single decomposition. To make the random scanning computationally more efficient a solver based on parallel scanning processes might be developed. We also emphasize that RAMUS with its current formulation, the proposed algorithm can be applied to reduce discretization errors not only with the present IAS MAP method but potentially for a variety of source reconstruction techniques.