M/EEG source localization for both subcortical and cortical sources using a convolutional neural network with a realistic head conductivity model

Electroencephalography (EEG) and magnetoencephalography (MEG) are noninvasive brain imaging techniques that measure the electrical and magnetic fields generated by neuronal electrical activity. In contrast, other noninvasive measurements, like near-infrared spectroscopy (NIRS) and functional magnetic resonance imaging (fMRI), provide a low temporal resolution on the order of seconds due to their reliance on hemodynamic responses.¹ M/EEG offers high temporal resolution on the millisecond order, directly reflecting the neuronal activity.² Nevertheless, M/EEG's spatial resolution is limited by the volume conduction effect.³ Enhancement of the spatial resolution of M/EEG will contribute to neurological applications, from diagnostics to understanding various brain functions.

Many initiatives have been undertaken to enhance the spatial resolution of M/EEG using electrophysiological source imaging techniques (ESI).⁴ ESI represents an inverse problem that is inherently ill-posed, as there is no unique solution: multiple neural source combinations can produce identical M/EEG signal distributions.⁵ Traditional approaches require a priori assumptions, such as energy minimization or spatial smoothness, to deduce a unique solution.^6,7 However, the complexity of brain activity makes setting optimal regularization priors a significant challenge.

Recent advances in machine learning-based ESI techniques have sparked interest due to their data-driven approach that circumvents the need for explicitly defining a regularization term.^8–10 These methods have demonstrated superior localization accuracy of estimated sources in cortical regions, compared to traditional ESI methods.^8–10 All the previous studies focused on estimating cortical activity, so the effectiveness of machine learning-based ESI in localizing subcortical sources still needs to be confirmed.

Subcortical structures are pivotal to numerous brain functions, as highlighted in recent literature.^11,12 Their dysfunction is implicated in various diseases, such as Parkinson's disease¹³ and epilepsy.¹⁴ Characterizing the electrical activity in subcortical regions is essential yet challenging due to the limitations of noninvasive techniques and the restrictions on invasive recordings to clinically relevant areas. Contrary to the long-held belief that deep brain structures are not discernible via M/EEG, recent studies have shown the potential of these techniques in detecting subcortical activity through simultaneous M/EEG and invasive electrode recordings.^15,16 These findings, along with the remarkable success of machine learning-based ESI methods in imaging cortical sources,^8–10 suggest that such methods might also be adept at accurately imaging subcortical activities.

In this study, we explored the accuracy of ESI using a deep learning approach for not only cortical but also subcortical sources (Fig. 1). Based on the usefulness of convolutional neural networks (CNNs) in various types of inverse problems such as x-ray image reconstruction, sub-surface shear wave velocity estimation, and ESI,^8,17,18 we developed a four-layer convolutional neural network (4LCNN) designed to estimate source activity from M/EEG topography, which was trained using a vast amount of training data from forward M/EEG simulations based on a realistic head volume conductor model. The efficacy of our 4LCNN model for ESI was tested using a validation dataset created by forward M/EEG simulations and two types of human experimental data: (1) somatosensory evoked potentials (SEP) recorded by EEG and (2) concurrent recordings from invasive brain electrodes and MEG signals.

The architecture of our 4LCNN model (Fig. 2) is based on ConvDip, a previously published CNN model for cortical source estimation.⁸ Our study aimed to estimate not only cortical but also subcortical source activity by enhancing the original CNN model to increase its expressive ability. We increased the resolution of the input M/EEG map, the number of kernels (filters), and expanded the convolutional layers beyond those in ConvDip. Additionally, we upgraded the EEG/MEG recordings and individual head conductivity modeling compared to the previous study⁸ to capture the subtle differences between deep brain region sources. We used a higher number of EEG/MEG electrodes (160 channels) compared to the previous study (64 channels) and employed a more sophisticated realistic head conductivity model using the state-of-the-art head tissue segmentation model.¹⁹ Then, we compared the performance of the 4LCNN model with ConvDip and two commonly utilized ESI methods—eLORETA²⁰ and the LCMV beamformer²¹—using the same set of simulations and experimental data.

In the simulation, while we tested various noise levels in the EEG signals [with signal-to-noise ratios (SNRs) ranging from −20 to 30 dB], our primary focus was on EEG signals with SNRs between 5 and 30 dB, based on values reported in EEG studies on human subjects in evoked potential protocols.^22–24  Figure 3 shows examples of estimated cortical source activity by 4LCNN, ConvDip, LCMV, and eLORETA. The simulated sources were located in the medial part of the right superior frontal gyrus [Fig. 3(a)] and the left middle temporal cortex [Fig. 3(b)]. Visual inspection suggests that both 4LCNN and ConvDip accurately estimated the location and spatial spread. Conversely, the LCMV generally accurately estimated the location of activity but exhibited large spatial spreads. Variations in noise level did not seem to cause significant changes. eLORETA's estimated activities were dispersed at an SNR of 10 dB, but at higher SNR conditions, the estimations appeared more focal and accurate. However, eLORETA and LCMV failed to differentiate the left and right sides of the medial part of the superior frontal gyrus separately [Fig. 3(a)].

Figure 4 and supplementary material Fig. 1 present examples of estimated subcortical source activity. The simulated sources were located in the right thalamus [Fig. 4(a)] and the left pallidum [supplementary material Fig. 1(a)]. The 4LCNN method accurately estimated the location and spatial spread of these examples [Fig. 4(b) and supplementary material Fig. 1(b)], with no activity observed in the cortical regions. This accurate estimation was consistent across SNR conditions, with the exception of the left pallidum under the 10 dB SNR condition [supplementary material Fig. 1(b)]. ConvDip estimated activities near the subcortical structures without cortical activity for both the right thalamus and the left pallidum [Fig. 4(c) and supplementary material Fig. 1(c)]. However, the maximum activity was erroneously estimated in adjacent structures that were not the targets. LCMV estimated maximum activities in the cortical regions for both the right thalamus [Fig. 4(d)] and left pallidum data [supplementary material Fig. 1(d)]. For the right thalamus simulation [Fig. 4(d)], while broad activities were estimated within deep brain regions, the maximum activity was incorrectly localized in the medial orbitofrontal or insular cortices. In the pallidum simulation [supplementary material Fig. 1(d)], the maximum activity was wrongly estimated in the left medial temporal lobe, with high activity scattered across several subcortical regions. Similar to the cortical source simulation examples (Fig. 3), the SNR condition did not significantly impact LCMV estimation. For eLORETA [Fig. 4(e) and supplementary material Fig. 1(e)], in lower SNR conditions, estimated activities were widely scattered across cortical and subcortical regions, with maximum values mistakenly placed in the cortical regions for the examples of the right thalamus and the left pallidum. At a high SNR condition (SNR = 30), eLORETA accurately located the maximum activity for both examples, but the estimated activities were broadly distributed over several subcortical and medial cortical regions [Fig. 4(e) and supplementary material Fig. 1(e)].

We calculated three metrics (localization error, spatial dispersion, and area under the precision-recall curve (AUPRC)) to quantitatively evaluate the ESI performance of the 4LCNN compared to the other methods (Fig. 5). The localization error (distance between the locations of the sources with the maximum value of the simulated and estimated data) of the 4LCNN was significantly lower than that of the other ESI methods across all SNR conditions for both cortical and subcortical sources [Fig. 5(a)] (p < 0.05, see the supplementary material Excel data file for detailed statistical values). For cortical sources, the mean error was below 10 mm in the 4LCNN at SNRs higher than 10 dB, reaching up to 6.6 mm at an SNR of 30 dB. ConvDip showed relatively low localization errors; however, eLORETA was more accurate than ConvDip in high SNR conditions. Despite this, the mean error for eLORETA was 9.3 mm at an SNR of 30 dB.

For subcortical sources [Fig. 5(a)], 4LCNN also showed superior performance regarding localization error compared to other methods (p < 0.05, see the supplementary material Excel data file for detailed statistical values) except at an SNR of 5 dB, where no significant difference from LCMV was found. The mean error in the 4LCNN was below 10 mm at SNRs higher than 15 dB and decreased to 5.9 mm at an SNR of 30 dB. ConvDip generally had lower localization errors than LCMV and eLORETA; however, eLORETA was more accurate than ConvDip at an SNR of 30 dB. On the other hand, at SNRs below 25 dB, the error for eLORETA was larger than LCMV and ConvDip.

Spatial dispersion (spatial extent of sources) was consistently lower in the 4LCNN compared to the other ESI methods in all conditions for both cortical and subcortical sources [Fig. 5(b)] (p < 0.05, see the supplementary material Excel data file for detailed statistical values). The 4LCNN tended to estimate source activity more focally than the other two ESI methods.

The AUPRC, a comprehensive metric assessing the central location and spatial extent of estimated sources, was significantly higher in the 4LCNN than the other ESI methods across all SNR conditions for both cortical and subcortical sources [Fig. 5(c)] (p < 0.05, see the supplementary material Excel data file for detailed statistical values).

These results (Fig. 5) were derived from 300 simulations based on head models from three participants for validation (100 for each). Figure 2 of the supplementary material displays all three metrics for each participant's head model, showing consistent results with the pooled data in Fig. 5.

We then examined the relationship between the simulated and estimated activation area sizes (Fig. 6). As SNR increased, a stronger positive correlation between simulated and estimated area sizes was observed with the 4LCNN. For cortical sources [Fig. 6(a)], a moderate to high positive correlation was noted (r = 0.43–0.68) at SNRs of 15 or above. For subcortical sources [Fig. 6(b)], weak to moderate positive correlations were observed (r = 0.32–0.41) at SNRs of 15 or above. In contrast, negligible to very weak negative correlations were observed (r = −0.19 to 0.14) with the other ESI methods for both cortical and subcortical sources.

In addition to the data with SNRs between 5 and 30 dB, Fig. 3 of the supplementary material shows localization error, spatial dispersion, and AUPRC for data with lower SNR conditions ranging from −20 to 0 dB. Compared to data with SNRs higher than 5 dB, all metrics for the lower SNR data were largely worse. Even the 4LCNN, which demonstrated localization errors of less than 20 and 30 mm for cortical and subcortical regions, respectively, at 5 dB, showed errors of around 50 mm or more at 0 dB. Based on these results, these ESI algorithms should be used with data of higher SNR.

In addition to the simulation data, we further validated the 4LCNN model using real human data. For this purpose, somatosensory evoked potential (SEP) data acquired with EEG from three healthy participants (participant IDs: HY1–HY3) was utilized. It is hypothesized that somatosensory information from median nerve stimulation travels through the cervical spinal cord, medulla, thalamus, and somatosensory cortex at approximately 11, 13–14, 16–17, and 20 ms, respectively [Fig. 7(a)].²⁵ SEP data from all participants exhibited clear peaks corresponding to the SEP components [Figs. 7(a), 7(c), and supplementary material Fig. 3]. From the SEP components, we focused on the latter three as they originate from the brain. Each SEP component presented distinct topographic maps [Fig. 7(d)], and these maps were consistent across participants [Fig. 7(d) and supplementary material Fig. 4]. The SNRs of the SEP components were 18.6 ± 0.8 dB (mean ± SD), 20.1 ± 1.3 dB, and 22.0 ± 2.3 dB for the 13–14, 16–17, and 20 ms components, respectively.

We estimated source activity based on the topographic distribution of amplitudes from SEP components. The estimation results for participants HY1, HY2, and HY3 are shown in Fig. 8 and supplementary material Figs. 5 and 6, respectively. Using the 4LCNN method [Fig. 8(b) and supplementary material Figs. 5(b) and 6(b)], the 13 ms component's activity was accurately localized in the medulla for all participants, aligning with assumptions. For the 17 ms component, activity was localized in the left thalamus for two participants (HY1 and HY2). For the third participant (HY3), although the maximum activity was found slightly above the left thalamus in the left caudate, significant activity was also present in the left thalamus. The 20 ms component was consistently localized in the hand region of the left primary somatosensory cortex for all participants. In summary, the estimated source activity of SEP components was generally in accordance with a neural pathway of the somatosensory system.²⁶ 

ConvDip [Fig. 8(c) and supplementary material Figs. 5(c) and 6(c)] estimated the maximum activity for the 13 ms component in the mesial lobe (HY1) and the brainstem (HY2 and HY3). For the 17 ms component, the maximum activity was misplaced in the right parietal cortex (HY1), the left amygdala (HY2), and the left ventral diencephalon (HY3), rather than the left thalamus. Estimated activity for the 20 ms component was observed in the left primary sensory area (HY1–HY3) and additionally in the left primary motor area for two participants (HY1 and HY3).

LCMV [Fig. 8(d), supplementary material Figs. 5(d) and 6(d)] incorrectly localized the maximum activity for the 13 ms component in the mesial lobe for all three participants, with either large (HY1) or minimal (HY2 and HY3) medullary activity. All participants showed activity across extensive areas of the cerebral cortex for the 17 ms component, with pronounced activity in the left sensorimotor area, mesial lobe, and temporal lobe. The estimated activity for the 20 ms component was similar to that for the 17 ms component.eLORETA [Fig. 8(e) and supplementary material Figs. 5(e) and 6(e)] showed weak (HY2) or strong (HY1 and HY3) activity in the brainstem for the 13 ms component, yet only HY1 correctly localized the maximum activity in the brainstem (pons, not medulla). The other two (HY2 and HY3) exhibited maximum activity in the temporal (HY2) and occipital cortices (HY3). For the 17 ms component, despite observing weak activity in the left thalamus (HY3 and HY2), eLORETA incorrectly estimated large activity areas in the left sensorimotor area. For the 20 ms component, the primary activity was noted in the left primary sensory area; however, unlike with 4LCNN, there was no distinct peak activity localized in the hand sensory region.

We also analyzed the spatial dispersion of the estimated source activity from SEP components (Fig. 9). As with the simulation data [Fig. 5(b)], spatial dispersion for the 4LCNN was significantly lower than that observed in the other ESI methods (4LCNN vs ConvDip: p = 0.0083, effect size (ES) = 0.43; 4LCNN vs LCMV: p = 0.0032, ES = 1.97; 4LCNN vs eLORETA: p < 0.001, ES = 1.61; ConvDip vs LCMV: p < 0.001, ES = 2.01; ConvDip vs eLORETA: p = 0.020, ES = 1.46; LCMV vs eLORETA: p = 0.76, ES = 0.14).

To further evaluate the 4LCNN model on the real data, we conducted ESI using the 4LCNN on MEG data and evaluated its accuracy with simultaneously recorded invasive brain signal recordings from SEEG and ECoG (Fig. 10) from three individuals with epilepsy (subject IDs: EP1–EP3). The recordings were carried out in a resting state without any epileptic seizures, lasting for 5 min and repeated four to six times per participant. Since the 4LCNN is designed for analyzing single source topographic maps, it cannot directly process continuously recorded M/EEG signals, which are mixed signals from multiple sources. Resting-state brain activity consists of concurrent signals from multiple sources. Independent component analysis (ICA) has been used to separate these mixed signals into individual source activities.²⁷ Most of the separated sources exhibit dipolarity,²⁷ which is thought to be generated by a single source in a localized brain region. Therefore, the sources separated by ICA are considered applicable for the 4LCNN. Following a methodology from Hnazaee et al.,²⁸ which investigated the localization accuracy of EEG sources using SEEG recordings, we applied ICA to separate the MEG signals into distinct sources and evaluated the source localization accuracy using ESI methods. Blind source separation (BSS) techniques, including ICA, are utilized to extract single source activities of interest such as cognitive processes, resting brain components, and epileptic discharges.^29–31 Each independent component (IC) includes a time course and a spatial pattern (topography), and each topography can be applied to source localization.^30,31 Before this validation, we verified the accuracy of the ESI methods using MEG simulation data (refer to supplementary material Figs. 7 and 8), which produced results comparable to those obtained with EEG simulation data (Figs. 5 and 6).

Source activities of ICs separated from MEG signals were estimated. If the activity time courses of the ICs correlated with signals from invasive electrodes, the Euclidean distance between the location with the maximum activity in the estimated source activity and the most correlated electrode was calculated. Examples from a participant (ID: EP1) are provided in Fig. 11. Activation time series of two ICs were correlated electrodes located at the temporal lobe (i.e., an ECoG electrode) and the hippocampus (i.e., a SEEG electrode), respectively. The topography maps of two correlated ICs are shown in Figs. 11(a) and 11(b). Distances between the locations of maximum activity and the most correlated electrodes in the temporal lobe and hippocampus ranged from 7.3 to 38.5 mm for all ESI methods (Fig. 11). Notably, in this example, the 4LCNN method estimated a source location within the hippocampus that was very close to the corresponding electrode [7.3 mm, Fig. 11(b)].

From the three participants, 14 ICs with correlated invasive electrodes were identified. The SNR of the ICs was 10.1 ± 2.2 dB (mean ± SD). The distances between the maximum activity locations of the estimated source activities and the corresponding most correlated electrodes did not show significant differences among the four ESI methods (4LCNN vs ConvDip: p = 0.51, ES = 0.15; 4LCNN vs LCMV: p = 0.51, ES = 0.28; 4LCNN vs eLORETA: p = 0.51, ES = 0.48; ConvDip vs LCMV: p = 0.51, ES = 0.10; ConvDip vs eLORETA: p = 0.51, ES = 0.27; LCMV vs eLORETA: p = 0.51, ES = 0.21) [Fig. 12(a)]. Although the mean distances were not different among the methods, focusing on sources which were estimated spatially close to the correlated electrodes by the 4LCNN, sources with the correlated electrodes placed in the deep brain, cuneiform nucleus, and frontal pole were accurately estimated (distance ranged from 7.3 to 10.1 mm). These accurately estimated regions were relatively distant from the craniotomy site used for implanting the invasive electrodes. Spatial dispersion of the estimated source activity was also evaluated [Fig. 12(b)]. The spatial dispersion in 4LCNN and ConvDip was significantly lower than in LCMV and eLORETA (4LCNN vs ConvDip: p = 0.23, ES = 0.43; 4LCNN vs LCMV: p < 0.001, ES = 1.97; 4LCNN vs eLORETA: p < 0.001, ES = 1.61; ConvDip vs LCMV: p < 0.001, ES = 2.01; ConvDip vs eLORETA: p = 0.0030, ES = 1.46; LCMV vs eLORETA: p = 0.76, ES = 0.14).

This study is the first to show the feasibility of deep learning-based ESI for subcortical sources. The inverse problem of M/EEG source reconstruction, inherently highly ill-posed, has been tackled with our 4LCNN, which leverages realistic head conductive models for estimating cortical and subcortical activity from M/EEG data. Compared to ConvDip, eLORETA, and LCMV, 4LCNN outperformed them in various metrics through M/EEG simulation. Its effectiveness was also confirmed by its generalizability to experimental data, aligning with somatosensory pathways and invasive electrode recordings. Although its localization accuracy was comparable to other methods in datasets with concurrent invasive recordings, this might be due to MRI's limitations in accounting for ECoG-related brain deformations post-electrode implantation, potentially undervaluing 4LCNN's performance.³² Our findings highlight 4LCNN's potential as a highly accurate method for source localization for both cortical and subcortical sources, suggesting future contributions in various neurological applications, from diagnostics to understanding cognitive and motor functions.³³ 

Because of the high ill-posed nature of the inverse problem of M/EEG source reconstruction, most of the commonly used ESI methods require a priori assumptions to obtain a unique solution by explicitly setting regularization terms. Given the complex nature of brain activity, it is quite a challenge to set optimal regularization priors.³³ Even widely utilized LORETA variants, despite their effectiveness, can lead to diffused estimations of activity.^8,9,34 A previous study comparing different ESI methods with various regularization priors indicated that choosing between these methods involves a trade-off: smaller localization errors (e.g., sLORETA) vs reduced spatial dispersion (e.g., MNE).³⁵ Therefore, setting optimal regularization priors that adequately reflect the features of M/EEG sources is a complex task.

Recent approaches have deployed artificial neural networks (ANNs) to surmount this inverse problem, adopting a data-driven approach that dispenses with regularization priors.^8–10 For instance, Hecker et al.⁸ showed that a shallow CNN-based ESI method, ConvDip, could accurately estimate source localizations than conventional ESI methods. However, the success of ANN-based ESI methods has so far been demonstrated only for cortical sources.^8–10 Therefore, we explored the applicability of ANN-based ESI to both cortical and subcortical sources using the 4LCNN, a model inspired by ConvDip but enhanced to augment its expressive power.

In this study, 4LCNN showed the best performance across all the four ESI methods for SNRs between 5 and 30 dB, which are typical SNRs in event-related potential (ERP) and SEP studies.^22–24 Regarding the localization error for maximum value sources [Fig. 5(a)], mean errors were below 10 mm for SNRs higher than 10 dB for both cortical and subcortical sources. Remarkably, errors were as low as approximately 6.6 and 5.9 mm for cortical and subcortical sources, respectively, at an SNR of 30 dB. The model's robustness to noise and its capacity to focus estimations more precisely than other ESI methods underscore its potential. This accurate estimation may be attributed to the nonlinearity introduced by ReLU layers.³⁶ Additionally, the augmented number of learnable parameters, due to enhanced input 2D topographic resolution, filters, and layers in 4LCNN compared to the preceding CNN model (ConvDip⁸), likely contributed to the improved accuracy in source localization estimation. Notably, even though the 4LCNN was trained with 20 dB SNR data, it can accurately localize EEG data across a wide range of SNR levels. Such robustness to noise has also been reported in other ANN-based ESI methods.^8,10 The generalizability can be attributed to batch normalization³⁷ and dropout layers.³⁸ 

The 4LCNN estimates source activities more focally than other ESI methods [Fig. 5(b)]. For instance, in a representative example [Fig. 3(a)], 4LCNN could distinguish activities between the left and right regions within the medial part of the superior frontal gyrus—a task that poses challenges for conventional ESI methods like eLORETA and LCMV. Conversely, ConvDip and 4LCNN achieved more focal estimations of source activities. This increased focality, first noted in studies introducing ConvDip. Such high focality probably stems from high model's expressive power, which is linked to its nonlinearity and complexity. While both ConvDip and 4LCNN benefit from the nonlinearity provided by ReLU layers, the 4LCNN is a more complicated architecture. Thus, 4LCNN provides more focal estimations, even when trained on the same M/EEG simulation data as ConvDip. The estimates derived by 4LCNN were not just focal, but were able to follow source extent that was impossible with other ESI methods (Fig. 6 and supplementary material Fig. 7).

The effectiveness of neural networks in ESI tasks heavily relies on the traits of the training data, making the M/EEG simulation's modeling a critical step for the successful application of ANN-based ESI models to actual M/EEG data. To this end, we employed the state-of-the-art brain segmentation method, CHARM, which automatically generates a highly accurate segmented MRI with delineation of ten different tissue types.¹⁹ Additionally, we incorporated MR-compatible fiducial markers to precisely co-register MRI and EEG data. Such efforts for making precise head model prevents deviations between simulation-generated training data and actual measured M/EEG data, and is thought to contribute to the generalization performance of the created 4LCNN to real M/EEG data.

Validating 4LCNN with SEP data (Fig. 8 and supplementary material Figs. 5 and 6) revealed its ability to accurately localize sources in line with established neural pathways of the somatosensory system.²⁶ Invasive depth electrode recordings have provided evidence that somatosensory information relays the medulla and thalamus from 14–16 ms and 16–18 ms, respectively.²⁵ The 14 ms component showed localized activity within the lower brainstem (medulla), although lateralization remained undetected by the 4LCNN. The 17 ms component was localized in the left thalamus for two participants (HY1 and HY2), while the third participant (HY3) also exhibited significant activity in this region, albeit with the peak activity detected in the nearby left caudate. The somatosensory information by the median nerve stimulation reaches around 20 ms poststimulation.²⁵ The 4LCNN successfully estimated high activity in the somatosensory cortex's hand area³⁹ for all participants for the 20 ms component. Thus, the 4LCNN, although trained by simulated data, is capable of providing highly accurate estimations of source locations even when estimating activities in real EEG data.

When applying 4LCNN-based ESI to time-series activities in M/EEG data that involve complex overlapping from multiple sources, it is challenging to directly input the raw signals into the 4LCNN, as the model is trained to localize a single source. Hence, we employed ICA to separate the M/EEG data into distinct signal sources⁴⁰ and assessed the source localization accuracy with ESI methods. The mean location accuracy estimated by the 4LCNN was the most precise at 24.1 mm; however, it was not significantly different from the other ESI methods [Fig. 12(a)]. Possible reason for the lack of superiority of the 4LCNN in this validation is brain deformations resulting from the implantation of ECoG electrodes. Such implantations can lead to deformations due to electrode thickness, hematoma, and tissue swelling, causing unpredictable and non-uniform brain shifts of up to 24 mm in cortical areas and up to 3 mm in deeper structures.³² Following previous studies,^28,41 we created the head models for ESI from pre-implantation MRI data, because implanted electrodes distort MRI data. Thus, structural differences between the created head model and the head for which the signal was measured (an example shown in Fig. 13) potentially declined ESI accuracy in all the methods, probably masking the differences in estimation accuracy among the ESI methods. Notably, among the four sources that were precisely localized by the 4LCNN [Fig. 12(a)], they were located deeper or far from the craniotomy site, where less deformation is anticipated. Although this finding is tentative, it suggests the precise localization capabilities of the 4LCNN when an ideal head model is used. To further validate the 4LCNN, datasets involving simultaneous SEEG and MEG recordings without ECoG would be preferable, because SEEG typically does not result in substantial brain deformation.⁴² In previous research employing the LAURA method, an advanced form of MNE, for ESI from concurrent SEEG and EEG recordings, the source localization error ranged between 14.8 and 23.5 mm. In comparison, the accuracy of the 4LCNN for the aforementioned four sources, which are likely less impacted by brain deformation, was superior (accuracy ranging between 7.3 and 10.1 mm). Nevertheless, additional validation is necessary to further prove the precise localization capabilities of 4LCNN due to the limited number of sources and differing experimental conditions.

In addition to the brain deformation, another possible reason for the lower performance of ESIs in the validation in the MEG experiment is that ECoG and SEEG electrodes likely recorded activity from multiple nearby sources. In this validation, we compared the estimated source locations with the most correlated electrode positions. Considering that the electrodes captured activity from several sources, even if the electrodes showed some correlation with the estimated source activity (r > 0.2), they were close enough to record the source but likely not positioned at the center of the target source. As a result, the localization performance of the MEG data in all ESI methods was lower than in the simulation data. In a recent study, the activity of MEG ICs and SEEG electrodes was compared, as in the current study, but the target was interictal epileptiform discharges. That study demonstrated stronger IC-SEEG correlations, suggesting that IEDs could be considered a more reliable ground truth for validation of ESI methods. In addition to the IEDs, simultaneous intracerebral stimulation and M/EEG recordings are a suitable method for the ESI validation with reliable ground-truth. An open dataset that simultaneously measures intracranial stimulation and high-density EEG is available,⁴³ although the stimulation sites are limited to cortical regions. Utilizing this dataset will help further validate our ESI methods in future studies.

Our model focally estimated the brain activity compared to the other ESIs [Figs. 5(b), 9, and 12(b)]. The focality would be derived from the model being trained to localize a single source cluster. This model would be suitable for experimental paradigms where active sources are assumed to be temporally separated such as the short latency SEP experiment. However, brain activities in most cases are more complicated and several sources are active at the same time. For example, resting-state brain activity, which was used in the MEG experiment in the present study, involves multiple brain regions, and they are connected as functional network.⁴⁴ Given the strong coupling between the regions, it is possible that even separated components by ICA do not reflect activity from a single brain region, which has a gap from the 4LCNN model trained by single source cluster data. Thus, the validation presented here has limitations due to the nature of the recordings. To make a more reliable validation between the estimated activity and intracranial recordings, we should either use experimental paradigms that assume single-source activity, such as auditory brainstem response (ABR).⁴⁵ Another way to overcome the limitation is to extend our models for multiple source clusters. We can generate training data for the multiple source model by modification of the forward simulation in this study. The deep learning-based ESI model for multiple sources should be a more complex architecture than the 4LCNN, considering the complexity of the brain activity.

Another limitation of the 4LCNN is that it was trained to work on single time points of M/EEG data. As a result, while the model can predict the location of sources, it is unable to estimate the temporal dynamics of brain activity. Recent studies have succeeded in accurately estimating both spatial and temporal brain activity of the cerebral cortex using an approach that combines neural mass models (NMMs), which can simulate realistic spatiotemporal brain dynamics, with deep learning models.^9,46 By incorporating such NMMs to generate more realistic training data, the 4LCNN could achieve higher performance when applied to real human data. Since CNNs are effective at aggregating spatial information, our 4LCNN could be developed into a model capable of utilizing spatiotemporal information, such as CNN-LSTM models that combine CNN with long short-term memory (LSTM) layers commonly used for time series data.⁴⁷ By integrating LSTM and NMMs with our CNN-based approach, we may enable a more accurate estimation of spatiotemporal brain activity both in cortical and subcortical regions. The advantage of this new approach is that it allows temporal dynamics to be analyzed, but it is expected that the use of temporal information tied to a specific brain region will also increase the spatial accuracy of the new approach compared to the 4LCNN.

One disadvantage of using 4LCNN is its high computational cost. Compared to conventional ESI methods such as LCMV and eLORETA, deep learning-based ESIs require model training time. For instance, training 4LCNN took about four hours, while ConvDip took around 50 min on an Nvidia RTX A6000 GPU. Additionally, our approach requires the preparation of training data through simulations, which took approximately six minutes to generate 450 000 simulated datasets on an Intel i7-12700k desktop processor. As for the estimation times for 100 test data, while there were differences in computation times across the ESIs, they were all short enough not to hinder the execution of studies: 0.11 s for 4LCNN, 0.01 s for ConvDip, 0.7 s for LCMV, and 6.1 s for eLORETA. Reducing the training time for 4LCNN is an important aspect for future research. One potential strategy is to modify the model's output from electrical sources to brain parcels based on anatomical MRI analysis, which could significantly reduce the solution space.

In conclusion, this study is the first to show the feasibility of deep learning-based ESI for a single source in subcortical regions. We have proposed a data-driven ESI framework with the 4LCNN combined with a realistic head conductivity model. The results presented from a series of accuracy validations involving both numerical simulations and real human data demonstrate the high feasibility of the 4LCNN for accurate M/EEG source localization of both cortical and subcortical sources. The demonstrated potentials of the 4LCNNs in accurate source localization in this study indicates their future contributions to a range of research areas, including clinical diagnosis, understanding the pathophysiology of various neuronal diseases, and basic brain functions.

For EEG data analysis, we obtained MRI and EEG data from three healthy adult males aged 23–25 years (Participant IDs: HY1−HY3). For MEG data analysis, MRI, MEG, Electrocorticographic (ECoG), and stereo-electroencephalographic (SEEG) data were collected from three patients suffering from intractable epilepsy undergoing pre-surgical clinical assessment (participant IDs: EP1–EP3, aged 10–29 years). Placements of invasive electrodes for EP1–EP3 are displayed in supplementary material Figs. 9–11, respectively. Structural MRI and M/EEG sensor position data were used for M/EEG simulation. M/EEG signals were utilized to evaluate source estimation for implementations with real data. ECoG and SEEG data were employed for the localization accuracy evaluation of estimated sources from the MEG signals.

The locations of EEG electrodes were digitized using an optical digitizing system (Brainsight, Rogue Research, Montreal, Canada). To achieve precise co-registration of MRI and EEG electrodes, we attached three donut-shaped multimodal (MR/CT) fiducial markers (PINPOINT, Beekley Corp., Bristol, CT) with a hole size of 1.27 mm onto the left and right pre-auricular points (LPA and RPA) and the forehead, and digitized their locations. The coordinate systems of the electrode and MRI data were aligned based on these marker positions using the Fieldtrip software.⁴⁸ 

T1 (TR = 2400.0 ms, TE = 2.2 ms, 1 mm isotropic resolution) and T2-weighted (TR = 2500.0 ms, TE = 272 ms, 1 mm isotropic resolution) structural MRIs were obtained using a Siemens Prisma 3 T scanner (Siemens Medical Systems, Erlangen, Germany). As mentioned previously, three donut-shaped multimodal fiducial markers were attached to the head during MRI data acquisition for precise co-registration with EEG electrodes.

To simulate M/EEG data realistically and solve the inverse solution, we created a forward model based on anatomical head MRI images and M/EEG electrode positions.

We utilized SimNIBS equipped with the complete head anatomy reconstruction method (CHARM), which is a state-of-the-art technique, to automatically generate a highly accurate segmented MRI image with ten tissues: white matter, gray matter, cerebrospinal fluid (CSF), scalp, eyeballs, compact bone, spongy bone, blood, muscle, and air pockets, based on individual T1 and T2 MRI images.¹⁹ A hexahedral mesh was generated directly from the segmented MRI image using the “prepare_mesh_hexahedral” function within the Fieldtrip toolbox. To avoid staircase artifacts, we created geometry-adapted hexahedral meshes, where mesh nodes at tissue boundaries were slightly shifted to provide a smoother representation of the boundaries.⁴⁹ Then, utilizing finite element method (FEM) approaches, we constructed highly realistic multicompartment volume conductor models for forward EEG simulations. The conductivity (unit: S/m) of each tissue type was set as follows based on a review article:⁵⁰ white matter: 0.22, gray matter: 0.47, CSF: 1.71, scalp: 0.41, eyeballs: 0.5, compact bone: 0.0046, spongy bone: 0.050, blood: 0.57, muscle: 0.32. The conductivity for eyeballs was adopted from Haueisen et al.⁵¹ due to the absence of this information in McCann et al.⁵⁰ 

We employed T1 MRI data to obtain cortical surfaces and subcortical anatomical regions for defining the locations and orientations of the elementary dipole sources. The sources were placed on cortical surfaces and within subcortical volumes in accordance with a previous study.⁵² Initially, cortical surfaces were extracted using FreeSurfer.⁵³ Based on the cortical surface data, we created a high-resolution cortical surface triangular mesh with the HCP Workbench⁵⁴ and resampled it to approximately 8000 sources using the iso2mesh toolbox.⁵⁵ The edge lengths of the mesh were about 3–4 mm. Subcortical parcellations were obtained by FreeSurfer. The sources within subcortical volumes were evenly distributed with a grid resolution of 3 mm using the “ft_prepare_sourcemodel” function in Fieldtrip. For cortical sources, dipoles were placed with orientations normal to the cortical surface mesh.⁵⁶ For subcortical sources, triplets of orthogonal dipoles were situated within the subcortical volumes.¹⁵ The leadfield matrix, which comprises vectors of potential amplitudes received by electrodes from each dipole source, was generated using the DUNEuro toolbox in Brainstorm.^57,58

Furthermore, we developed a second forward model for validation of the 4LCNN. We modified the forward model used for training M/EEG data to create the second forward model for the validation, considering potential differences in brain-electric signal volume conduction between the actual head and our model, which cannot be precisely known a priori. We modified the electrode positions by adding Gaussian noise (mean: 0 mm, variance: 2 mm) to the three-dimensional coordinates of each electrode. Additionally, we altered the tissue conductivity in the FEM solution by scaling the original values by 1.2 or 0.8: white matter: 0.17 (original: 0.22), gray matter: 0.56 (original: 0.47), CSF: 1.36 (original: 1.71), scalp: 0.24 (original: 0.41), eyeballs: 0.6 (original: 0.5), compact bone: 0.0037 (original: 0.0046), spongy bone: 0.060 (original: 0.050), blood: 0.46 (original: 0.57), muscle: 0.38 (original: 0.32). All simulations in the validation section were performed using the validation model, whereas the calculation of the inverse solutions was based on the original model. The 4LCNN training data were simulated using only the original model. In this study, the original and modified models are referred to as the training model and validation models, respectively.

We generated two different sets of simulated EEG signals: (1) EEG data for training the 4LCNN model using the training model and (2) validation data to test the constructed 4LCNN model against the validation model.

For training the 4LCNN model, we simulated 100 000 and 350 000 samples for cortical and subcortical sources, respectively, using the training model. For model evaluation, 100 samples were simulated for each cortical and subcortical source estimation using the validation model.

In this study, we utilized the 4LCNN to estimate source activity from M/EEG signals by addressing the M/EEG inverse problem. The 4LCNN model was specifically designed to estimate a single clustered source from the spatial distribution map of M/EEG signals.

As a preprocessing step for EEG signals, we standardized the amplitude of each spatial distribution map, comprising M/EEG signals from electrodes at a given time instance, to have a zero mean and unit variance. These maps were then transformed into input data for a 2D CNN by interpolating onto a 2D image of size 26 × 26. From this 2D matrix, the 4LCNN model estimates the activity of distributed dipoles within the brain. The design and training of the 4LCNN model were performed using MATLAB's Deep Learning Toolbox on an NVIDIA RTX A6000 GPU.

The architecture of our 4LCNN model (Fig. 2) is based on ConvDip, a CNN model for cortical source estimation.⁸ Our study aimed to estimate both cortical and subcortical source activity by enhancing the original CNN model to increase its expressive ability. We augmented the resolution of the input M/EEG map and the number of kernels (filters) in the CNN and expanded the convolutional layers beyond those in ConvDip. The first layer in the 4LCNN was a convolutional layer with 32 filters, each of size 3 × 3, employing zero-padding. Each convolutional layer was followed by a batch normalization (BN) layer³⁷ to accelerate training and mitigate internal covariance shift problems. A rectified linear unit (ReLU) function followed each BN layer as the activation function. After four sets of convolution, BN, and ReLU layers, a dropout layer with a dropout rate of 30% was introduced to enhance the model's generalizability. The final two layers of the 4LCNN were a fully connected layer and an output layer.

The 4LCNN model parameters were optimized using adaptive momentum estimation (ADAM)⁶² with a learning rate = 0.001, β1 = 0.5, and β2 = 0.999. The loss function utilized was the root mean squared error (RMSE).

We assessed the quality of inverse source estimation using three performance metrics:

• Localization error: We calculated the Euclidean distance between the locations of the sources with the maximum value between the simulated and estimated data. This metric quantifies the accuracy of the estimated source center.⁸ 

• Spatial Dispersion: To assess the spatial extent of sources, we calculated the spatial dispersion metric as follows:³⁴ 

  $Spatial Dispersion=∑j=1majrpeak−rj2∑j=1maj2,$

  where m is the number of sources, a_j and r_j are the activity intensity and location of source j, respectively, $rpeak$ is the location of source, which has the maximum intensity across all sources, $rpeak−rj$ is the Euclid distance between $rpeak$ and $rj$⁠.

• Area under the precision-recall curve (AUPRC): We computed the AUPRC to assess the accuracy considering both location and spatial extent of sources. The AUPRC is more suitable than other metrics, such as the receiver operating characteristic (ROC) curve, for evaluating performance in imbalanced datasets like ours, which has a large imbalance between negatives and positives.⁶³ In our case, since simulated active sources are limited in a cluster, AU-PRC was considered to be a suitable method. The AUPRC quantifies the performance of a binary classification system as its discrimination threshold T is varied. The estimated sources are classified as active or inactive based on whether their amplitude exceeds the threshold T. We then tallied the true positives (TP), true negatives (TN), false positives (FP), and false negatives (FN) relative to the true sources. Precision and recall were computed as follows:

  $Precision=TPTP+FP,$

  $Recall=TPTP+FN.$

To validate the 4LCNN model with real EEG data, we employed SEP experiments due to the accumulation of evidence supporting the neural origin of SEP components. EEG and depth electrode recordings have confirmed that somatosensory information from median nerve stimulation propagates sequentially through the cervical spinal cord, medulla, thalamus, and somatosensory cortex at approximate intervals of 11, 13–14, 16–17, and 20 ms, respectively.²⁵ We administered electrical monophasic square wave constant current pulses for a duration of 200 ms from a clinical constant current stimulator (DS7A, Digitimer Ltd., Welwyn Garden City, United Kingdom) to stimulate the right median nerve. The electrode pair was placed on the wrist of the right hand. The current intensity was set to the sum of the motor and sensory thresholds for each participant, in line with the International Federation of Clinical Neurophysiology's recommendations.⁶⁴ Stimulation was applied 2000 times at a repetition rate of 4.98 Hz, with the entire session repeated four times, totaling 8000 stimulations. To maintain the participants' attention to the stimulus, pauses were included approximately every 45–90 s. During each session, stimulation was paused randomly for about 2 s between 4 and 7 times. The duration of each session was approximately 7 min, pauses included. Participants were instructed to count the number of pauses and report their count at the session's end. All participants accurately reported the number of pauses for all sessions.

EEG signals were recorded from 159 channels using an EEG amplifier (ActiCHamp Plus, Brain Products, Germany) with active electrodes (ActiCap Slim, Brain Products GmbH, Germany) at a sampling rate of 5000 Hz during the SEP experiment. Electrode impedances were maintained below 10 kΩ, significantly lower than the recommended impedance (<50 kΩ) for high-impedance EEG amplifiers. The same three-dimensional electrode coordinates employed for the head model creation were used, as the electrode position data were measured immediately prior to the SEP experiment. The EEG signals were bandpass filtered with a 4th order Butterworth filter at a cutoff frequency of 100–250 Hz.^65,66 This high-frequency component was selected to minimize baseline fluctuation and to separately evaluate components occurring closely at 13–14, 16–17, and 20 ms, corresponding to the medulla, thalamus, and primary somatosensory cortex components, respectively. SEPs from all participants exhibited distinct peaks corresponding to these components. We estimated the source activity of these components from the EEG spatial distribution map at the peak times using the 4LCNN model. The noise level of SEPs was evaluated using the SNR. The signal size was defined as the RMS of the amplitude across all electrodes at the peak timing of each component, while the noise level was defined as the mean RMS value during the pre-stimulus period (from −50 to −5 ms before the stimulation).

To further evaluate the 4LCNN model on the real data, we conducted ESI using the 4LCNN on MEG data and evaluated its accuracy with simultaneously recorded invasive brain signal recordings from SEEG and ECoG (Fig. 10). This study included three patients diagnosed with epilepsy.

We simultaneously recorded MEG, SEEG, and ECoG signals from the participants in a resting state (with eyes closed, relaxed, possibly asleep) for 5 min. This recording was repeated four to six times per participant, accumulating 20–30 min of data for each participant. In addition, we recorded “empty-room data” (with no subject present) for 5 min to assess the background noise in the MEG signals.

MEG data were recorded using a 160-channel whole-head MEG system (MEG Vision NEO, Yokogawa Electric Corporation, Tokyo, Japan). Participants lay supine in a magnetically shielded room with five head-marker coils attached to their face. The positions of these head-marker coils were measured before and after each MEG recording to determine the relative position and orientation of the MEG sensors to the participant's head. The MEG signals were sampled at 2000 Hz with a bandpass filter setting of 0.1–500 Hz. To co-register MEG data with individual MRI data, the three-dimensional locations of 50 points on each participant's scalp and facial surface were digitized (FastSCAN, Polhemus, Colchester, Vermont, USA).

SEEG and ECoG electrodes were implanted, and their locations in each participant are shown in supplementary material Figs. 9–11. We collected data from 24 to 54 SEEG contacts and 42–70 ECoG contacts, recorded at 2000 Hz by a digital EEG system (EEG-1200, Nihon Koden, Tokyo, Japan). All subcortical electrodes were referenced to average potential of two electrodes placed subcutaneously or on the cortex during the recordings. Initially, signals were filtered using a bandpass filter ranging from 1 to 60 Hz.

Pre- and post-implantation T1-weighted MRI scans were conducted with Discovery MR750 (GE Healthcare, Chicago, U.S.A.; pre: TR = 7.3 ms, TE = 2.8 ms, 1 mm isotropic resolution; post: TR = 7.7 ms, TE = 3.5 ms, 0.95 mm isotropic resolution) for EP1 and SIGNA Architect (GE Healthcare; pre: TR = 2300.0 ms, TE = 3.0 ms; post: TR = 10.0 ms, TE = 4.0 ms, 1 mm isotropic resolution) for EP2 and EP3.

MRI data were defaced using the AnonyMI algorithm⁶⁷ to reduce the risk of identification while conserving geometrical information before stored in a processing workstation in accordance with the ethical policy. After electrode implantation, a CT scan was also performed to ascertain the positions of the implanted electrodes with the Aquilion Precision scanner (Toshiba Medical Systems, Tokyo, Japan).

We created 4LCNN models for estimating MEG source activity using the same methodology as for EEG source estimation, with minor modifications tailored for MEG. The spatial alignment of MEG sensors and MRI data were based on the coordinates of the head surfaces digitized during the MEG experiment and the head surfaces extracted from T1-weighted MRI images. This alignment was achieved through rigid-body transformations to minimize the distance of these points to the head surface using the Brainstorm software. MRI and CT data were also aligned using the “ft_convert_coordsys” function in Fieldtrip. ECoG and SEEG electrode locations were then determined using MRI and CT data via a pipeline designed for merging anatomical images with intracranial electrode placement in Fieldtrip.⁶⁸ 

MEG signals were bandpass filtered within a frequency range of 1–60 Hz, and then re-referenced with a common average reference (CAR) technique.²⁸ Subsequently, independent component analysis (ICA) was performed using the “runica” function in the EEGLAB toolbox.⁶⁹ Correlation coefficients (r) between the time courses of the independent components (ICs) and the signals from SEEG and ECoG electrodes were calculated. For each IC, if r values exceeded 0.2 for any electrode, the electrode with the highest r was assumed to be most proximal to the source of the IC. We then computed the Euclidean distance between the estimated location of the IC and the corresponding electrode location.

We compared the performance of the 4LCNN model with ConvDip,⁸ a previously published CNN model for ESI, and two commonly utilized ESI methods—eLORETA²⁰ and the LCMV beamformer²¹—using the same set of simulations and experimental data. ConvDip was trained with the same forward M/EEG simulation data as the 4LCNN using MATLAB's Deep Learning Toolbox. The eLORETA and LCMV analyses were conducted using the Fieldtrip software. Both eLORETA and LCMV require noise covariance information among sensors. We utilized an identity matrix for the noise covariance for the simulated EEG data because white Gaussian noise was used in the EEG simulation. For MEG ICs, an identity matrix was also employed since the noise covariance among sensors within a single IC was indeterminate, following a recent study.³⁰ For SEP data, covariance matrices were derived from signals during the pre-stimulation period (−50 to −5 ms) in accordance with previous studies on evoked responses.^8,70 The regularization parameter for the noise covariance matrix in eLORETA and LCMV was set to SNR⁻¹ for the simulation and SEP data,⁵⁶ and a default value of 0.05, as per Fieldtrip's settings, was used for the single ICs of MEG. The noise level of individual ICs was evaluated using the SNR. The signal size was defined as the RMS of the activity time courses of the ICs extracted from the MEG signals using the unmixing vector obtained through ICA, while the noise level was defined as the RMS of the activity time courses derived from the empty-room MEG data using the same unmixing vector applied for the signal size calculation.

The normality of the datasets was evaluated using the Lilliefors test. For all ESI methods, estimated activity was normalized to the peak value across sources. With respect to the three performance metrics—localization error, spatial dispersion, and AUPRC—for the accuracy validation with simulation data, the differences between the ESI methods were statistically compared using a permutation test, a non-parametric method,⁷¹ across each SNR condition and metric. The permutation test was selected because the normality was not observed in all cases. P-values derived from the permutation test were adjusted for multiple comparisons using the false discovery rate (FDR) control method.⁷² 

We also explored the relationships of activation areas between the simulated and estimated data in the simulation dataset. The activation area was defined as the number of active sources with a value exceeding 0.5.

In the ESI analysis of SEP data from healthy participants, the spatial dispersion metric was computed, as this metric can be determined without true source activity information. Spatial dispersion was compared among the ESI methods using a paired t-test with FDR correction since normality was confirmed.

In the ESI analysis of resting state MEG data from individuals with epilepsy, both localization error and spatial dispersion were compared among ESI methods using permutation tests with FDR correction due to the lack of normality in the datasets.

See the supplementary material for detailed statistical results, individual data, and supportive additional data.

The work was supported by a Grant-in-Aid for Scientific Research (B) from the Japan Society for the Promotion of Science (JSPS) to H.Y. (No. 21H03340), ACT-X from the Japan Science and Technology Agency (JST) to H.Y. (No. JPMJAX22AN), and the JST-MOONSHOT program to K.N. and T.Y. (No. JPMJMS2012).

The authors have no conflicts to disclose.

Ethics approval for experiments reported in the submitted manuscript on animal or human subjects was granted. All participants provided written informed consent for participation in the study, which was conducted in accordance with the Declaration of Helsinki and approved by the Ethics Review Committee of The University of Tokyo (reference number 701-2) and Osaka University (reference number 14353).

Hikaru Yokoyama: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Funding acquisition (equal); Investigation (equal); Methodology (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Naotsugu Kaneko: Data curation (equal); Writing – original draft (equal); Writing – review & editing (equal). Noboru Usuda: Data curation (equal); Formal analysis (supporting); Methodology (supporting); Visualization (equal); Writing – review & editing (equal). Tatsuya Kato: Data curation (equal); Methodology (equal); Writing – review & editing (equal). Hui Ming Khoo: Data curation (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal). Ryohei Fukuma: Data curation (equal); Methodology (equal); Writing – review & editing (equal). Satoru Oshino: Data curation (equal); Methodology (equal); Writing – review & editing (equal). Naoki Tani: Data curation (equal); Methodology (equal); Writing – review & editing (equal). Haruhiko Kishima: Conceptualization (equal); Supervision (equal); Writing – review & editing (equal). Takufumi Yanagisawa: Methodology (equal); Supervision (equal); Writing – review & editing (equal). Kimitaka Nakazawa: Funding acquisition (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.