NeuroImage 277 (2023) 120257

Contents lists available at ScienceDirect

NeuroImage
journal homepage: www.elsevier.com/locate/ynimg

Sensor array design of optically pumped magnetometers for accurately
estimating source currents
Yusuke Takeda a,b,∗, Tomohiro Gomi c, Ryu Umebayashi c, Sadamu Tomita c, Keita Suzuki b,
Nobuo Hiroe b, Jiro Saikawa c, Tatsuya Munaka c, Okito Yamashita a,baComputational Brain Dynamics Team, RIKEN Center for Advanced Intelligence Project, 2-2-2 Hikaridai, Seika-cho, Soraku-gun, Kyoto, 619-0288, Japan
Department of Computational Brain Imaging, ATR Neural Information Analysis Laboratories, 2-2-2 Hikaridai, Seika-cho, Soraku-gun, Kyoto, 619-0288, Japanc
Technology Research Laboratory, Shimadzu Corporation, 3-9-4 Hikaridai, Seika-cho, Soraku-gun, Kyoto 619-0237, Japanba r t i c l e

i n f oKeywords:
Optically pumped magnetometer (OPM)
Magnetoencephalography (MEG)
Sensor array
Source imaging
Resolution matrix

a b s t r a c t
An optically pumped magnetometer (OPM) is a new generation of magnetoencephalography (MEG) devices thatis small, light, and works at room temperature. Due to these characteristics, OPMs enable ﬂexible and wearable
MEG systems. On the other hand, if we have a limited number of OPM sensors, we need to carefully design their
sensor arrays depending on our purposes and regions of interests (ROIs). In this study, we propose a method thatdesigns OPM sensor arrays for accurately estimating the cortical currents at the ROIs. Based on the resolution
matrix of minimum norm estimate (MNE), our method sequentially determines the position of each sensor to
optimize its inverse ﬁlter pointing to the ROIs and suppressing the signal leakage from the other areas. We callthis method the Sensor array Optimization based on Resolution Matrix (SORM). We conducted simple and realistic
simulation tests to evaluate its characteristics and eﬃcacy for real OPM-MEG data. SORM designed the sensor
arrays so that their leadﬁeld matrices had high eﬀective ranks as well as high sensitivities to ROIs. Although
SORM is based on MNE, the sensor arrays designed by SORM were eﬀective not only when we estimated thecortical currents by MNE but also when we did so by other methods. With real OPM-MEG data we conﬁrmed its
validity for real data. These analyses suggest that SORM is especially useful when we want to accurately estimate
ROIs’ activities with a limited number of OPM sensors, such as brain-machine interfaces and diagnosing brain
diseases.

1. Introduction
Our brain activities produce very weak magnetic ﬁelds outside ofour heads. By measuring these magnetic ﬁelds, we can examine brain
activities with excellent temporal resolution on the millisecond order.
Conventionally, the magnetic ﬁelds produced by our brains have beenmeasured by magnetoencephalography (MEG) systems equipped with
Superconducting QUantum Interference Devices (SQUIDs), such as the
Elekta Neuromag Vectorview 306 system (Helsinki). These systems arelarge, placing more than 200 sensors over the whole head, and liquid
helium is housed around the sensors to cool them.
Recently, a new device was introduced (Shah and Wakai, 2013), an
optically pumped magnetometer (OPM), which has attracted great interest from the MEG community (Boto et al., 2016; 2017; Brookes et al.,2022; Ilmoniemi and Sarvas, 2019; Tierney et al., 2019). An OPM sensor is small, light, and works at room temperature without cryogenic
cooling. Due to these characteristics, OPMs enable ﬂexible and wearable MEG systems, which allow us to measure MEGs under such severe∗

situations that include the brain activities of children and during movements (Seymour et al., 2021). Furthermore, we can put OPM sensors on
the head, yielding high sensitivity to brain activities (Boto et al., 2017;
Iivanainen et al., 2017).
On the other hand, because this technology is still in its early stage,well-established OPM-MEG measurement systems are currently unavailable. After buying OPM sensors, we ourselves need to construct an MEG
measurement system, such as a coil system for nulling the environmental magnetic ﬁeld (Holmes et al., 2018; 2019; Iivanainen et al., 2019)
and a cap or helmet to hold the sensors (Hill et al., 2020). If we only
have a few OPM sensors (a common case for many OPM users), we mustdesign the sensor array carefully depending on our purposes and the regions of interests (ROIs). For example, if we want to extract the brain
activity at the motor cortex with ﬁve OPM sensors for a brain-computerinterface (BCI), the sensors should be placed around the motor cortex.
However, the accuracy of this expectation is unclear. Even if it is true,
we do not know at which positions the sensors should be put.

Corresponding author.
E-mail address: takeda@atr.jp (Y. Takeda).

https://doi.org/10.1016/j.neuroimage.2023.120257.
Received 8 November 2022; Received in revised form 21 April 2023; Accepted 28 June 2023
Available online 29 June 2023.
1053-8119/© 2023 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license(http://creativecommons.org/licenses/by-nc-nd/4.0/)

Y. Takeda, T. Gomi, R. Umebayashi et al.

NeuroImage 277 (2023) 120257

Recently, Beltrachini et al. (2021) proposed a method to design OPM
sensor arrays for minimizing source localization error. It sequentially determines the sensor positions to minimize the Cramér-Rao bound of thesource localization error. However, because their scheme was developed
for dipole methods (Muravchik et al., 2000), which approximate brain
activities by a small number of current dipoles (Aine et al., 2000; Hari,
1991; Mosher et al., 1992), it cannot be used for distributed source methods, such as Minimum Norm Estimates (MNE) (Hämäläinen et al., 1993;
Hämäläinen and Ilmoniemi, 1994; Ilmoniemi and Sarvas, 2019) or a Linearly Constrained Minimum Variance (LCMV) beamformer (Ilmoniemi
and Sarvas, 2019; Van Veen et al., 1997). Distributed source methods
assume distributed currents in the brain and are widely used in the MEG
literature.
In this study, we propose a method that designs OPM sensor arrays
for accurately estimating distributed currents. From a set of possible
sensor positions, our method selects some positions so that the estimated cortical currents have accurate time-series at the ROIs. To achieve
this aim, we focus on the MNE’s resolution matrix (Hauk et al., 2022;
Grave de Peralta-Menendez and Gonzalez-Andino, 1998), which represents the relationship between the true and estimated currents. The
larger the diagonal elements in the resolution matrix are, the more accurate the time-series of the estimated current becomes. Therefore, we
derived a recurrence formula of the resolution matrix and based on it,
sequentially selected sensor positions so that the resolution matrix haslarge diagonal elements at the ROIs. We call this method a Sensor array
Optimization based on Resolution Matrix (SORM).
We conducted simple and realistic simulation tests to evaluate its
characteristics and eﬃcacy for real OPM-MEG data. SORM selected sensor positions so that their leadﬁeld matrices had high eﬀective ranks
as well as high sensitivities to ROIs. The sensor arrays designed by
SORM were eﬀective not only when we estimated the cortical currents
by MNE but also when we did so by other methods: an LCMV beamformer (Van Veen et al., 1997) and the hierarchical Variational Bayesian
method (hVB) (Sato et al., 2004; Yoshioka et al., 2008). Furthermore,
we applied SORM to real OPM-MEG data to conﬁrm its validity for realdata. Based on these analyses, we believe that SORM is especially usefulfor accurately estimating ROIs’ activities with a limited number of OPM
sensors, such as BCI and the diagnosis of brain diseases.

Fig. 1. Relationship between sensor array and signal leakage. (A) Assumed sensor arrays. Red circles indicate sensor positions. (B) ROI (a vertex in left somatosensory cortex). (C) Signal leakage into ROI from other areas. For eachsensor array in A, we calculated MNE’s resolution matrix and showed its absolute row vector values at ROI. Values over 0.3 of their maximum value areshown.

where 𝑳 is an 𝑁 × 𝑀 matrix, called an inverse ﬁlter. In the case of MNE,
𝑳 is calculated by
𝑳 = [𝑮T 𝑮 + 𝜆𝑰 ]−1 𝑮T ,where 𝜆 is a regularization constant and 𝑰 is the 𝑁 × 𝑁 identity matrix.
By substituting Eq. (1) for 𝑩 in Eq. (2), we obtain𝑱̂ = 𝑳𝑮𝑱 + 𝑳𝑬 .

𝑹 = 𝑳𝑮(5)

and taking the expectation of 𝑬 , Eq. (4) is rewritten asE[𝑱̂ ] = 𝑹𝑱 .
𝑹 is an 𝑁 × 𝑁 matrix, called a resolution matrix (Hauk et al., 2022;
Grave de Peralta-Menendez and Gonzalez-Andino, 1998), which represents the relationship between the true and estimated currents.
If 𝑹 is an identity matrix, the current is perfectly estimated. However, 𝑹 is strongly degenerate because its rank is 𝑀 << 𝑁, and therefore
it cannot be an identity matrix with rank 𝑁. The non-diagonal elementsof 𝑹 represent the signal leakage that occurs when estimating the current (Brookes et al., 2012; Colclough et al., 2015; Palva et al., 2018;
Sato et al., 2018). Speciﬁcally, 𝑹(𝑎, 𝑏), where 𝑏 ≠ 𝑎, represents the signal leakage from vertex 𝑏 to 𝑎. Therefore, to accurately estimate thecurrent at vertex 𝑎, 𝑹(𝑎, 𝑎) needs to be larger than 𝑹(𝑎, 𝑏).
In the case of MNE, the resolution matrix is calculated from
Eqs. (3) and (5) by

In this section, we propose our SORM algorithm.
2.1. Assumption and purpose
We assume that we know the ROIs and that the possible OPM sensor
locations have already been determined. Hereafter, for brevity, the set of
possible OPM sensor locations is called the candidates. Given these ROIs
and candidates, SORM selects the sensor positions so that the estimatedcortical currents at the ROIs have accurate time-series.

𝑹 = [𝑮T 𝑮 + 𝜆𝑰 ]−1 𝑮T 𝑮.

2.2. Basic strategy(6)This is a symmetric matrix; that is, the signal leakage from vertex 𝑏 to𝑎 and from 𝑎 to 𝑏 is identical. Because this equation does not include
MEG data, we do not need to conduct MEG experiments to calculate theresolution matrices of MNE.
Figure 1 shows sample resolution matrices. We assumed two sensor
arrays (Fig. 1A, red circles). For each, we calculated the resolution matrix of MNE by Eq. (6) (see Section 3 for detailed descriptions). We set
ROI to a vertex in the left somatosensory cortex (Fig. 1B). To examine
the signal leakage into the ROI from other areas, we show the absolutevalues of the row vector at the ROI in the resolution matrix (Fig. 1C):a so-called cross-talk function (Hauk et al., 2022). The signal leakage
largely changed depending on the sensor arrays.
These considerations and results lead to SORM’s basic strategy: determine the sensor positions so that the resolution matrix of MNE has
larger diagonal elements than the non-diagonal elements at the ROIs.

For this purpose, we focus on the resolution matrix (Hauk et al.,2022; Grave de Peralta-Menendez and Gonzalez-Andino, 1998) of MNE(Hämäläinen et al., 1993; Hämäläinen and Ilmoniemi, 1994). We explain this by starting from the basics of cortical current estimation.
In distributed source methods, MEG data are expressed as(1)

where 𝑩 is the 𝑀 × 𝑇 MEG data matrix, 𝑀 is the number of channels,
𝑇 is the number of time points, 𝑮 is the 𝑀 × 𝑁 leadﬁeld matrix, 𝑁 is
the number of current dipoles (vertices), 𝑱 is the 𝑁 × 𝑇 cortical current
matrix, and 𝑬 is the 𝑀 × 𝑇 noise matrix.
Cortical current 𝑱 is estimated by𝑱̂ = 𝑳𝑩 ,(4)

By putting

2. Designing sensor arrays based on resolution matrix

𝑩 = 𝑮𝑱 + 𝑬 ,

(3)

(2)
2

Y. Takeda, T. Gomi, R. Umebayashi et al.

NeuroImage 277 (2023) 120257

2.3. Procedure to select sensor positions

3.1.2. Simulating brain activities
For simulating brain activities, we downloaded the 16 subjects’ T1 images from the multi-subject, multi-modal neuroimaging
dataset (OpenNEURO ds000117-v1.0.1) created by Wakeman and Henson (2015). From these images, we made brain models, which deﬁne
10,003 vertices on their cortical surfaces.
For the motor task, we selected two vertices from the left and rightprimary motor cortices. These selected vertices are referred to as sourcesin the simulation tests. We assigned currents at the sources and spatiallysmoothed them using a Gaussian ﬁlter whose standard deviation (SD)was 6 mm. For the auditory task, we selected two vertices, which arealso referred to as sources in the simulation tests, from the left and rightprimary auditory cortices, assigned currents at the sources, and spatiallysmoothed them using the Gaussian ﬁlter. To highlight the eﬀect of thesignal leakage, the time-series of the currents were generated to be orthogonal between the sources so that the signal leakage between themaﬀected their estimation accuracy. The amplitudes of the simulated currents were determined based on that of a current estimated from actual
SQUID-MEG data.

To identify the relationship between the sensor positions and MNE’s
resolution matrix, we derived how the resolution matrix changes byadding a new sensor. Here we assume that each sensor measures a magnetic ﬁeld along a single axis; each sensor has a single channel.
Suppose that we already set some sensors whose leadﬁeld and resolution matrices of MNE are 𝑮 and 𝑹. By adding a new sensor, 𝑹 ischanged to
𝑹+ = 𝑹 +𝜆𝑪 −1 𝒉T 𝒉𝑪 −1,
1 + 𝒉𝑪 −1 𝒉T(7)where 𝑪 = 𝑮T 𝑮 + 𝜆𝑰 and 𝒉 is a new row of the leadﬁeld matrix corresponding to the new sensor [see Appendix A for the derivation of
Eq. (7)]. Hereafter, a row of a leadﬁeld matrix is called a sensitivity vector because it represents a sensor’s sensitivity to each vertex. By deﬁning
a gain vector
𝒙 = 𝑪 −1 𝒉T = [𝑮T 𝑮 + 𝜆𝑰 ]−1 𝒉T ,(8)

Eq. (7) is rewritten as
𝑹+ = 𝑹 +3.1.3. Regions of interest
We set the ROIs based on the cortical parcellations deﬁned in the
HCP MMP 1.0 (Glasser et al., 2016). For the motor task, we set them tothe primary motor cortices that are labeled “L_4_ROI” and “R_4_ROI” in
the HCP MMP 1.0. For the auditory task, we set the ROIs to the early auditory cortices that are labeled “L_A1_ROI”, “L_MBelt_ROI”, “R_A1_ROI”,
and “R_MBelt_ROI” in the HCP MMP 1.0.

𝜆𝒙𝒙T.
1 + 𝒙T 𝑪 𝒙

Because 𝜆∕(1 + 𝒙T 𝑪 𝒙) is scalar, we focus on 𝒙𝒙T , which is an 𝑁 × 𝑁
matrix. The power of its diagonal elements at ROIs can be calculatedby 𝒙(𝒏)T 𝒙(𝒏), where 𝒏 = [𝑛1 , ⋯ , 𝑛𝑣 ] represents the elements corresponding to ROIs. If 𝒙(𝒏)T 𝒙(𝒏) is large, the diagonal elements at ROIs in 𝑹+
increase.
Based on these results, we propose the following procedure that selects the sensor positions from their candidates (Algorithm 1):

3.1.4. Candidates of sensor positions
We set the candidates of the OPM sensor positions based on those of
the 102 magnetometers of an Elekta Neuromag Vectorview 306 system(Helsinki). From the multi-subject, multi-modal neuroimaging dataset,
we derived the sensor positions of the magnetometers relative to thebrains, moved them 12 mm away from the scalp surfaces, and used themas the candidates of the OPM sensor positions.

Algorithm 1 Select K sensor positions.
1: Set 𝜆 to a ﬁxed value
2: Initialize 𝑮T 𝑮 by 𝑮T 𝑮 = zeros(𝑁 , 𝑁 )3: for 𝑘 = 1 ∶ 𝐾

Calculate 𝒙 = [𝑮T 𝑮 + 𝜆𝑰 ]−1 𝒉T for each remaining sensor position
Select the sensor position that maximizes 𝒙(𝒏)T 𝒙(𝒏)5:6:Update 𝑮T 𝑮 by 𝑮T 𝑮 = 𝑮T 𝑮 + 𝒉T 𝒉
7: end4:

3.1.5. Making leadﬁeld matrices
For each subject, we made a leadﬁeld matrix of all the candidate
sensor positions. From the T1 image, we constructed a 1-shell (cerebrospinal ﬂuid [CSF]) head conductivity model. Based on it, we made aleadﬁeld matrix by solving the Maxwell equations with a boundary element method (BEM). In this study, we assumed one-dimensional current
dipoles normal to the cortical surface.

By default, regularization constant 𝜆 was set to𝜆0 =tr (𝑮Tall 𝑮all )𝑁

× 0.1,(9)3.1. Simple simulation test

3.1.6. Selecting sensor positions
We selected the sensor positions by three methods: uniform, leadﬁeld norm, and SORM. For all the methods, to select 𝐾 sensor positions,
we ranked all the candidate sensor positions and selected the ﬁrst 𝐾
positions.
In the uniform method, we manually ranked the sensor positions sothat the sensors were uniformly distributed over the brain. This rank
was ﬁxed regardless of the subjects and the tasks.
In the leadﬁeld norm method, we ranked the sensor positions basedon their sensitivities to the ROIs. From the leadﬁeld matrix of all the candidate sensor positions, we extracted the leadﬁeld vectors corresponding
to the ROIs, calculated their norms for the individual sensor positions,and ranked them in descending order of the norms.
For the SORM, we ranked the sensor positions by the order of theselected positions.

3.1.1. Simulated experiment
Our simulation posited that 16 healthy adults participated in an
OPM-MEG experiment. In it, they performed two tasks: motor and auditory, each of which consisted of 100 trials. During the tasks, we recorded
the magnetic ﬁelds normal to their scalp surfaces (Z-axis) with OPM sensors.

3.1.7. Simulating the OPM-MEG data
We set the number of OPM sensors to 1, ⋯ , 50. For each number, we
simulated OPM-MEG data by Eq. (1). Observation noise 𝑬 was madeusing Gaussian white noise whose SD was 2 × 10−13 T.
Figure 2 shows the simulated brain activities and sample OPM-MEG
data at all the candidate positions. The signal-to-noise ratios (SNRs) of

where 𝑮all is the leadﬁeld matrix of all the candidate sensor positions.
3. Material and methods
We examined SORM’s performance with two simulation tests: simple
and realistic. The simple simulation test characterizes SORM’s behavior, and the realistic simulation test evaluates its eﬃcacy for real OPMMEG data. Furthermore, we applied SORM to real OPM-MEG data. These
analyses were conducted using VBMEG (https://vbmeg.atr.jp/), which
is a MATLAB toolbox to estimate cortical currents from MEG and/or
electroencephalography (EEG) data (Takeda et al., 2019).

3

Y. Takeda, T. Gomi, R. Umebayashi et al.

NeuroImage 277 (2023) 120257

Fig. 2. Simulated brain activities and OPM-MEG data in simple simulation test. (A) Simulated brain activities. Top ﬁgures show amplitude distributions of simulated
cortical currents averaged across time. Values over 0.3 of their maximum value are shown. Bottom ﬁgures show time-series of simulated cortical currents in left andright hemispheres. (B) Stimulus-triggered average of sample OPM-MEG data at 102 candidate positions.

the OPM-MEG data were −8.40 ± 1.57 dB and −14.08 ± 1.05 dB (mean± SD) for the motor and auditory tasks. The SNRs were calculated by
10 log10 𝑃𝑠 ∕𝑃𝑛 , where 𝑃𝑠 and 𝑃𝑛 are the powers of 𝑮𝑱 and 𝑬 in Eq. (1),assuming that 102 OPM sensors were put at all the candidate positions.

In this study, LCMV beamformer indicates a scalar-type beamformer because we assumed a one-dimensional current dipole at each vertex.
hVB provides a framework to use the prior information on the current variance as soft constraints on the variance Sato et al. (2004);
Yoshioka et al. (2008). It can use a variety of prior information, suchas the uniform distribution, the fMRI activity, and a meta-analysis offMRI studies (Suzuki and Yamashita, 2021). We used the uniform distribution as prior information on the current variance to set a conﬁdence
parameter, “bayes_parm.prior_weight,” to 0.0001. This parameter controls the conﬁdence in the prior information relative to the amount ofdata samples (ranging from 0 to 1) (Takeda et al., 2019).
The estimated currents were averaged across the trials.

3.1.8. Estimating the cortical currents
From the simulated OPM-MEG data, we estimated the cortical currents by three methods: MNE (Hämäläinen et al., 1993; Hämäläinen andIlmoniemi, 1994), a LCMV beamformer (Van Veen et al., 1997), and hVB(Sato et al., 2004; Yoshioka et al., 2008).
MNE estimated the current by
𝑱̂ = [𝑮T 𝑮 + 𝜆MNE 𝑰 ]−1 𝑮T 𝑩 ,where

3.1.9. Evaluating the accuracy of estimated source currents
For each condition (task, current estimation method, sensor array,and subject), we quantiﬁed the time-series accuracy of the estimatedcurrents. At the sources, we calculated the correlation coeﬃcients between the true and estimated currents and averaged them between
hemispheres.

tr (𝑮T 𝑮)
𝜆MNE =
× 0.1.
𝑁
The LCMV beamformer estimated the current by
𝒋̂𝑛 = 𝒍𝑛 𝑩 ,where 𝒋̂𝑛 is the estimated current at vertex 𝑛 and 𝒍𝑛 is the 1 × 𝑀 inverse
ﬁlter. The inverse ﬁlter was calculated bỹ
𝒍𝑛 = [𝒈T𝑛 𝑫−1̃𝒈𝑛 ]−1 𝒈T𝑛 𝑫−1,

3.1.10. Statistical test
We compared the time-series accuracy of SORM with those of theuniform and leadﬁeld norm methods. For each task, sensor array, andcurrent estimation method, we performed a two-tailed Wilcoxon signedrank test under the null hypothesis where the medians of the correlationcoeﬃcients were identical between SORM and the uniform/leadﬁeld
norm methods.
This is the multiple comparison problem, which we solved by controlling the false discovery rate (FDR). FDR manages the expected proportion of false positive ﬁndings among all the rejected null hypotheses(Benjamini and Hochberg, 1995). We estimated the 𝑞-values by Storey
and Tibshirani’s method (2003). From the distribution of the 𝑝-values,
we ﬁrst estimated the proportion of null 𝑝-values 𝜋0 , and based on 𝜋0 weconverted the 𝑝-values to 𝑞-values. The FDRs were controlled at 0.01.

(10)

where
̃ = 𝑫 + 𝜆BF 𝑰 ,
𝑫
𝜆BF =

tr 𝑫
× 0.1.
𝑀

𝒈𝑛 is the 𝑀 × 1 leadﬁeld vector for vertex 𝑛. 𝑫 is the 𝑀 × 𝑀 covariancematrix computed between all the sensor pairs. We imposed regularization term 𝜆BF 𝑰 , assuming that in real data analyses covariance matrix 𝑫
often becomes rank deﬁcient due to the noise reduction by independent
component analysis (ICA) and a common average referencing for EEG.
4

Y. Takeda, T. Gomi, R. Umebayashi et al.

NeuroImage 277 (2023) 120257

3.1.11. Calculating contribution ratios in estimated source currents
To identify the factors responsible for the time-series accuracy, we
examined the contribution ratios in the estimated source currents. We
deﬁned three kinds of contribution ratios: target, non-target, and noise.
The target contribution ratio is that of the component originated fromthe left/right source in the estimated current at the left/right source;
it represents the purity of the estimated currents. The non-target contribution ratio is that of the component originated from the right/leftsource in the estimated current at the left/right source; it represents theamount of signal leakage into a source from another source. The noise
contribution ratio is that of the components originated from the noise[𝑬 in Eq. (1)].
To calculate the contribution ratios, we decomposed the simulated
OPM-MEG data into components originated from the left source, theright source, and the noise. We separately converted these components
into cortical currents using the inverse ﬁlters from Section 3.1.8. For theleft source, the target contribution ratio was calculated by
𝑐𝑟target =3.3. Realistic simulation test
We conducted a realistic simulation test in the same way as the simple one except for procedures to simulate the brain activities and the
OPM-MEG data. In this test, we simulated the former using the SQUIDMEG+EEG data during the auditory task and the latter using the SQUIDMEG+EEG data during the resting-state task and the empty-room OPM
data. These procedures are described below.
3.3.1. Simulating brain activities
To simulate realistic brain activities, we used the cortical currents
estimated from the SQUID-MEG+EEG data during the auditory task by
hVB. We extracted currents at the left and right auditory cortices froma subject’s cortical current. At the same sources as the simple simulation test, we assigned the currents and spatially smoothed them using
a Gaussian ﬁlter whose SD was 6 mm. The amplitudes of the simulated
currents were determined so that the reconstructed SQUID-MEG datafrom the simulated currents had the same amplitudes as the original
one.

𝑦𝑙
,
𝑦𝑙 + 𝑦𝑟 + 𝑦𝑛

where 𝑦𝑙 , 𝑦𝑟 , and 𝑦𝑛 are respectively the powers of the components originated from the left source, the right source, and the noise. The non-target
contribution ratio was calculated by
𝑐𝑟non-target =3.3.2. Simulating the OPM-MEG data
To simulate realistic OPM-MEG data, we used the cortical currents
estimated from the SQUID-MEG+EEG data during the resting-state task
and the empty-room OPM-MEG data. The OPM-MEG data was generatedby𝑦𝑟.
𝑦𝑙 + 𝑦𝑟 + 𝑦𝑛

The noise contribution ratio was calculated by

𝑩 = 𝑮𝑱 + 𝑮𝑱 𝑠 + 𝑬𝑦𝑛
𝑐𝑟noise =.
𝑦𝑙 + 𝑦𝑟 + 𝑦𝑛

where 𝑱 is the simulated brain activity produced in Section 3.3.1 and𝑱 𝑠 is the spontaneous brain activity. 𝑱 𝑠 was extracted from the corticalcurrents estimated from the SQUID-MEG+EEG data during the restingstate task by MNE. 𝑬 was generated from the empty-room OPM-MEG
data. We generated surrogate data which had the same amplitude and asimilar power spectrum with the preprocessed empty-room OPM-MEG
data and used it as 𝑬 .
Figure 3 shows the simulated brain activity and the sample OPMMEG data at all the candidate positions. The SNRs of the OPM-MEG
data were −18.52 ± 1.40 dB.

In the same way, the contribution ratios were also calculated for theright source. Then, they were averaged between the left and rightsources.
3.2. Analyzing real OPM-MEG, SQUID-MEG, and EEG data
For the realistic simulation test and applying SORM to real OPMMEG data, we created and analyzed the OPM-MEG, SQUID-MEG, andEEG (OSE) dataset (https://vbmeg.atr.jp/nictitaku209/), which contains four subjects’ OPM-MEG, SQUID-MEG, and EEG data during fourtasks: auditory, motor, somatosensory, and resting-state. It also contains
empty-room OPM-MEG data.
The four subjects participated in three recording sessions: OPMMEG, SQUID-MEG and EEG, and MRI. The OPM-MEG data wererecorded with 15 sensors (10 QZFM Gen-2 and ﬁve QZFM Gen-3, QuSpin Inc., U.S.), each of which recorded a magnetic ﬁeld along twoaxes (Y and Z) (Supplementary material, SFig. 1). The SQUID-MEG
and EEG data were simultaneously recorded with a whole-head 400channel system (210-channel Axial and 190-channel Planar Gradiometers; PQ1400RM; Yokogawa Electric Co., Japan) and a whole-head 63channel system (BrainAmp; Brain Products GmbH, Germany). Electrooculograms (EOGs) were also recorded simultaneously.
We preprocessed the OPM-MEG, SQUID-MEG, and EEG data for estimating the cortical currents. Because the OPM-MEG data contain large
ﬂuctuation below 8 Hz, we applied a lowpass ﬁlter at 40 Hz and a highpass ﬁlter at 8 Hz to all the data. Then we regressed out the EOG components from the data. For the EEG data, we took a common averagereference and made the averages of the EEG data across the channels to
0.
From T1 images, we constructed 1-shell (CSF) and 3-shell (CSF, skull,and scalp) head conductivity models for the MEG and EEG data. Basedon the models, we made leadﬁeld matrices by solving the Maxwell equations with the BEM. Using the leadﬁeld matrices, we estimated the cortical currents from the preprocessed OPM-MEG data. We also estimatedthe cortical currents from the SQUID-MEG and EEG (SQUID-MEG+EEG)
data (Takeda et al., 2019).

3.4. Applications to real OPM-MEG data
Finally, we evaluated the eﬃcacy of SORM using real OPM-MEG data
during the auditory task. The OPM-MEG data were recorded with 30
channels (15 sensors × 2 axes) (Supplementary material, SFig. 1). Here
we regarded the 30 channels as the candidates of pseudo sensors. Fromthem, we selected sensors by the uniform, the leadﬁeld norm, and the
SORM methods and estimated the cortical currents from the OPM-MEG
data of the selected sensors. We quantiﬁed the time-series accuracy ofthe estimated currents by comparing them with the cortical currents estimated from the identical subjects’ SQUID-MEG+EEG data during theauditory task. For each subject and vertex, we calculated the correlation coeﬃcient of the cortical currents between the OPM and SQUIDMEG+EEG and averaged them within the ROIs (the early auditory cortices). We did not conduct statistical tests because of the small sample
size (four subjects).
4. Results
4.1. Simple simulation test
To characterize the SORM behavior, we conducted a simple simulation test. We assumed that the 16 subjects performed two tasks: motor
and auditory. During them, magnetic ﬁelds normal to their scalp surfaces were recorded with OPM sensors. We set the ROIs to the primary
motor cortices for the motor task and the early auditory cortices for theauditory task (Fig. 4A).
5

Y. Takeda, T. Gomi, R. Umebayashi et al.

NeuroImage 277 (2023) 120257

Fig. 3. Simulated brain activity and OPM-MEG
data in realistic simulation test. (A) Simulated
brain activity. Top ﬁgures show amplitude distributions of simulated cortical current averaged
across time. Values over 0.3 of their maximum
value are shown. Bottom ﬁgures show timeseries of simulated cortical current in left andright hemispheres. (B) Stimulus-triggered average of sample OPM-MEG data at 102 candidate
positions.

Fig. 4. Sensor arrays designed by uniform, leadﬁeld
norm, and SORM with 10 sensors. (A) ROIs. (B) Sensor arrays. Red circles indicate selected sensor positions. A subject’s results are shown.

4.1.1. Designed sensor arrays
Figure 4B shows the sensor arrays designed by the uniform, leadﬁeld
norm, and SORM methods with 10 sensors. Compared with the uniform
method, SORM prioritized the positions around the ROIs. On the other
hand, compared with the leadﬁeld norm method, SORM prioritized positions far from the ROIs.
To characterize the behaviors of these three methods, for each sensor
array we examined the property of the leadﬁeld matrix by two indices:sensitivity to the ROIs and the eﬀective rank. To obtain sensitivity to
the ROIs, we extracted from the leadﬁeld matrix the submatrix corresponding to the ROIs and calculated its Frobenius norm. The eﬀective
rank represents the diversity of the sensitivity vectors across the selected
sensor positions. We estimated it based on the singular values of theleadﬁeld matrix (Roy and Vetterli, 2007).
Figure 5 shows the sensitivity to the ROIs and the eﬀective rank foreach sensor array. The uniform method showed high eﬀective ranks but
low sensitivities to the ROIs (Fig. 5, blue lines). In contrast, the leadﬁeld
norm method showed high sensitivities to the ROIs but low eﬀective
ranks (Fig. 5, red lines). SORM showed both high sensitivities and higheﬀective ranks (Fig. 5, green lines). These results indicate that SORM
selected sensor positions to increase both the sensitivity and the rank,
while the other methods only increased one of them.
To comprehend the rationale of this result, we expressed 𝒙T 𝒙 using
the singular value decomposition (SVD) (Appendix B). It suggests that
SORM tends to select a sensor position whose sensitivity vector has large
values at ROIs and is in the null space of the existing sensors’ leadﬁeld
matrix, heightening the eﬀective rank of the new sensor array. This suggestion is consistent with Fig. 5. These results are intuitively interpreted
as follows: a new sensor should be sensitive to the ROIs and represent
outside of them to eﬀectively suppress the signal leakage from the outside.

Fig. 5. Properties of leadﬁeld matrices. (A) Sensitivity to ROIs. (B) Eﬀective
rank. Thick lines and shaded areas represent averages and SDs of values acrosssubjects.

4.1.2. Accuracy of estimated source currents
For each sensor array, we simulated the OPM-MEG data and estimated the cortical current with MNE, the LCMV beamformer, and hVB.
We evaluated the time-series accuracy of the estimated currents by the6

Y. Takeda, T. Gomi, R. Umebayashi et al.

NeuroImage 277 (2023) 120257

Fig. 6. Time-series accuracy of estimated source currents in simple simulation test. For each subject, we
calculated correlation coeﬃcients between true and estimated currents at left and right sources and averaged
them between hemispheres. Thick lines and shaded
areas represent averages and SDs of correlation coefﬁcients across subjects. Asterisks indicate that correlation coeﬃcients of SORM were signiﬁcantly higher
than those of uniform and leadﬁeld norm methods(𝑞 < 0.01).

correlation coeﬃcients between the true and estimated currents at thesources.
Figure 6 shows the time-series accuracy of the estimated currents.
For all the current estimation methods, SORM outperformed the uniform
and leadﬁeld norm methods, especially when we used few sensors (< 15sensors). This result indicates the eﬃcacy of SORM for various current
estimation methods.
On the other hand, when we used MNE and SORM, the time-series
accuracy decreased as the number of sensors increased (Figs. 6A andD, green lines). This result indicates that using more sensors was notalways better for MNE.

4.2. Realistic simulation test

4.1.3. Contribution ratios in estimated source currents
To interpret the time-series accuracy, we calculated the contribution
ratios in the estimated source currents during the auditory task. Figure 7
shows the contribution ratios averaged across the subjects. The target
contribution ratios highly correlated with the time-series accuracy (theircorrelation coeﬃcients were 0.94 ± 0.03) (Fig. 7, blue areas). This result
suggests that the time-series accuracy mainly reﬂected the purity of theestimated currents.
When we used MNE, SORM drastically decreased the non-target contribution ratios as the number of sensors increased (Fig. 7G, red area).
This result indicates that SORM eﬀectively suppressed the signal leakage and increased the time-series accuracy, suggesting the validity of its
algorithm.
On the other hand, when we used MNE and SORM, the time-series
accuracy decreased as the number of sensors increased (Figs. 6A andD, green lines). The contribution ratios suggest that this was due to theincrease of the noise contribution ratios (Fig. 7G, green areas). As thenumber of sensors increased, SORM added sensor positions where thedata had low SNR. This increased the noise contribution ratios and decreased the time-series accuracy.
On the contrary, when we used the LCMV beamformer and hVB, thetime-series accuracy did not greatly decrease (Figs. 6B, C, E, and F).
The contribution ratios suggest that the LCMV beamformer suppressed
the noise contribution ratios (Figs. 7B, E, and H, green areas), hVB suppressed the non-target and noise contribution ratios (Figs. 7C, F, and I,
red and green areas), and as a result their time-series accuracy did not
greatly decrease.

4.3. Dependency on regularization constantTo evaluate the eﬃcacy of SORM for real OPM-MEG data, we conducted a realistic simulation test. In it, we used the real data to simulate
the brain activities and noise. We used the same sensor arrays from thesimple simulation test (Fig. 4B) because we set the ROIs to the same areas (the early auditory cortices) as the simple simulation test (Fig. 4A).
Figure 8 shows the time-series accuracy of the estimated currents.
SORM outperformed the other methods, especially when we used the
LCMV beamformer and just a few sensors (< 15 sensors). This tendency
held unless the SNR was too low (Supplementary material, SFig. 2), indicating the eﬃcacy of SORM for real OPM-MEG data.

SORM has a hyperparameter: the regularization constant [𝜆 inEq. (8)], which was set to 𝜆0 [Eq. (9)] by default. We examined thedependency of SORM on the regularization constant by changing theregularization constant to 100𝜆0 , 𝜆0 , and 0.01𝜆0 . For each constant, we
designed the sensor array for the auditory task by SORM, simulated the
OPM-MEG data during the auditory task, estimated the current by MNE,
and evaluated the time-series accuracy.
Figure 9 shows the sensor arrays and the time-series accuracy foreach regularization constant. They were almost the same across the constants, indicating that SORM did not largely depend on the regularization constant.
4.4. Application to real OPM-MEG data
Finally, we applied SORM to the real OPM-MEG data during the auditory task. We regarded the 30 OPM channels as the candidates of thepseudo sensors and from them selected some sensors (Fig. 10A). From
the OPM-MEG data of the selected sensors, we estimated the cortical currents (Fig. 10B). The estimated currents at the ROIs (the early auditory
cortices) were compared with those estimated from the same subjects’
SQUID-MEG+EEG data during the auditory task (Fig. 10C).
Figure 10D shows the similarity of the estimated currents. Even when
we used just a few pseudo sensors (< 10 sensors), SORM showed highsimilarity comparable with the leadﬁeld norm method. That is, evenwith few pseudo OPM sensors, the leadﬁeld norm method and SORM
7

Y. Takeda, T. Gomi, R. Umebayashi et al.

NeuroImage 277 (2023) 120257

Fig. 7. Contribution ratios in estimated source currents of auditory task in simple simulation test. They
were averaged across subjects.

Fig. 8. Time-series accuracy of estimated source currents in realistic simulation test. For each subject, we calculated correlation coefﬁcients between true and estimated currents at left and right sources
and averaged them between hemispheres. Thick lines and shaded areas represent averages and SDs of correlation coeﬃcients across subjects. Asterisks indicate that correlation coeﬃcients of SORM weresigniﬁcantly higher than those of uniform and leadﬁeld norm methods (𝑞 < 0.01).

enabled us to obtain the ROIs’ currents resembling those estimated froma massive amount of data (400-ch SQUID-MEG and 63-ch EEG data).
This result also indicates the eﬃcacy of SORM for real OPM-MEG data.

5.1. Methodological considerations
SORM sequentially selects the sensor positions from their candidates
so that the MNE’s resolution matrix has large diagonal elements at the
ROIs. This is a combinatorial optimization problem, which we solved bya greedy algorithm. For example, when we select the 10 optimal sensor
positions from the 102 candidates, we need to ﬁnd one solution from
13 combinations. Note that SORM does not ﬁnd the op102 C10 > 2 × 10
timal solution; it provides a sub-optimal one. By using the recurrence
formula of the resolution matrix [Eq. (7)], SORM approximately solves
this problem quickly. Indeed, it selected the 10 sensor positions from the
102 candidates within one minute using our computer [Intel(R) Xeon(R)CPU E5-2680 v2 @ 2.80GHz, 20 cores] and the selected positions outperformed those of the uniform and leadﬁeld norm methods (Figs. 6
and 8). Since the SORM algorithm is very easy to implement, SORM ispractical for a wide range of OPM users.
To accurately estimate ROIs’ activities, SORM focuses on the diagonal elements at them in the MNE’s resolution matrix and tries to increase
them. However, this strategy also increased the non-diagonal elements
and the diagonal elements (Supplementary material, SFig. 3), indicating that SORM enlarged the signal leakage into the ROIs as well as the

5. Discussion
In this study, we proposed a SORM algorithm to design an OPM sensor array for accurately estimating the activity of ROIs. It sequentiallyselects sensor positions so that the resolution matrix of MNE has large
diagonal elements at the ROIs. This strategy made the leadﬁeld matrix
have a high eﬀective rank as well as high sensitivity to the ROIs (Fig. 5,
green lines). Indeed, our mathematical consideration provided the rationale behind this result (Appendix B). Simulation tests showed that
the SORM yielded higher time-series accuracy than the other methods,
especially when we used a few sensors (< 15 sensors) (Figs. 6 and 8).
Although SORM is based on MNE, the sensor arrays designed by SORM
worked well regardless of the current estimation methods (Fig. 6). We
also applied SORM to real OPM-MEG data and conﬁrmed its eﬃcacy forreal data. These results suggest the eﬃcacy of SORM when we want toaccurately estimate ROIs’ activity with a limited number of OPM sensors.
8Y. Takeda, T. Gomi, R. Umebayashi et al.

NeuroImage 277 (2023) 120257

ROIs’ activities. Nevertheless, SORM worked well in our simulation tests
and applications to the OPM-MEG data (Figs. 6–10). As the number of
sensors increased, the diagonal elements increased faster than the nondiagonal elements at the ROIs (Supplementary material, SFig. 3). As aresult, the time-series of the estimated currents at the ROIs became moreaccurate as the number of sensors rose.
5.2. Spatial accuracy of estimated cortical currents
We focused on the temporal accuracy of the estimated cortical currents at ROIs and ignored their spatial accuracy. This is because SORM
assumes that we have ROIs, meaning that we already know the information of the source locations. In this case, since evaluating the spatial
accuracy is not very meaningful, we focused on the temporal accuracy
and optimized the SORM algorithm to increase it.
On the other hand, the spatial accuracy of estimated currents has
been widely evaluated in source imaging studies by localization accuracy and spatial extent (e.g. Hauk et al., 2022; Sekihara et al., 2005.
Following these studies, we also evaluated the spatial accuracy of theestimated currents. SORM showed comparable spatial accuracy with theuniform and leadﬁeld norm methods (Supplementary material, SFig. 4).
Fig. 9. Dependency of SORM on regularization constant. For each regularization constant, sensor array designed with 10 sensors (A) and time-series accuracyof source currents estimated by MNE (B) are shown for auditory task. Thick lines
and shaded areas represent averages and SDs of correlation coeﬃcients acrosssubjects.

5.3. Contribution ratios
To identify the factors responsible for the time-series accuracy of theestimated source currents, we examined the contribution ratios in theFig. 10. Application to real OPM-MEG
data. (A) Sensor array of 10 pseudo sensors designed by SORM. (B) Cortical current
estimated from OPM-MEG data at selected
pseudo sensors (A). (C) Cortical current estimated from SQUID-MEG+EEG data. In
B and C, cortical currents during auditory task estimated by hVB are shown.
Top ﬁgures show their amplitude distributions averaged across time. Values over 0.3
of their maximum value are shown. Bottom ﬁgures show time-series of currents at
ROIs (early auditory cortices). In A–C, results of a subject are shown. (D) Similarity of currents at ROIs between OPM andSQUID-MEG+EEG data. For each subject
and vertex, we calculated correlation coeﬃcient of currents between OPM and SQUIDMEG+EEG data and averaged them within
ROIs. Thick lines and shaded areas represent
averages and SDs of correlation coeﬃcients
across subjects.

9

Y. Takeda, T. Gomi, R. Umebayashi et al.

NeuroImage 277 (2023) 120257

Data and code availability

currents (Fig. 7). The contribution ratios suggest that the time-series
accuracy (Fig. 6) mainly reﬂected the purity of the estimated currents(Fig. 7, blue areas). They also indicate that SORM eﬀectively suppressed
the signal leakage in the currents estimated by MNE (Fig. 7G, red area)
and increased their time-series accuracy (Fig. 6D, green line), suggesting
the validity of the SORM algorithm.
On the other hand, the contribution ratios revealed diﬀerences acrossthe current estimation methods. For example, although the LCMV beamformer and hVB showed comparable time-series accuracy (Figs. 6E andF, green lines), the contribution ratios identiﬁed diﬀerent strategies. The
LCMV beamformer achieved high accuracy by suppressing the noise
contribution ratio (Fig. 7H, green area), and hVB achieved it by suppressing the signal leakage (Fig. 7I, red area). In this way, the contribution ratios highlighted the diﬀerences across the current estimation
methods and helped us select the best method for our purpose.

The data used in this study are publicly available from the multisubject, multi-modal neuroimaging dataset (OpenNEURO ds000117v1.0.1,
https://openneuro.org/datasets/ds000117/versions/1.0.1/)
and the OSE dataset (https://vbmeg.atr.jp/nictitaku209/). The
codes for performing SORM are available at https://github.com/
RIKEN- AIP- CBDteam/Takeda2023_SORM.
Declaration of Competing Interest
The authors declare no competing interests.
Credit authorship contribution statement
Yusuke Takeda: Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Software, Validation, Visualization,
Writing – original draft, Writing – review & editing. Tomohiro Gomi:
Formal analysis, Investigation, Methodology, Validation, Writing – review & editing. Ryu Umebayashi: Formal analysis, Investigation,
Methodology, Validation, Writing – review & editing. Sadamu Tomita:
Formal analysis, Investigation, Methodology, Supervision, Validation,
Writing – review & editing. Keita Suzuki: Data curation, Formal analysis, Investigation, Resources, Software, Writing – review & editing.
Nobuo Hiroe: Data curation, Investigation, Resources, Writing – review & editing. Jiro Saikawa: Project administration, Supervision, Validation, Writing – review & editing. Tatsuya Munaka: Project administration, Supervision, Validation, Writing – review & editing. Okito
Yamashita: Conceptualization, Formal analysis, Funding acquisition,
Project administration, Resources, Supervision, Writing – review & editing.

5.4. Extended usages
SORM uses MNE’s resolution matrix because it can be calculated
without MEG data [Eq. (6)]; that is, we do not need to conduct an MEG
experiment in advance to design the sensor arrays. On the other hand,
we can calculate the resolution matrices for other current estimation
methods, such as the LCMV beamformer. Such adaptive ﬁlter methods
need MEG data to calculate their resolution matrices. For example, the
LCMV beamformer needs a covariance matrix computed between all the
sensor pairs, and we need to estimate it by conducting a preparatory experiment or from existing data, such as open data. Once the covariancematrix is obtained, we can apply the same greedy algorithm as SORM
(Section 2.3) except for steps 4 and 5. At these steps, instead of calculating 𝒙(𝒏)T 𝒙(𝒏), we must calculate Eqs. (10) and (5) for each remaining
sensor position to ﬁnd the position that maximizes the diagonal elements
at the ROIs in the resolution matrix.
SORM assumes that each OPM sensor measures a magnetic ﬁeld
along a single axis. However, current commercial OPM sensors can simultaneously measure a magnetic ﬁeld along two or three axes. Indeed,
the OPM-MEG data of the OSE dataset were recorded with the dual-axial
sensors. With tri-axial OPM sensors, we can estimate 3D magnetic ﬁeld
vectors at the sensors. Furthermore, using tri-axial OPM sensors has theoretical advantages for current estimation (Brookes et al., 2021). SORM
can be easily extended for tri-axial OPM sensors. In this case, each sensor has three channels. Using SORM, we can rank all the channels andthe sensor positions based on the best channel in each sensor positionand select the top 𝐾 sensor positions.

Acknowledgments
This research was supported by a contract with the National Instituteof Information and Communications Technology entitled “Analysis ofmulti-modal brain measurement data and development of its applicationfor BMI open innovation” (209). It was also supported by JSPS KAKENHI(Grant Number 21K07269).
Appendix A
In this section, we derive Eq. (7), which represents how the resolution matrix of MNE changes by adding a new sensor. Here we assume
that each sensor has a single channel.
We use the following notations and equations:5.5. Applicability
Since SORM designs the sensor arrays for accurately estimating ROIs’
activity, it is suitable when we know the ROIs, such as BCI and diagnosing brain diseases. Recent proof-of-concept studies suggested the usefulness of OPM for non-invasive BCI (Wittevrongel et al., 2021) anddiagnosing epilepsy (Vivekananda et al., 2020). For such purposes, simple MEG measurement systems with a few OPM sensors are practically
helpful. For BCI, by setting the ROIs to the motor cortices, SORM candesign a sensor array to accurately estimate a subject’s motor intension.
For the diagnosis of brain diseases, by setting ROIs to such aﬀected regions as epileptic foci, SORM could design a sensor array to accurately
monitor their activities.
Since SORM only relies on leadﬁeld matrices [Eq. (8)], it can be applied to various problems where leadﬁeld matrices (or forward models)
exist, such as designing sensor arrays of EEG and near-infrared spectroscopy (NIRS). SORM can be used not only for the brain but also forthe heart, muscles, and peripheral nerves, etc. For example, it could design an eﬀective OPM sensor array for reconstructing cardiac activity
from a magnetocardiography.

1. 𝑀: number of existing channels (sensors)
2. 𝑁: number of vertices
3. 𝑮: 𝑀 × 𝑁 leadﬁeld matrix of existing sensors
4. 𝒉: 1 × 𝑁 sensitivity vector of a new sensor
5. 𝑮+ : concatenated matrix of 𝑮 and 𝒉 (𝑮+ = [𝑮; 𝒉] in MATLAB notation)
6. 𝑮T+ 𝑮+ = 𝑮T 𝑮 + 𝒉T 𝒉
7. 𝑪 = 𝑮T 𝑮 + 𝜆𝑰
−1 T 𝒉𝑪 −1
8. [𝑪 + 𝒉T 𝒉]−1 = 𝑪 −1 − 𝑪 𝒉 −1
[Derived from Eq. (11) in Faul and
1+𝒉𝑪 𝒉T
Tipping (2001)
9. 𝑹 = 𝑪 −1 𝑮T 𝑮: resolution matrix of existing sensors
Using the above notations and equations, the resolution matrix after
adding a new sensor can be written:[]−1
𝑹+ = 𝑮T+ 𝑮+ + 𝜆𝑰 𝑮T+ 𝑮+[]−1 ( T)
𝑮 𝑮 + 𝒉T 𝒉(Using the 6-th equation)= 𝑮T 𝑮 + 𝜆𝑰 + 𝒉T 𝒉
[
] ( T
)
T −1
T
= 𝑪 +𝒉 𝒉
𝑮 𝑮+𝒉 𝒉(Using the 7-th equation)10

Y. Takeda, T. Gomi, R. Umebayashi et al.

NeuroImage 277 (2023) 120257())
𝑪 −1 𝒉T 𝒉𝑪 −1 ( T
𝑮 𝑮 + 𝒉T 𝒉(Using the 8-th equation)
𝑪 −1 −−1 T
1 + 𝒉𝑪 𝒉
𝑪 −1 𝒉T 𝒉𝑪 −1 𝑮T 𝑮 𝑪 −1 𝒉T 𝒉𝑪 −1 𝒉T 𝒉
= 𝑪 −1 𝑮T 𝑮 + 𝑪 −1 𝒉T 𝒉 −−1 + 𝒉𝑪 −1 𝒉T
1 + 𝒉𝑪 −1 𝒉T
−1 T
−1
−1 T
𝑪
𝒉
𝒉
𝑪
𝑪
−
𝜆𝑰
(
)
𝑪 𝒉 𝒉𝑪 −1 𝒉T 𝒉
= 𝑹 + 𝑪 −1 𝒉T 𝒉 −−−1 T
1 + 𝒉𝑪 𝒉
1 + 𝒉𝑪 −1 𝒉T(Using the 9-th equation)

Beltrachini, L., von Ellenrieder, N., Eichardt, R., Haueisen, J., 2021. Optimal design of onscalp electromagnetic sensor arrays for brain source localisation. Hum. Brain Mapp.
42, 4869–4879. doi:10.1002/hbm.25586.
Benjamini, Y., Hochberg, Y., 1995. Controlling the false discovery rate: A practical andpowerful approach to multiple testing. J. R. Stat. Soc. Series B Stat. Methodol. 57,
289–300.
Boto, E., Bowtell, R., Krüger, P., Fromhold, T.M., Morris, P.G., Meyer, S.S., Barnes, G.R.,
Brookes, M.J., 2016. On the potential of a new generation of magnetometers for
MEG: a beamformer simulation study. PLoS One 11, e0157655. doi:10.1371/journal.pone.0157655.
Boto, E., Meyer, S.S., Shah, V., Alem, O., Knappe, S., Kruger, P., Fromhold, T.M.,
Lim, M., Glover, P.M., Morris, P.G., Bowtell, R., Barnes, G.R., Brookes, M.J.,
2017. A new generation of magnetoencephalography: room temperature measurements using optically-pumped magnetometers. Neuroimage 149, 404–414.
doi:10.1016/j.neuroimage.2017.01.034.
Brookes, M.J., Boto, E., Rea, M., Shah, V., Osborne, J., Holmes, N., Hill, R.M.,
Leggett, J., Rhodes, N., Bowtell, R., 2021. Theoretical advantages of a triaxial optically
pumped magnetometer magnetoencephalography system. Neuroimage 236, 118025.
doi:10.1016/j.neuroimage.2021.118025.
Brookes, M.J., Leggett, J., Rea, M., Hill, R.M., Holmes, N., Boto, E., Bowtell, R.,
2022. Magnetoencephalography with optically pumped magnetometers (OPM-MEG):
the next generation of functional neuroimaging. Trends Neurosci. 45, 621–634.
doi:10.1016/j.tins.2022.05.008.
Brookes, M.J., Woolrich, M.W., Barnes, G.R., 2012. Measuring functional connectivity in
MEG: a multivariate approach insensitive to linear source leakage. Neuroimage 63,
910–920. doi:10.1016/j.neuroimage.2012.03.048.
Colclough, G.L., Brookes, M.J., Smith, S.M., Woolrich, M.W., 2015. A symmetric multivariate leakage correction for MEG connectomes. Neuroimage 117, 439–448.
doi:10.1016/j.neuroimage.2015.03.071.
Faul, A.C., Tipping, M.E., 2001. Analysis of sparse bayesian learning. Advances in Neural Information Processing Systems 14. https://proceedings.
neurips.cc/paper/2001/ﬁle/02b1be0d48924c327124732726097157-Paper.pdf
Glasser, M.F., Coalson, T.S., Robinson, E.C., Hacker, C.D., Harwell, J., Yacoub, E., Ugurbil, K., Andersson, J., Beckmann, C.F., Jenkinson, M., Smith, S.M., Van Essen, D.C.,
2016. A multi-modal parcellation of human cerebral cortex. Nature 536, 171–178.
doi:10.1038/nature18933.
Hämäläinen, M., Hari, R., Ilmoniemi, R.J., Knuutila, J., Lounasmaa, O.V., 1993. Magnetoencephalography - theory, instrumentation, and applications to noninvasive studiesof the working human brain. Rev. Modern Phys. 65, 413–497.
Hämäläinen, M., Ilmoniemi, R.J., 1994. Interpreting magnetic ﬁelds of the brain: minimum norm estimates. Med. Biol. Eng. Comput. 32, 35–42. doi:10.1007/bf02512476.
Hari, R., 1991. On brain’s magnetic responses to sensory stimuli. J. Clin. Neurophysiol. 8,
157–169. doi:10.1097/00004691-199104000-00004.
Hauk, O., Stenroos, M., Treder, M.S., 2022. Towards an objective evaluation of
EEG/MEG source estimation methods - the linear approach. Neuroimage 255, 119177.
doi:10.1016/j.neuroimage.2022.119177.
Hill, R.M., Boto, E., Rea, M., Holmes, N., Leggett, J., Coles, L.A., Papastavroum, M.,
Everton, S.K., Hunt, B.A.E., Sims, D., Osborne, J., Shah, V., Bowtell, R.,
Brookes, M.J., 2020. Multi-channel whole-head OPM-MEG: helmet design
and a comparison with a conventional system. Neuroimage 219, 116995.
doi:10.1016/j.neuroimage.2020.116995.
Holmes, N., Leggett, J., Boto, E., Roberts, G., Hill, R.M., Tierney, T.M., Shah, V.,
Barnes, G.R., Brookes, M.J., Bowtell, R., 2018. A bi-planar coil system for nullingbackground magnetic ﬁelds in scalp mounted magnetoencephalography. Neuroimage
181, 760–774. doi:10.1016/j.neuroimage.2018.07.028.
Holmes, N., Tierney, T.M., Leggett, J., Boto, E., Mellor, S., Roberts, G., Hill, R.M., Shah, V.,
Barnes, G.R., Brookes, M.J., Bowtell, R., 2019. Balanced, bi-planar magnetic ﬁeld andﬁeld gradient coils for ﬁeld compensation in wearable magnetoencephalography. Sci.
Rep. 9, 14196. doi:10.1038/s41598-019-50697-w.
Iivanainen, J., Stenroos, M., Parkkonen, L., 2017. Measuring MEG closer tothe brain: performance of on-scalp sensor arrays. Neuroimage 147, 542–553.
doi:10.1016/j.neuroimage.2016.12.048.
Iivanainen, J., Zetter, R., Grön, M., Hakkarainen, K., Parkkonen, L., 2019. On-scalp MEG
system utilizing an actively shielded array of optically-pumped magnetometers. Neuroimage 194, 244–258. doi:10.1016/j.neuroimage.2019.03.022.
Ilmoniemi, R.J., Sarvas, J., 2019. Brain signals: physics and mathematics of MEG and EEG.
The MIT Press.
Mosher, J.C., Lewis, P.S., Leahy, R.M., 1992. Multiple dipole modelling and localization from spatio-temporal MEG data. IEEE Trans. Biomed. Eng. 39, 541–557.
doi:10.1109/10.141192.
Muravchik, C., Bria, O., Nehorai, A., 2000. EEG/MEG error bounds for a dynamic dipole
source with a realistic head model. Methods Inf. Med. 39, 110–113.
Palva, J.M., Wang, S.H., Palva, S., Zhigalov, A., Monto, S., Brookes, M.J., Schoffelen, J.M., Jerbi, K., 2018. Ghost interactions in MEG/EEG source space: anote of caution on inter-areal coupling measures. Neuroimage 173, 632–643.
doi:10.1016/j.neuroimage.2018.02.032.
Grave de Peralta-Menendez, R., Gonzalez-Andino, S.L., 1998. A critical analysis of linear
inverse solutions to the neuroelectromagnetic inverse problem. IEEE Trans. Biomed.
Eng. 45, 440–448. doi:10.1109/10.664200.
Roy, O., Vetterli, M., 2007. The eﬀective rank: a measure of eﬀective dimensionality. EURASIP. https://www.eurasip.org/Proceedings/Eusipco/Eusipco2007/Papers/
a5p-h05.pdf
Sato, M., Yamashita, O., Sato, M., Miyawaki, Y., 2018. 10.1371/journal.pone.0198806,
Information spreading by a combination of MEG source estimation and multivariate
pattern classiﬁcation. PLoS One 13, e0198806.

== 𝑹 + 𝑪 −1 𝒉T 𝒉 −= 𝑹 + 𝑪 −1 𝒉T 𝒉 −= 𝑹+

𝑪 −1 𝒉T 𝒉 − 𝜆𝑪 −1 𝒉T 𝒉𝑪 −1 + 𝒉𝑪 −1 𝒉T 𝑪 −1 𝒉T 𝒉
1 + 𝒉𝑪 −1 𝒉T()
1 + 𝒉𝑪 −1 𝒉T 𝑪 −1 𝒉T 𝒉 − 𝜆𝑪 −1 𝒉T 𝒉𝑪 −1
1 + 𝒉𝑪 −1 𝒉T

𝜆𝑪 𝒉 𝒉𝑪.
1 + 𝒉𝑪 −1 𝒉T
−1 T

−1

Appendix B
To understand which kind of sensor position maximizes 𝒙T 𝒙, we express it by the singular vectors of the leadﬁeld matrix of the existing
sensors. We use the same notations and equations as in Appendix A. For
simplicity, we set the ROIs to the whole brain.
SVD decomposes the leadﬁeld matrix of existing sensors 𝑮 as
𝑮 = 𝑼 𝚺𝑽 T ,where 𝑼 is an 𝑀 × 𝑀 orthogonal matrix, 𝚺 is an 𝑀 × 𝑁 rectangulardiagonal matrix, and 𝑽 is an 𝑁 × 𝑁 orthogonal matrix. By using 𝑮 =𝑼 𝚺𝑽 T , 𝑪 and 𝑪 −1 can be written as
𝑪 = 𝑮T 𝑮 + 𝜆𝑰
= 𝑽 (𝚺T 𝚺 + 𝜆𝑰 )𝑽 T
and𝑪 −1 = 𝑽 (𝚺T 𝚺 + 𝜆𝑰 )−1 𝑽 T .
Gain vector 𝒙 can be written as
𝒙 = 𝑪 −1 𝒉T= 𝑽 (𝚺T 𝚺 + 𝜆𝑰 )−1 𝑽 T 𝒉T .
Consequently, we obtain𝒙T 𝒙 = 𝒉𝑽 (𝚺T 𝚺 + 𝜆𝑰 )−2 𝑽 T 𝒉T
(
)2
𝑁
∑
𝒉𝒗𝑛
=
,
𝜎𝑛2 + 𝜆𝑛=1where 𝒗𝑛 is the 𝑛-th column vector in 𝑽 , 𝒉𝒗𝑛 is the dot product of 𝒉and 𝒗𝑛 , 𝜎𝑛 is the 𝑛-th diagonal element in 𝚺, and 𝜎𝑛 = 0 if 𝑛 > 𝑀. 𝒗𝑛 and𝜎𝑛 represent the orthogonal basis and its scale in 𝑮. Because 𝜎𝑛 = 0 for
𝑛 > 𝑀, 𝒗𝑀+1 , ⋯ , 𝒗𝑁 represent the null space of 𝑮.
The above equation suggests that 𝒙T 𝒙 increases in the following twocases:••𝒉 is large; that is, a sensor position has high sensitivity to the ROIs.
𝒉𝒗𝑛 is large when 𝜎𝑛 = 0; that is, the sensitivity vector of a sensor
position is in the null space of the leadﬁeld matrix of the existing
sensors.

Supplementary material
Supplementary material associated with this article can be found, inthe online version, at 10.1016/j.neuroimage.2023.120257
References
Aine, C., Huang, M., Stephen, J., Christner, R., 2000. Multistart algorithms for MEG empirical data analysis reliably characterize locations and time courses of multiple sources.
Neuroimage 12, 159–172. doi:10.1006/nimg.2000.0616.
11

Y. Takeda, T. Gomi, R. Umebayashi et al.

NeuroImage 277 (2023) 120257

Sato, M., Yoshioka, T., Kajihara, S., Toyama, K., Goda, N., Doya, K., Kawato, M., 2004.
Hierarchical bayesian estimation for MEG inverse problem. Neuroimage 23, 806–826.
doi:10.1016/j.neuroimage.2004.06.037.
Sekihara, K., Sahani, M., Nagarajan, S.S., 2005. Localization bias and spatial resolution ofadaptive and non-adaptive spatial ﬁlters for MEG source reconstruction. Neuroimage
25, 1056–1067. doi:10.1016/j.neuroimage.2004.11.051.
Seymour, R.A., Alexander, N., Mellor, S., O’Neill, G.C., Tierney, T.M., Barnes, G.R.,
Maguire, E.A., 2021. Using OPMs to measure neural activity in standing, mobile participants. Neuroimage 244, 118604. doi:10.1016/j.neuroimage.2021.118604.
Shah, V.K., Wakai, R.T., 2013. A compact, high performance atomic magnetometer for biomedical applications. Phys. Med. Biol. 58, 8153–8161.
doi:10.1088/0031-9155/58/22/8153.
Storey, J.D., Tibshirani, R., 2003. Statistical signiﬁcance for genomewide studies. Proc.
Natl. Acad. Sci. U S A. 100, 9440–9445. doi:10.1073/pnas.1530509100.
Suzuki, K., Yamashita, O., 2021. MEG current source reconstruction using a meta-analysis
fMRI prior. Neuroimage 236, 118034. doi:10.1016/j.neuroimage.2021.118034.
Takeda, Y., Suzuki, K., Kawato, M., Yamashita, O., 2019. MEG source imaging and group
analysis using VBMEG. Front. Neurosci. 13, 241. doi:10.3389/fnins.2019.00241.
Tierney, T.M., Holmes, N., Mellor, S., López, J.D., Roberts, G., Hill, R.M., Boto, E.,
Leggett, J., Shah, V., Brookes, M.J., Bowtell, R., Barnes, G.R., 2019. Optically pumped

magnetometers: from quantum origins to multi-channel magnetoencephalography.
Neuroimage 199, 598–608. doi:10.1016/j.neuroimage.2019.05.063.
Van Veen, B.D., van Drongelen, W., Yuchtman, M., Suzuki, A., 1997. Localization of brain
electrical activity via linearly constrained minimum variance spatial ﬁltering. IEEE
Trans. Biomed. Eng. 44, 867–880. doi:10.1109/10.623056.
Vivekananda, U., Mellor, S., Tierney, T.M., Holmes, N., Boto, E., Leggett, J., Roberts, G.,
Hill, R.M., Litvak, V., Brookes, M.J., Bowtell, R., Barnes, G.R., Walker, M.C., 2020.
Optically pumped magnetoencephalography in epilepsy. Ann. Clin. Transl. Neurol. 7,
397–401. doi:10.1002/acn3.50995.
Wakeman, D.G., Henson, R.N., 2015. A multi-subject, multi-modal human neuroimaging
dataset. Sci. Data 2, 150001. doi:10.1038/sdata.2015.1.
Wittevrongel, B., Holmes, N., Boto, E., Hill, R., Rea, M., Libert, A., Khachatryan, E.,
Van Hulle, M.M., Bowtell, R., Brookes, M.J., 2021. Practical real-time MEG-basedneural interfacing with optically pumped magnetometers. BMC Biol. 19, 158.
doi:10.1186/s12915-021-01073-6.
Yoshioka, T., Toyama, K., Kawato, M., Yamashita, O., Nishina, S., Yamagishi, N., Sato, M.,
2008. Evaluation of hierarchical bayesian method through retinotopic brain activities reconstruction from fMRI and MEG signals. Neuroimage 42, 1397–1413.
doi:10.1016/j.neuroimage.2008.06.013.

12