Imaging somatosensory cortex responses measured by OPM-MEG: Variational free energy-based spatial smoothing estimation approach.

In recent years, optically pumped magnetometer (OPM)-based magnetoencephalography (MEG) has shown potential for analyzing brain activity. It has a flexible sensor configuration and comparable sensitivity to conventional SQUID-MEG. We constructed a 32-channel OPM-MEG system and used it to measure cortical responses to median and ulnar nerve stimulations. Traditional magnetic source imaging methods tend to blur the spatial extent of sources. Accurate estimation of the spatial size of the source is important for studying the organization of brain somatotopy and for pre-surgical functional mapping. We proposed a new method called variational free energy-based spatial smoothing estimation (FESSE) to enhance the accuracy of mapping somatosensory cortex responses. A series of computer simulations based on the OPM-MEG showed better performance than the three types of competing methods under different levels of signal-to-noise ratios, source patch sizes, and co-registration errors. FESSE was then applied to the source imaging of the OPM-MEG experimental data.

Introduction

The topographic organization of the somatosensory cortex was first established by Penfield and Boldrey using an invasive technique of direct electrical stimulation of the human cortical surface (Penfield and Boldrey, 1937). Their study revealed that the primary somatosensory in the posterior central cortex is organized in a one-to-one representation with the contralateral body surface, which is known as the somatosensory “homunculus.” With the development of brain functional imaging techniques including magnetoencephalography (MEG), electroencephalography, and functional magnetic resonance imaging, somatosensory evoked responses recorded outside the brain enable noninvasive studies of the somatosensory cortex (Nakamura et al., 1998; Baumgartner et al., 1993; Willoughby et al., 2020). Among these techniques, MEG provides a direct measure of weak extracranial neuromagnetic fields produced by thousands of neurons that act synchronously (Hämäläinen et al., 1993). It provides sufficient resolution to localize the precise source location and capture the time-frequency dynamics of evoked responses. Somatosensory evoked fields (SEFs) have been used to localize the clinical-style primary somatosensory cortex (SI) in pre-surgical neuroimaging (Solomon et al., 2015), to study the reorganization of the somatosensory cortex (Papadelis et al., 2018), and assess plastic changes in various physiological and pathophysiological states (Badura-Brack et al., 2015).

Superconducting quantum interference device (SQUID)-based MEG has developed rapidly over the last four decades. However, conventional SQUID-MEG is thought to be cumbersome and has a fixed sensor array, which limits its further use, such as in studying the neurodevelopment of children in whom the head circumference rapidly changes with time. In addition, its sensors need to work at a low temperature, which prevents it from getting close to the subject’s scalp to obtain higher signal-to-noise ratio (SNR) data. Fortunately, new types of sensors, including optically pumped magnetometers (OPMs) (Osborne et al., 2018) and high-temperature SQUIDs (Faley et al., 2017) make on-scalp MEG possible (Schneiderman et al., 2019). In particular, OPM-MEG has attracted the attention of researchers owing to its wearable and flexible sensor configuration (Boto et al., 2018; Iivanainen et al., 2017). OPM-MEG is currently undergoing technical exploration. Researchers from various backgrounds have focused their energy on OPM-MEG system construction (Borna et al., 2020), active shielding (Iivanainen et al., 2019), helmet design (Hill et al., 2020), sensor array design (Beltrachini et al., 2021), co-registration with MRI (Zetter et al., 2019; Cao et al., 2021), and applications (Iivanainen et al., 2020; Boto et al., 2021; Wittevrongel et al., 2021). In this study, we constructed a 32-channel wearable OPM-MEG system and aimed to use it to analyze the SEFs of four subjects under median and ulnar nerve electrical stimulations.

For the source analysis, various magnetic source imaging (MSI) techniques have been developed, such as the classic and powerful methods, namely the minimum norm estimation (MNE) method (Liu et al., 2002) and its variant. However, these methods are believed to blur the source extent (Grech et al., 2008); that is, they overestimate the source extent. This blurred estimate is because source signals reconstructed at two different locations may be sensitive to brain activity in the same location; namely, the existence of cross talk and point spread between sources (Samuelsson et al., 2021; Liu et al., 2002). Although the sources can be separated from the background by subjectively adjusting the imaging threshold, the adjustment lacks objective evidence (Sohrabpour and He, 2021). To realize the full potential of OPM-MEG, it is necessary to improve the accuracy of the source extent estimation. In recent years, solving the MEG inverse problem under the Bayesian framework has attracted significant attention (Cai et al., 2019; Costa et al., 2017; Liu et al., 2020). The Bayesian approach converts the MSI problem to the determination of the source covariance and enables the comparison of different imaging results. In this study, under the Bayesian framework, we propose a new MSI algorithm called variational free energy-based spatial smoothing estimation (FESSE). Simulations at different levels of the SNR, source patch size, and co-registration error show that it has better performance in estimating the source location and its extent compared to the benchmark algorithms, including MNE, beamforming—LCMV (Van Veen et al., 1997) and the multiple sparse priors (MSP) algorithm (Friston et al., 2008a). FESSE is applied to the source analysis of the measured OPM-MEG data and it shows better imaging results than the other methods.

Results

FESSE evaluation via simulations

A series of computer simulations using the OPM-MEG system were conducted to evaluate the FESSE performance. FESSE was compared with benchmark algorithms, including beamforming—LCMV, MSP, and MNE. LCMV was provided by the software Fieldtrip (Oostenveld et al., 2011) while MSP and MNE were based on SPM software (Litvak et al., 2011). Figure 1 schematically depicts the simulation protocol.

Figure 1

Simulation process of the OPM-MEG data

The location of the source patch was randomly chosen in the designed source space (shown as the red mesh), and its orientation was taken perpendicular to the local mesh. The source was a bandpass signal. Using the OPM-MEG sensor configuration, the data were computed by solving a forward problem. Different levels of noise could be added to the sensor data.

The aim of the proposed FESSE method was to enhance the ability of MSI to map somatosensory responses. Thus, we first investigated the performance of FESSE and the benchmark methods when the source patch was located in Brodmann area 3b (BA3b-L) of the SI. Figure 2 shows the reconstruction results of the source patch with varying sizes, that is, FWHM = 5, 10, and 30 mm. It can be observed that FESSE has better fidelity (higher AUC; AUC, area under the curve) for imaging the focal or extended source compared with the benchmark algorithms. As mentioned previously, the MNE shows blurred reconstruction results. LCMV and MSP provide focal imaging results, with the localization results of LCMV being more accurate than those of MSP.

Figure 2

Example reconstruction results

The source patches with varying spatial sizes were located in Brodmann area BA3b-L. The first column shows the simulated source, while the other columns show the reconstruction results of benchmark algorithms and FESSE. The imaging results are normalized using the maximum value of the estimated source amplitude for each method and shown on the inflated cortex.

Influence of the source patch size

To quantitatively illustrate the influence of the source patch size, we simulated data using candidate sources with different FWHM values to measure the effect of the source patch size on reconstruction accuracy. The SNR was set at 10 dB. Figure 3 shows the AUC and dipole location error (DLE) values of the four approaches. For all the source configurations considered, FESSE could reconstruct the actual source patch more accurately than the other methods. In particular, when the size of the source patch increased, its ability to reconstruct the source spatial extent also improved accordingly. FESSE also performed well in terms of localization accuracy. In our simulations, LCMV showed better performance for imaging the focal source than for imaging the extended source, which was manifested by its decreasing AUC and increasing DLE. The AUC values of the MSP and MNE decreased as the size of the source patch increased.

Figure 3

Influence of the source patch size

(A) Boxplots of AUC values of all methods at different sizes of the source patch.

(B) Boxplots of DLE values at each size of the source patch.

Influence of the SNR

Different levels of Gaussian noise were added to each group of simulated data generated by each candidate source with FWHM = 20 mm. The Gaussian noise was independent across sensors. We measured the performance of the MSI methods mentioned above for a wide range of SNRs. The simulation results are shown in Figure 4. FESSE had the lowest DLE compared with other methods under all SNR conditions and it showed better performance in estimating the source extent (AUC close to one except for SNR = −20 dB). The performances of LCMV and MNE improved with increasing SNR, as indicated by the increasing AUC and decreasing DLE. In other words, these two methods are susceptible to low SNRs. The ability of the MSP to estimate the source extent was quite stable with increasing SNR and it was almost better than MNE and LCMV; however, its DLE was not.

Figure 4

Influence of the SNR

(A) Boxplots of AUC results for all methods at varying SNRs.

(B) Boxplots of DLE results for all methods at varying SNRs.

It can be seen that the DLE of MSP has no steady tendency with SNR and it is almost higher than that of the other methods. For MSP, it assumes possible patch locations as candidate source priors and iteratively searches for the simplest source distribution that explains the most data. This causes a problem that the localization error cannot be reduced if the source priors are not given precisely. The current version of the MSP in the SPM toolbox adopts a strategy to assume multiple groups of random source patches to invert the same data several times and choose the best results as the final source estimate. This could reduce the source localization error to some extent; however, it cannot resolve this problem completely. This might explain why although the spatial priors are given in the MSP, its localization error is not better than that of the conventional methods (MNE and LCMV).

Influence of co-registration errors

In practice, co-registration errors are inevitable. Here, we studied the effect of co-registration errors on source reconstruction. The influence was assessed by adding different levels of rotation and translation errors to the sensor array when solving the inverse problem. According to our previous study on the accuracy of co-registration for OPM-MEG (Cao et al., 2021), we designed four levels of rotation and translation errors. Figure 5 shows the simulation results under the co-registration error when the source patch is 10 mm and the SNR is 10 dB. The reconstruction results of all MSI methods worsened with increasing co-registration error, which was indicated by the gradually decreasing AUC values and increasing DLE values. Although LCMV showed similar localization accuracy (see DLE values) with FESSE, its ability to recover the source extent (see AUC values) was the worst when co-registration errors existed simultaneously. Considering these two complementary metrics comprehensively, that is, AUC and DLE, FESSE provided better results (higher AUC and lower DLE) among all MSI methods. In addition, it should be noted that the co-registration error has a significant impact on the source reconstruction, especially the rotation error. For a rigid helmet, the relative locations of the sensors are known accurately when designing the helmet, so the co-registration error manifests as a rotation and translation relative to the subject’s head. When the sensors are rotated at the same angle, those far away from the rotation axis will have a larger displacement and greater location errors. Hence, it can be inferred that larger the helmet, greater is the effect of the rotation error. Owing to the unequal location errors, it is not surprising that the rotation error has a great impact.

Figure 5

Influence of co-registration errors

(A) The performance metric results under different levels of rotation error of sensor locations and orientations. The results of Rot = 0° are used as a contrast group without co-registration error.

(B) The performance metric results under different levels of translation error of sensor locations and orientations. The results of Tran = 0 mm are used as a contrast group without co-registration error.

Results of OPM-MEG experimental data analysis

As shown in Figure 6A, the OPM-MEG system was composed of second generation OPM sensors (QZFM, QuSpin Inc., United States), a magnetically shielded room (MSR), a three-dimensional (3D)-printed helmet, a data acquisition system, and an electrical stimulator. The co-registration of OPM-MEG and MRI was performed for each subject, and the co-registration results for one subject are shown in Figure 6B.

Figure 6

Optically pumped magnetometer-magnetoencephalography (OPM-MEG) system

(A) The diagram of the OPM-MEG system.

(B) The co-registration results of one subject.

We performed a sensor-level analysis for all the subjects. First, we evaluated the SNR of each dataset using post-stimulus and baseline data. The SNRs are all greater than 5 dB, which indicates the high quality of the measured SEFs. Figure 7 presents representative results of the somatosensory responses of Subject 1. The results for the same type of nerve stimulated in the left and right hands showed analogous SEF waveforms. The differences between the median and ulnar nerve responses were in the value of the maximum amplitude and its occurrence time. The maximum amplitude of the median nerve was higher than that of the ulnar nerve, which occurred at approximately 35 ms (M35), while that for ulnar stimulation occurred at 60–70 ms (M60). The TFRs showed that the power of the activity for all stimulus types was concentrated in the 20–100 ms period and it had a broad frequency band ranging from 10 to 40 Hz at 20–30 ms. For all stimulation types, a dipolar pattern was visible in the topography at the time of maximum magnetic amplitude. In addition, the dipole patterns of the left and right stimuli showed opposite polarities.

Figure 7

Representative results of the sensor-level somatosensory responses of Subject 1

(A) Somatosensory evoked fields (SEFs) of Subject 1 under left median electrical stimulation. The left panel shows the averaged SEFs with optically pumped magnetometer-magnetoencephalography (OPM-MEG) data in black and global field power (GFP) data in green. The time-frequency representation of the channel with the maximum of the magnetic field and the topography of the component with the maximum amplitude are shown in the upper and lower portions of the right panel, respectively.

(B) Somatosensory evoked fields (SEFs) of Subject 1 under right median electrical stimulation.

(C) Somatosensory evoked fields (SEFs) of Subject 1 under left ulnar electrical stimulation.

(D) Somatosensory evoked fields (SEFs) of Subject 1 under right ulnar electrical stimulation.

The SEFs of other subjects provided similar results to those of Subject 1, as shown in Figures S1, S3, and S5. After a comprehensive consideration of the occurrence time of the peak of the GFP, as well as the existence of the dipolar pattern in the topography, we assessed the peak stages of two nerve responses for all subjects after each type of stimulus, with the results presented in Table 1. In the post-stimulus period of 100 ms, two peak stages were recognized: M35–40 and M60–70. In particular, the M35 component was present for each stimulus. Thus, M35 was regarded as the major activity in our study, which was mapped to the cortex using the MSI methods mentioned above.

Table 1

Peak stages of brain activity for each subject

	Left median	Right median	Left ulnar	Right ulnar
Subject 1	M20, M35, M65	M20, M35, M65	M35, M65	M35, M60
Subject 2	M20, M35, M65	M20, M35, M70	M35, M60	M35, M70
Subject 3	M20, M35, M60	M20, M35, M60	M35, M60	M35, M50
Subject 4	M20, M35, M60	M20, M35, M60	M35, M60	M35, M60

The time window is within 100 ms of the post-stimulus period.

For the source-level analysis, we mapped the source activity at 35-ms post-stimulus (M35) using FESSE and benchmark methods. Using the measured OPM-MEG data, we successfully mapped the source of the somatosensory response. Representative image results for Subject 1 are shown in Figure 8. For all methods, the estimated source of each stimulus was localized near the somatosensory cortex; that is, around the post-central gyrus and it was on the contralateral hemisphere of the stimulus side, which is consistent with previous research results based on SQUID-MEG (Del Gratta et al., 2002; Schulz et al., 2004). In our study, FESSE provided a clear and concentrated source. From the results of FESSE (Figure 8), it could be observed that the source locations were all located at the anterior wall of the posterior central gyrus for all stimulus types. The estimated source extent for each stimulus was similar, and the average FWHM of the estimated source was 20 mm. In contrast, the source location and the spatial extent under the left/right median nerve stimulus were close to the cortex surface, while those under ulnar nerve stimulation were close to the central sulcus. The source imaging results of the other subjects are shown in Figures S2, S4, and S6. These results indicate that the median and ulnar nerve evoked responses could be discriminated in the source-level analysis. It can also be observed that the locations of the reconstructed sources were slightly different for different subjects, indicating individual differences in the cortex.

Figure 8

Source reconstruction results of all stimulus types for Subject 1

Each column shows the reconstruction results of each method. The normalized estimated source power is shown on the cortex of Subject 1. The source imaging results of LCMV and MNE are shown at a threshold of 50% of the maximum power value.

Discussion

In this study, we constructed a 32-channel OPM-MEG system and used it to perform median and ulnar nerve stimulation experiments. To accurately map somatosensory cortex responses, we proposed a new MSI method called FESSE. We verified its better performance compared to the benchmark algorithm via simulations and then applied it to the source reconstruction of SEFs. The measured OPM-MEG data were analyzed at the sensor and source levels.

Based on the Bayesian framework, different types of MSI methods can be derived from different constraints of the source covariance matrix (Wipf and Nagarajan, 2009; López et al., 2014). MNE and LCMV assume that the source prior is the identity matrix and the inverse of the data covariance matrix, respectively (Liu et al., 2002; Van Veen et al., 1997), while MSP and FESSE provide a spatial prior constraint, which can better constrain the spatial expansibility of the sources. For MSP, it assumes a limited number of candidate source patches of the same size in the source space and attempts to use GS to search a subset of source patches highly correlated with the observation as the inverse solution (Friston et al., 2008a). However, the estimated results of the MSP are highly dependent on the selected source patches, and the source patches are randomly chosen in the entire source space increasing the randomness of the source localization error. FESSE uses the variational free energy as model evidence to choose the appropriate spatial smoother for all candidate sources within the limited source space estimated by MNE or LCMV. In this sense, FESSE incorporates the advantages of the MNE and MSP methods.

Simulation results under different levels of SNR (Figure 4) and source patch size (Figure 3) show the effectiveness and superiority of FESSE over the competing methods MNE, LCMV, and GS-based MSP. The co-registration of OPM-MEG is more complicated than that of SQUID-MEG on account of its flexible sensor arrays (Cao et al., 2021; Zetter et al., 2018). Moreover, its close distance to the scalp and an increased SNR makes OPM-MEG sensitive to the forward model error. Therefore, it is necessary to analyze the impact of co-registration errors on source reconstruction. In our work, we evaluated the performance of MSI methods under co-registration errors. FESSE provides better source reconstruction results than the others (Figure 5). It should also be noted that the impact of the co-registration error is significant, and it is strongly recommended to use devices and methods with high co-registration accuracy. For FESSE, the smoothness level is a parameter that must be set previously. To reduce the computational burden, it can be set at intervals to produce spatial priors from the focus source to the extended source. The size of the spatial prior under different smoothness levels can be estimated by calculating the FWHM of each spatial prior.

Using the constructed wearable OPM-MEG, we have shown that the SEFs under median and ulnar nerve stimulations can be reliably detected by OPMs with high quality (maximum amplitude reaches 1 pT). Both the sensor-level and source-level analyses have demonstrated substantial differences in cortical responses to different stimulations of peripheral nerves, which is consistent with the findings of previous studies (Theuvenet et al., 2006; Huttunen et al., 2006). In particular, FESSE provides more reasonable source imaging results than competing methods without adjusting the imaging threshold (Figure 8). It can distinguish the cortex responses under the two types of electrical sensory stimulation.

Based on the prior information provided by studies conducted previously on the somatosensory cortex (Nakamura et al., 1998; Baumgartner et al., 1993; Willoughby et al., 2020), we distributed the sensors in our experiments to primarily cover the partial brain functional area (bilateral parietal and temporal lobes of the brain) instead of the whole scalp. Simulation results on the influence of the sensor layouts on the FESSE performance have demonstrated that different sensor layouts affect the source estimate, and the source reconstruction using the local distribution sensor layout is better than using a global distribution layout for detecting sources in the somatosensory cortex (see Figures S7 and S8). Our experimental results also verified the ability of the locally distributed OPM array to successfully map the source activity of the somatosensory cortex. In addition, we also preliminarily simulated a decreasing number of sensors on the FESSE performance. Simulation results (see Figure S8) show that when the number of sensors is greater than 20, FESSE can provide a source localization error of less than 5 mm. Currently, OPM sensor technology is at the initial exploration stage; therefore, few sensors are available for researchers. With a limited number of sensors, covering the functional area being studied can allow for convenient research without loss of important information.

Limitations of the study

It should be noted that the source imaging results of FESSE are based on the limited source space obtained from the imaging results of, for example, the MNE method, which makes FESSE unable to reconstruct the source areas that are not detected by the MNE method. FESSE is suitable for estimating the source and its range in the region of interest of the brain.

In our work, we used the same sensor arrays across different subjects without considering individual differences in the human cortex. Relative positions of the sensor array and heads of the subjects were slightly different because of their personalized head structures. The measured signal for each subject might not have the best quality because the sensors could not accurately cover the brain functional area of each subject. In other words, the sensor layout for detecting the somatosensory responses for each subject was different. Our simulations have shown that different sensor layouts influence source reconstruction. One of the outstanding advantages of OPM-MEG is its flexible sensor configuration, which is not present in conventional SQUID-MEG. To further improve the detection accuracy, optimization of the sensor array with a limited number for a specific subject is required accordingly. Research has been performed to consider the sensor configuration design from the perspective of source reconstruction accuracy (Beltrachini et al., 2021; Duque-Muñoz et al., 2019). However, for practical application, designing and fabricating a personalized helmet for each subject is complex and time-consuming and many studies still use a general helmet for different subjects (Iivanainen et al., 2020). Therefore, creation of helmet designs (flexible or rigid) that would enable the sensor array to be adjusted freely, timely, and accurately for each subject would be beneficial in the future.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Xiaolin Ning (ningxiaolin@buaa.edu.cn).

Materials availability

This study did not generate new unique reagents.

Experimental model and subject details

Human subjects

Four right-handed healthysubjects (two males and two females, 25–27 years of age) participated in this study. This study was reviewed and approved by the Ethics Committee of Beihang University. All subjects provided written informed consent for the experimental procedure, in accordance with the Declaration of Helsinki.

Method details

Variational FESSE

Estimation of the source extent is of interest to researchers for mapping the somatosensory cortex; however, it is also a major challenge for MSI algorithms. As widely described in the literature (Sohrabpour and He, 2021; Becker et al., 2015), MNEs generally lead to blurred source localization results. In recent years, MSI algorithms based on the Bayesian framework have been developed for source reconstruction. Here, we propose an MSI algorithm based on the Bayesian approach for source reconstruction.

Description of the inverse problem

The sources of MEG can be equivalent to a large number of current dipoles distributed throughout the cortical surface. The relationship between the magnetic fields measured at a set of N sensors for N time points and N dipoles can be described by a linear model:where is the additive noise and is the lead field matrix, which can be obtained from the analytical or numerical solutions of the forward problem in the quasi-static approximation (Mosher et al., 1999; Nolte, 2003). When and are provided, the inverse problem is to estimate the source activity . In a Bayesian framework, this estimation is transformed into maximizing the posterior distribution of conditioned on the observed data , that is, , which can be computed via the Bayes rule:

With proper spatial whitening based on the estimated noise covariance matrix (Ramírez et al., 2011), the noise can be assumed to present a zero-mean Gaussian distribution, that is, . Subsequently, the likelihood associated with (2) can be written as . Given the dataset, p() is a constant value. Thus, the estimation of relies on the determination of the source prior p(). Generally, Gaussian process models are used to represent p() as a Gaussian process prior (Friston et al., 2008b; Harrison et al., 2007), . This allows us to express the posterior distribution aswhere is the transpose operator. In this step, the maximization of is transformed into Taking the derivative of and setting it to zero yields (Dale and Sereno, 1993):

which is the solution to the inverse problem under the Bayesian framework. This equation can also be derived from other methods, such as the Tikhonov regularization method (see Liu et al., 2002, Appendix).

The source prior is , where is the initial prior of and h is the regularization parameter. The essence of different inverse methods can be described by the different constraints on the prior covariance For example, if there is no available spatial prior about the sources, the MNE solution always assumes that (Hincapié et al., 2017), and subsequently,

However, MNEs cannot precisely estimate the extent of the source. To obtain an accurate estimation of the source extent, the brain source structure within the limited source range can be modeled into the source prior, which is what our method attempts to accomplish.

Graph laplacian spatial smoother

Based on the MNE solution, the active sources are limited to The purpose of this step is to accurately estimate the location and size of the source within a limited source range. We represent the region of the non-zero values of where N is the number of candidate sources within . The location of each candidate source is denoted as To estimate the size of the source, graph Laplacian spatial smoothers with different smoothness scalars are performed on the cortex surface mesh. The graph Laplacian spatial smoother is defined based on an adjacency matrix , which is calculated using the vertices and faces defined in the original cortex mesh. If the vertices i and j are connected on a face element, the element in the i-th row and j-th column of matrix is 1, and 0 otherwise. The graph Laplacian spatial smoother is expressed as

where , and is the calculation of the sum of each row of . , denoting the diagonal matrix composed of a vector . is the factor controlling the smoothness level. The corresponding column in defines its smoothness for each candidate source in . By changing the smoothness , multiscale spatial priors are obtained. The larger the , the more extensive is the imaging region of the spatial prior.

Variational free energy as model evidence

For each spatial prior, the source covariance prior can be modeled as:where is the i-th column of , which represents the spatial prior; and is its hyperparameter. The change of the hyperparameter in (5) to guarantees the hyperparameter positive and applies a Gaussian assumption to it. Considering the sensor noise covariance, the source prior can be projected onto the sensor space.where denotes the hyperparameter of the sensor-noise covariance. Thereafter, the hyperparameters with the Gaussian hyperpriors need to be estimated, where and are the mean and covariance, respectively.

The advantage of the Bayesian approach for solving the inverse problem is that it enables the comparison of different models using the model evidence ; that is, in our method, and estimates hyperparameters from the data. However, it is difficult to directly calculate the model evidence. Fortunately, the variational free energy F provides a low bound to the log evidence (Friston et al., 2007). In other words, the maximization of free energy can approximate the model evidence. Considering the posterior densities of parameters presenting a Gaussian distribution: with the mean and covariance , the free energy can be expressed as (Friston et al., 2008a)where is the data covariance. We use to denote the determinant operator and to denote the trace of a matrix. Thefree energy is used as the maximum likelihood function to estimate the hyperparameters . With the expectation maximization (EM) algorithm, the hyperparameters can be estimated iteratively (Friston et al., 2008b).

For different spatial priors, it is natural to consider using free energy to provide evidence for model comparison. The spatial prior corresponding to the maximum free energy is the final source prior, which provides the imaging region of the estimated source. The source prior at the sensor space can be expressed aswhere and are estimated by the EM algorithm, and subsequently is taken into (4) to obtain the estimated source activity .

Evaluating FESSE performance by simulations

Simulation protocol

The cortex derived from MRI (MRI) data of Subject 1 (one of the subjects in our OPM-MEG experiment) was down-sampled to a mesh with 15,002 vertices to be regarded as the source model. The 32-channel OPM sensor configuration covering the bilateral parietal and temporal lobes of the brain of Subjectone was obtained accordingly. Radial field components were measured using an OPM array. In the simulations, the length of the data was 700 ms, starting from −200 ms to 500 ms, with a sampling frequency of 1 kHz. The source was activated only within 0–500 ms. We used full width at half maxima (FWHM) to define the size of each source patch to simulate sources with different active spatial extents and adopted the Gaussian spatial model to generate each source patch (Friston et al., 2008a). Each source patch was generated by selecting dipoles whose geodesic distances (Dijkstra, 1959) from the target dipole location were less than a given FWHM. The amplitudes of the dipoles within the source patch were set as where g is the geodesic distance and A is the intensity of the current source. A was set to 10 nA∙m in the simulations. The candidate source space was the cortex surface covered by OPM sensors. The orientation of the sources was perpendicular to the local cortical surface (Wipf et al., 2010). The source signals were bandpass orthogonal Gaussian signals (within 10–30 Hz). The corrected-sphere model (Nolte, 2003) was used as the forward model and used for data generation.

We evaluated the performance of all the source imaging methods to be compared under simulation situations of varied source patch size, SNR, and even co-registration errors. Under each simulation condition, we randomly selected candidate source locations on the designed source space to perform Monte Carlo simulations and only one source patch was activated in each dataset. The source patch size was chosen as FWHM = 10, 15, 20, 25, and 30 mm. The SNR was defined as , where is the root-mean-square of the matrix. We calculated the SNRs of the measured OPM-MEG data, which were within the range of 5–15 dB. To better explore the performance of the MSI methods, we set the SNRs in the simulations to a wider range, SNR = −20 dB, −10, 0, 10, and 20 dB. For our OPM-MEG system, we utilized a rigid helmet whose sensor locations and orientations were relatively fixed. In this case, the co-registration errors between the OPM sensor array and MRI were systematic errors, that is, the transformation errors of the estimated sensor array coordinate system and its real coordinate system (Hill et al., 2020). The transformation errors can be represented by the rotation and translation errors of the two coordinate systems. Thus, when solving the inverse problem, we added rotation and translation errors separately on the sensor array to evaluate the effects of the co-registration error on the performance of each MSI method. We chose the direction from the origin to the top of the head, that is, the z axis, as the rotation axis and the translation direction where the origin was defined as the midpoint of the left and right pre-auricular points. To analyze the effect of the rotation error, the whole sensor array was rotated with Rot = 0.3°, 0.6°, 0.9°, and 1.2°, about the z axis. Similarly, the whole sensor array was translated as Tran = 1, 2, 3, and 4 mm, along the z axis to evaluate the influence of the translation error.

For LCMV, we employed its current version in the Fieldtrip toolbox. The whole data, including the pre- and post-stimulus data, were used to estimate the common spatial filters, and the regularization parameter was set as lambda = 5% (the default setting of the ft_inverse_lcmv Fieldtrip function) to replace the data covariance matrix as its regularized version to guarantee a good estimate of the when the data was a low-rank matrix. Then the common filters were applied to the pre- and post-stimulus data separately. To avoid the centralization of source estimation, the estimated results were the estimated power of the post-stimulus data normalized by that of the pre-stimulus data, which is the so-called neural activity index. The implementation of MNE and MSP used the standard process provided by the SPM software. Their regularization parameters were set as Qe0 = 0.01 (the default setting of the spm_eeg_invert_classic SPM function) to bring the noise covariance matrix to a reasonable scale . For the MSP, a greedy search (GS) strategy was adopted. For FESSE, eight smooth levels of source patch were determined on the cortex of each subject, and the size of the source prior varied from 5 mm to 40 mm with a 5 mm increment.

Quantitative metrics

Two complementary metrics are used to quantitatively measure and compare the performance of the FESSE and benchmark algorithms. One is the dipole location error (DLE), which calculates the Euclidean distance between the global maximum of the estimated source activity and the simulated source. The DLE only uses the global maximum of the reconstructed activity to evaluate the accuracy of source localization; however, the local estimated activity and the corresponding spatial extent of the sources should also be considered. The other index, the improved area under the curve (AUC) (Grova et al., 2006) provides an appropriate measure of detection accuracy for source extent estimation. The curve is the receiver operating characteristic curve (Metz, 1986), which can be built based on the binary classification results of the normalized energy of the reconstructed and simulated activity under different choices of the imaging threshold (see Grova et al., 2006 for detail). The advantage of AUC is that it incorporates the imaging threshold into the statistical index and considers the source spatial change under different threshold selections. With these two metrics, the detection accuracy of MSI can be measured, such that the smaller the DLE and the closer the AUC is to 1, the better the effect of the source imaging.

Experimental data acquisition

Experimental paradigm

The subjects were presented with median and ulnar nerve stimulation at the right or left wrist with a bipolar electrode. The electrical stimulation was a square-wave electrical pulse with a duration of 200 μs, and was provided by a commercial stimulator (DS7A, Digitimer Inc., United Kingdom). Stimulus intensity was tailored individually and was set as the threshold value for each subject when a small movement of the thumb or little finger was visible (Antonakakis et al., 2019; Del Gratta et al., 2002). The inter-trial interval was about 1.3 s with a ±0.1 s time jitter. The subjects were instructed not to move their head and fixate on a point in front of them. A total of 300 trials were conducted for each stimulus. In total, four events, namely the right/left median and right/left ulnar nerve stimuli, were recorded for each subject, and the synchronizing trigger for each stimulus was also acquired.

OPM-MEG data acquisition

OPM-MEG system

During the experiment, the OPM sensors were operated inside the MSR (ambient-field amplitude and drift typically below 13 nT and 40 pT/h), and the subject was seated throughout the recording. A plastic board placed behind the seat was used to support the weight of the cables, which connected theOPM sensors to their electronic modules. The data acquisition system was composed of two pieces of the NI DAQ device PXI-4499, which enabled 32 channel data collection. Because of the need to use one data acquisition channel to record the synchronizing trigger of each stimulus, 31-channel OPM sensors (QZFM, QuSpin Inc.,United States) were used in the experiments. An array of 31 second generation OPM sensors was inserted into the sensor slots on a 3D-printed helmet and distributed symmetrically on both sides of the helmet to ensure coverage of the somatosensory functional area. The OPM sensors measured the radial components of the magnetic field. For different subjects, the same helmet was used, and the sensors were inserted into the same slots. The depths of the sensors for each subject were adjusted manually to ensure that the OPM sensors were as close to the scalp as possible. The electrical modules of the sensors and commercial stimulator were placed outside the MSR to avoid electromagnetic interference. All OPM-MEG data were acquired at a sampling frequency of 1 kHz.

Co-registration

After data acquisition for each subject, the insertion depth of each sensor was measured manually. Then the OPMs were removed from the helmet. The subject wearing the helmet was 3D-scanned immediately to obtain 3D superfine stereo images for the co-registration of OPM-MEG and MRI data. According to our previous studies on the co-registration comparison of three commonly used devices, we chose to use a laser scanner (HSCAN Prince 775, Scantech Inc., China), which has the best accuracy. Co-registration was accomplished in two steps: helmet match and face match (see Cao et al., 2021 for details of the co-registration method). Co-registration provided the sensor configuration for source imaging. The sensor locations and orientations were corrected according to the recorded depth. Through the co-registration step, the locations and orientations of the sensors relative to the head of each subject were obtained.

Anatomical MRI

MRI data of all subjects were collected using a Siemens MAGNETOM Prisma 3T MR system. We acquired T1-weighted MRI scans using an MPRAGE sequence (TR, 2,300 ms; TE, 3.03 ms; TI, 1,100 ms; FA, 8°; field of view, 256 × 256 × 192 mm; voxel size = 1 × 1 × 1 mm3). MRI data were pre-processed and segmented using Freesurfer software (Fischl, 2012) to obtain the scalp surface and cortex.

OPM-MEG experimental data analysis

The measured data were bandpass filtered between 2 and 40 Hz. We manually screened to identify bad channels and noisy segments. The data were then segmented into epochs based on the recorded triggers. For each group of data, the statistical metrics of each trial, including the variance, minimum, maximum, range, kurtosis, and z-value were computed. All trials were visually inspected and those with outliers of the statistical metrics were removed. With this step, the apparent artifacts caused by eye blinks and eye movements were rejected. After rejecting the noisy trials, for each stimulation type (left/right median nerve stimulation and left/right ulnar nerve stimulation), there were 214, 179, 201, and 79 trials left for Subject 1; 162, 178, 140, and 176 trials left for Subject 2; 162, 149, 108, and 133 trials left for Subjectthree; and 193, 173, 174, and 165 trials left for Subject 4. Compared to gradiometer-based SQUID-MEG, OPM-MEG is more susceptible to background magnetic field disturbances. Thus, we used the magnetic field mapping method (Mellor et al., 2021) to further correct the low-frequency artifacts caused by the natural head movement of the subject.

The pre-processed data were averaged over the time domain to obtain the averaged SEFs. The global field power (GFP), defined as the Frobenius norm of the averaged sensor data at each time point, was computed accordingly. In addition, the time-frequency representations (TFRs) were also calculated using the multi-taper-method convolution and the Hanning taper with a frequency-dependent window length. The time-frequency analysis was performed within the time interval from −200 to 500 ms and for a frequency from 2 to 40 Hz for each epoch, whichwere then averaged across epochs. An analysis of GFP and TFR provided significant information for estimating the main active period of the source activity.

The brain sources were modeled as distributed current dipoles on the down-sampled cortical sheet, with their locations on the vertices and orientations perpendicular to the local cortical surface (Wipf et al., 2010). The co-registration for each subject was performed to obtain the sensor configuration using the scanned 3D object and segmented scalp mesh from the MRI. The averaged SEFs within the estimated active period of the source were used to map the somatosensory responses to the median and ulnar stimulations. The lead field was computed using a corrected-sphere head model. Thereafter, FESSE and benchmark MSI methods were used for source imaging. The pipeline for each method was the same as that used in the simulations.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Software and algorithms
Fieldtrip v20201023	Oostenveld et al. (2011)	https://www.fieldtriptoolbox.org/download/
SPM12	The Wellcome Center for Human Neuroimaging	https://www.fil.ion.ucl.ac.uk/spm/
FreeSurfer v7.1.0	Fischl (2012)	https://surfer.nmr.mgh.harvard.edu/
MATLAB vR2019b	MathWorks	https://www.mathworks.com/products/matlab.html
FESSE algorithm	This paper	https://doi.org/10.5281/zenodo.5804077
Other
Human subjects	This paper	N/A