Optimal design of on‐scalp electromagnetic sensor arrays for brain source localisation

Abstract Optically pumped magnetometers (OPMs) are quickly widening the scopes of noninvasive neurophysiological imaging. The possibility of placing these magnetic field sensors on the scalp allows not only to acquire signals from people in movement, but also to reduce the distance between the sensors and the brain, with a consequent gain in the signal‐to‐noise ratio. These advantages make the technique particularly attractive to characterise sources of brain activity in demanding populations, such as children and patients with epilepsy. However, the technology is currently in an early stage, presenting new design challenges around the optimal sensor arrangement and their complementarity with other techniques as electroencephalography (EEG). In this article, we present an optimal array design strategy focussed on minimising the brain source localisation error. The methodology is based on the Cramér‐Rao bound, which provides lower error bounds on the estimation of source parameters regardless of the algorithm used. We utilise this framework to compare whole head OPM arrays with commercially available electro/magnetoencephalography (E/MEG) systems for localising brain signal generators. In addition, we study the complementarity between EEG and OPM‐based MEG, and design optimal whole head systems based on OPMs only and a combination of OPMs and EEG electrodes for characterising deep and superficial sources alike. Finally, we show the usefulness of the approach to find the nearly optimal sensor positions minimising the estimation error bound in a given cortical region when a limited number of OPMs are available. This is of special interest for maximising the performance of small scale systems to ad hoc neurophysiological experiments, a common situation arising in most OPM labs.


| INTRODUCTION
Optically pumped magnetometers (OPMs) are revolutionising the way we measure magnetic brain signals. Since their introduction in the brain imaging field in 2018 (Boto et al., 2018), they have provided insights into cerebral activity difficultly achieved by the wellestablished SQUID-based magnetoencephalography (MEG) sensors.
The reasons are simple: OPMs are portable and can be attached directly on the scalp. This allows (i) to reduce the distance between the sensors and the brain current generators by approximately 15 mm, with a consequent gain in signal-to-noise ratio (SNR) and (ii) the subject to move more freely during the acquisition, expanding the spectrum of MEG experiments to those previously achievable only by electroencephalography (EEG). Moreover, they do not require cryogenic cooling for functioning, making their running costs negligible. All these advantages are making OPM technology particularly attractive in a number of studies, from epilepsy (Vivekananda et al., 2020) to virtual reality .
As is the case with any new technique, scientists have been interested in quantifying the benefits over its predecessor. With regard to OPM-MEG, this means to compare their performance with conventional SQUID-MEG for characterising sources of brain activity. This task is of particular importance at this early developmental stage for determining optimal array designs that allow minimal source reconstruction errors. Such a comparison has been assessed through a number of metrics, for example, the SNR (Hill et al., 2020), the free energy (Duque-Munoz et al., 2019), the total information (TI) conveyed by the array (Iivanainen, Stenroos, & Parkkonen, 2017), and the spatial information density, a generalisation of the TI (Riaz, Pfeiffer, & Schneiderman, 2017). However, none of these works have studied how the OPM arrangement impacts directly on the variance of the reconstructed source considering arbitrary arrays and source estimation algorithms, which can consequently inform optimal sensor positioning systems.
In this article, we investigate the advantages of OPM-MEG systems by means of the Cramér-Rao bound (CRB). This metric provides a lower bound on the variance of the estimated source parameters achievable by any unbiased estimator, that is, regardless of the algorithm used (as long as it is unbiased). Furthermore, the CRB is an asymptotically tight bound, achievable by certain algorithms if the number of samples is large enough (Van Trees & Bell, 2013). These properties have made the CRB a standard tool in performance analysis and engineering design in a variety of fields, ranging from sonar to diffusion MRI (Alexander, 2008;Van Trees, 2002), and E/MEG in particular (Beltrachini, Von Ellenrieder, & Muravchik, 2011;Dogandzic & Nehorai, 2004;Hochwald & Nehorai, 1997;Muravchik & Nehorai, 2001;von Ellenrieder, Muravchik, & Nehorai, 2005). Here, we computed the CRB to quantify the impact of instrument and modelling noise in the variance of the position of the localised source, and employed it in three relevant problems. First, to assess the performance of wholehead standard OPM-based MEG arrays based on existing on-scalp sensor positioning setups. Although such large arrays are currently scarce (Hill et al., 2020), the CRB allows to quantify the improvements on source localisation compared with SQUID-based montages resembling commercially available systems, as well as with high density EEG (hdEEG). Secondly, we exploited the complementarity found between EEG and OPM-MEG to design optimal OPM and hybrid OPM/EEG arrays minimising the localisation error throughout the brain. This represents a major engineering problem in the area to which the CRB presents suitable and optimal solutions. Finally, we utilised the CRB for solving the more practical problem of finding optimal sensor positions when a specific brain region is targeted. We show that the CRB is a robust tool to design the optimal array for any given number of sensors and particular application, enhancing the flexibility of available equipment.

| The Cramér-Rao bound
The main objective in electromagnetic brain mapping is to characterise sources of electrical activity based on a set of measurements and a head and sensor models. Under certain assumptions, such sources can be expressed mathematically by a function of parameters representing their position and strength. Unfortunately, the estimation of these parameters from real data, usually referred to as the inverse problem in E/MEG, does not have a unique solution. This has led the community to present a myriad of methods based on different hypotheses and simplifications, and consequently leading to different results (Grech et al., 2008). For this reason, a comprehensive analysis on the impact of sensor type and position on source estimates would demand an immense effort.
The CRB solves this problem by establishing the optimal performance achievable by any unbiased estimator regardless of the algorithm used (Van Trees & Bell, 2013). Although the CRB is not the tightest deterministic bound (a property attributable to the Barankin bound [Barankin, 1949]), it has been shown to provide tight results achievable by existing algorithms (e.g., Beltrachini et al., 2011;Beltrachini et al., 2013;Fernández-Corazza, Beltrachini, von Ellenrieder, & Muravchik, 2013). Let ϑ ℝ N be a vector with elements ϑ i , i ¼ 1, …, N, defining the source of brain activity (both intensity and location). As we aim to characterise the endogenous current source from a given set of measurements, we look to estimate ϑ. In this context, the CRB states that where the inequality means that the difference between matrices on the left and right hand sides is positive semidefinite. In particular, the diagonal elements of CRB ϑ ð Þ establish the minimum variance on each source parameter ϑ i .
The calculation of the CRB depends on a number of ingredients, including the source and signal models, the sensor type and positions, and the numerical methodology employed to discretise the E/MEG governing equations (also known as E/MEG forward problem).
Here, we based our analysis on Muravchik and Nehorai's study (Muravchik & Nehorai, 2001), which assumed unconstrained dipolar sources of electrical activity and electromagnetic signals corrupted by zero-mean Gaussian noise with standard deviation σ b and σ e for MEG and EEG, respectively, representing instrument noise. By unconstrained sources we mean current generators whose locations are not assumed restricted to a surface. In this case, the parameter vector describing the source contains six elements, the first three defining the source position, and rest describing its moment. As for the numerical solution of the forward problem, we used the boundary element method (BEM).
We direct the interested reader to (Muravchik & Nehorai, 2001) for further details.
Once the CRB is computed for any given current source, we can use it to calculate meaningful error metrics related to the minimum variance achievable of the estimated position. In particular, we can obtain the volume of the 95% probability concentration ellipsoid, defined as the volume of the ellipsoid that encloses the mean of an unbiased estimate with 95% probability, and given by where c ¼ 2:8 and R is the 3 Â 3 submatrix of the CRB corresponding to the source location variables (Muravchik & Nehorai, 2001;Van Trees & Bell, 2013). For simplicity, we define the equivalent uncertainty radius, r 95 , as the radius of the sphere with volume equal to V 95 , Since r 95 defines the minimum volume containing the estimated source with 95% probability, we aim it to be as small as possible. In other words, the smaller r 95 is, the better the source location estimate can be.

| Head model and discretisation
We generated a head model based on the Colin27 high resolution MRI segmentation of the Montreal Neurological Institute (Aubert-Broche, Evans, & Collins, 2006). Images were employed to obtain three layers representing the scalp, skull, and brain. The innermost layer was tessellated in 40,960 triangular elements to avoid numerical errors for eccentric sources (Haueisen et al., 1997), whereas the other two consisted on 5,120 triangles. In the case of MEG, we employed the single shell model (Hämäläinen & Sarvas, 1989) and used linear elements (Ferguson, Zhang, & Stroink, 1994). In the case of EEG, BEM matrices were computed using the three layered model, linear elements (de Munck & Peters, 1993), and the isolated skull approach (Meijs et al., 1989). The adopted electrical conductivities were 0.011 S/m for the skull and 0.33 S/m for the scalp and brain (McCann, Pisano, & Beltrachini, 2019). Thirty-three thousand dipolar sources were uniformly located on each lobe's pial surface as obtained with FreeSurfer (Dale, Fischl, & Sereno, 1999). It is worth noting that sources were assumed unconstrained to any surface for their estimation, and their localisation on the pial surface was done just for clarity in the presentation of results.

| Sensor arrays
We considered five standard acquisition setups, all shown in Figure 1.
In the case of SQUID-based sensors, we utilised two commonly found 64 and 160 montages allowed to study the performance of medium and high density whole-head arrays that may be financially and ergonomically feasible with the OPM technology as it stands (the minimum distance between sensors was larger than 1.6 cm for ABC 160 in the Colin27 model). The ABC 160 system was also used to generate a discrete set of predefined locations where to place the on-scalp magnetometers to maximise the sensitivity to a particular brain region, as well as a basis for the hybrid OPM/EEG array (see Section 2.4). OPMs were modelled as magnetometers measuring the magnetic field in the F I G U R E 1 Sensor setups employed in the study. From left to right: 275 axial gradiometers (only coils closest to scalp are shown), 102 axial magnetometers and 204 planar gradiometers, ABC 64, 160 and 256 montages. Red dots indicate gradiometer coils, whereas yellow dots represent magnetometers or EEG electrodes (whenever corresponds) axial and one tangential direction (randomly chosen in this study; see Section 4 and Figure S1) and located at a distance of 6 mm from the scalp, emulating the second generation sensors provided by QuSpin Inc. (Colorado). The ABC 256 system was used to calculate the CRB for hdEEG, which has been shown to provide results competitive with MEG (e.g., Hedrich et al., 2017). Sensors were located utilising fiducial markers and employing instructions available in the vendor's website (as in the case of the ABC layouts). Point sensor and electrode models were used.

| Experiments
We conducted three experiments based on the CRB framework.
Firstly, we used it to compare the performance of whole-head OPMbased MEG systems with existing SQUID-MEG and hdEEG arrays.
This was done by calculating the equivalent uncertainty radius for sources located on the cortex and the arrays introduced in Figure 1.

Such a comparison is important for quantifying the pros and cons of
OPMs over established systems from a source localisation perspective, presenting the basis for engineering optimal prototypes.
In the second experiment, we employed the CRB to design opti- Finally, in the third experiment, we found the optimal arrangement for an arbitrary number of sensors minimising the lower bound on the source localisation variance for brain current generators in specific cortical regions. This is of particular interest for most MEG In all experiments, we adopted the MEG noise standard deviation,  On the other hand, all MEG arrays (regardless of the sensor type) showed, compared to hdEEG, better performance for eccentric and tangentially oriented sources, but worst performance for radially oriented sources (to which MEG is known for having lower sensitivity to) and generators located in deep brain regions.
In Figure 2b-d, we display a set of metrics that provide a quantitative means for comparing the performance of such acquisition systems. Figure

| Design of optimal whole-head systems
Based on the previous results, we further evaluated the complementarity between EEG and OPM systems for source localisation considering the same electrodes/sensors layouts. Figure 3a

| Optimal ad hoc arrays
Results to the third experiment are presented in Figure 5 for both pri-  Although hypothetical, the OPM ABC 160 concept is not difficult to envisage, with arrays composed of almost 50 sensors being constructed only 2 years after their introduction in the field (Hill et al., 2020).
It has been known for long time that EEG and MEG are not competing but complementary imaging techniques due to their dissimilar spatial sensitivities (e.g., Hari & Puce, 2017;Hunold et al., 2016  Results from the first experiment allowed to realise the importance of integrating EEG and OPMs for generating optimal portable whole-head systems. This led us to design two nearly optimal on-scalp electromagnetic arrays sensitive to signals generated throughout the brain: one built from OPMs (as it is currently the gold standard [Hill et al., 2020]) and another complemented with EEG electrodes. As expected, the more OPM sensors used, the better results the OPM array achieved. However, this was not the case for the hybrid system, for which we found a compromise between the number of magnetometers and electrodes for reaching the minimum of both mode and maximum of the equivalent uncertainty radius. This optimal relation is clearly dependent on the layout utilised, and resulted in 100 OPMs and 60 EEG electrodes for the ABC 160 setup. The results shown in The latter factor was found much more difficult to reduce to the level of the commercial array, requiring 47 and more than 100 to reach the same value as the MEG 275 GRAD system for the motor and temporal regions, respectively. Then, although OPM arrays can generate results comparable to SQUID MEG with only a limited number of sensors, one must remain cautious regarding the actual limitations of ad hoc systems. Even so, these findings are expected to impact on a number of studies centred on particular brain regions, such as brain computer interface (Clerc, Bougrain, & Lotte, 2016), focal cortical dysplasia (Vivekananda et al., 2020), and emotion processing (Ibáñez et al., 2011).
The utilisation of EEG positioning setups to constrain the location of the OPMs presents a number of benefits, such as the use of a common naming reference that aid in the dissemination and reproducibility of the results, as well as the consideration of a discrete and low dimensional space where to look for the optimal solution. Moreover, since all sensors need a corresponding holder, this option seems more practical than assuming free sensor positioning. Here, we used the CRB to find nearly optimal sensor locations within the ABC 160 layout for sources placed in cortical regions by sequentially adding the OPM that minimised the overall metric. Consequently, the resulting positions are not optimal by definition, since the algorithm assumes that any sensor arrangement will be also selected for any larger optimal configuration. Nevertheless, this resulted in a quicker implementation, which is needed for building an accessible tool with manageable computing requirements. Even more, this approach seems better suited for practical applications in which the number of sensors available at any moment may vary due to technological glitches, as we have personally experienced with state-of-the-art OPMs. The use of a more thorough methodology, such as combinatorial multilevel optimisation (Eichardt et al., 2019), is planned for future development.
Additionally, the use of EEG montages as the base setup for arbi- task to make them a reality. This led us to use the ABC 160 system, whose inter electrode/sensor distance allows for a full OPM array.
Additional computations ( Figure S1) have shown that the effect of rotating the tangential axis is negligible for whole-head systems, in line with other studies (Eichardt et al., 2019). This flexibility in the sensor orientation facilitates the assemble of high density OPM arrays, and therefore worth exploiting. Besides, the supplementary figure highlights the theoretical advantages of measuring the complete vector field at each sensor, which is now being explored in some research groups (Labyt et al., 2019). Notwithstanding, it is important to acknowledge constraints other than the space availability and sensor closeness, such as those related to sensor heating, weight, and crosstalk. In this regard, equidistant sensor systems used for dry EEG may present a valuable alternative (Fiedler et al., 2015).
It is worth noting that the CRB framework employed here does not come without limitations. First, as already discussed, the signal model assumed uncorrelated Gaussian noise representing mostly instrument and (up to some degree) geometrical modelling perturbations. Therefore, the methodology did not incorporate variability in the estimated parameters due to other error sources such as background activity de Munck, Vijn, & Lopes da Silva, 1992), sensor mislocation (Beltrachini et al., 2011) and tilting (Hill et al., 2020). All these perturbations are likely to increase the CRB metrics and consequently worth exploring in more detail. What is more, although this noise model is widely adopted and accepted, it is an idealised version of the real noise found in OPMs (and sensors in general), and therefore more work is needed to represent noise depending on the sensor technology (Eichardt et al., 2019;Oelsner et al., 2019). Second, we utilised the point electrode and sensor models, which are generally adopted in E/MEG, and in OPM research the omission of the CSF compartment (Piastra et al., 2021). Third, we modelled sources of brain activity as unconstrained current dipoles.
The use of such a simple model is required in the CRB context to avoid complicated expressions that may result difficult (or even impossible) to derive analytically. Consequently, this analysis gives an idea of the maximum sensitivity that can be expected from OPM and hybrid arrays for a single source. The separability and effect of crosstalk of multiple sources, which could impede to achieve the sensitivity bound employed here, is not addressed by the presented methodology. Further work should focus on analysing this particular aspect of the arrays obtained, which could be noticeable on practical designs.
Several source models were introduced to describe source generators more realistically, as those based on the multipolar expansion (Beltrachini, 2019;Jerbi, Mosher, Baillet, & Leahy, 2002). Nevertheless, dipolar current generators are generally the first choice for brain source characterisation, most prominently in epilepsy (Rampp et al., 2019;Vivekananda et al., 2020), as well as the basis for most E/MEG-IP algorithms. Lastly, the algorithm for constructing the hybrid array may be optimised in several ways. One option is to incorporate variable weighting for the MEG signals depending on the source direction, resulting minimal for radial and deep sources, and maximal for tangential. This may lead to a more distributed pattern between electrodes and sensors. Another is to incorporate both EEG and MEG simultaneously in the optimisation, differently to the work here presented where OPMs were prioritised. Further work will analyse the impact that these changes may have in the resulting arrays. Open access funding enabled and organized by Projekt DEAL.