Multimodal neuroimaging with optically pumped magnetometers: A simultaneous MEG-EEG-fNIRS acquisition system

Multimodal neuroimaging plays an important role in neuroscience research. Integrated noninvasive neuroimaging modalities, such as magnetoencephalography (MEG), electroencephalography (EEG) and functional near-infrared spectroscopy (fNIRS), allow neural activity and related physiological processes in the brain to be precisely and comprehensively depicted, providing an eﬀective and advanced platform to study brain function. Noncryogenic optically pumped magnetometer (OPM) MEG has high signal power due to its on-scalp sensor layout and en- ables more ﬂexible conﬁgurations than traditional commercial superconducting MEG. Here, we integrate OPM-MEG with EEG and fNIRS to develop a multimodal neuroimaging system that can simultaneously measure brain electrophysiology and hemodynamics. We conducted a series of experiments to demonstrate the feasibility and robustness of our MEG-EEG-fNIRS acquisition system. The complementary neural and physiological signals si- multaneously collected by our multimodal imaging system provide opportunities for a wide range of potential applications in neurovascular coupling, wearable neuroimaging, hyperscanning and brain-computer interfaces.


Introduction
Multimodal neuroimaging uses the strengths of individual modalities and can thus depict brain activity from different perspectives based on various natures of multimodal signals, such as combining electrophysiological and hemodynamic signals to study neurovascular coupling. In practice, multimodal neuroimaging requires the appropriate fusion of compatible and complementary modalities to develop an effective platform for multidisciplinary applications in basic scientific research and clinical use.
Both magnetoencephalography (MEG) and electroencephalography (EEG) are noninvasive techniques that directly measure brain electrophysiological activity by recording either magnetic fields or electric potentials generated by neuronal currents from the scalp with millisecond temporal resolution. MEG source imaging is less affected by distortions caused by heterogeneous conductivity and geometric approximations also reduced the distance between sensors and the scalp ( Pfeiffer et al., 2020 )). The depth of the OPM sensors in the helmet slot can be adjusted to fit the unique head shapes of individual subjects. Moreover, the flexible sensor configuration improves the compatibility of OPM-MEG sensors with other imaging modalities. Specifically, because they are not constrained by the rigid helmet at the bottom end of the dewar, as in SQUID-MEG sensors, the layout of the OPMs and sensors of other modalities can be optimized for data acquisition (e.g., interlaced on a cap or in a helmet customized for an individual subject's head shape). However, like SQUID-MEG, OPM-MEG also needs to operate in a magnetically shielded environment (e.g., inside a magnetically shielded room (MSR) Hill et al., 2022 ) or a compact cylindrical shield ( Borna et al., 2020( Borna et al., , 2017He et al., 2019 ;Nardelli et al., 2020 )) to avoid confounding noises (e.g., geomagnetic field), which requires additional consideration of compatibility when combined with other imaging modalities.
To date, MEG and EEG have been separately used in conjunction with functional near-infrared spectroscopy (fNIRS), which is a noninvasive and portable method to measure cortical hemodynamics caused by neural activity and to investigate the neurovascular coupling associated with brain functions ( Czeszumski et al., 2020 ;Kassab et al., 2018 ;Mackert et al., 2008Mackert et al., , 2004Ou et al., 2009 ;Sander et al., 2007 ;Seki et al., 2012 ). Although both MEG and EEG directly measure the neural activity in the brain, their sensitivity to neuronal currents varies at different orientations and depths ( Lopes da Silva, 2013 ). Therefore, the combination of the complementary information provided by MEG and EEG is expected to provide a more accurate source localization of neural activity in the brain ( Chowdhury et al., 2015 ;Ding and Yuan, 2013 ;Henson et al., 2009Henson et al., , 2011. This can also facilitate their comparison and integration with fNIRS data, which is mainly contributed by the regional signals near the detector ( ∼ 1 cm resolution). However, to date, no work that combines OPM-MEG, EEG and fNIRS to simultaneously measure electrophysiological and hemodynamic signals has been reported.
In this study, we designed and implemented our homemade OPM-MEG, EEG and fNIRS devices and integrated these devices into a multimodal neuroimaging system to simultaneously record both brain electrophysiology and hemodynamics. We demonstrate the feasibility and robustness of our multimodal OPM-MEG, EEG and fNIRS acquisition system through a series of experiments. Finally, we discuss the advantages, potential applications, and future directions of our multimodal neuroimaging technique.

Multimodal system overview
A schematic representation of our prototype MEG-EEG-fNIRS multimodal system is shown in Fig. 1 . A cylindrical magnetic shield was used to generate the OPM-MEG working environment. OPM sensors were placed in a 3D-printed rigid helmet (shown in Supplementary Materials Fig. 1 ), which enables flexible sensor placement and an individualized close-to-scalp configuration. The OPM electronics, including the multichannel control and analog output systems, were installed beneath the scanning bed and connected to the sensors through the drag chain of the sliding bed. The EEG cap worn on the subject's head inside the rigid helmet was fitted with optical fibers of the fNIRS system to deliver near-infrared light for measuring hemodynamic responses. Two analog input modules (PXIe-4309, National Instruments Corp., Austin, TX, USA) were used to record and synchronize the signals from the OPM, EEG and stimulus channels and confounding bioelectrical signals of the eye blink/movement and cardiac activity. The stimulus PC was used to control the stimulation and feedback system (SR-DP 100 M & SA-9800E, Shenzhen Sinorad Medical Electronics Co., Ltd., Shengzhen, Guangdong, China) and deliver stimulus triggers.

OPM-MEG
The OPM-MEG system consisted of a single-ended cylindrical shield, OPM sensors and associated electronics. The cylindrical shield was based on our previous design ( He et al., 2019 ) and built for wholebrain coverage OPM-MEG. The shield included four layers of 1.5 mm thick permalloy (model 1J85, Advanced Technology & Materials Co., Ltd., Beijing, China) and a 10 mm thick aluminum alloy shell. The shield was oriented with its longitudinal axis perpendicular to the geomagnetic field to achieve an optimal DC shielding performance.
A schematic diagram of the OPM sensor is shown in Fig. 2 . In total, 25 sensors were prepared for the current study. The sensors were developed based on our previous two-beam configuration ( Sheng et al., 2017 ), with an optimized vapor cell and structure ( Allred et al., 2002 ;Griffith et al., 2010 ;Zhang et al., 2018 ). Moreover, our sensors are alloptical designed and consequently free from cross-talk effects (due to magnetic modulation in neighboring sensors). A borosilicate glass cubic cell with an inner edge length of 4 mm contained a small droplet of rubidium (Rb). The cell was heated to approximately 150 °C to enter the spin-exchange relaxation-free (SERF) regime. One compact diode module emitted a circularly polarized pumping beam that was tuned to the D1 line of Rb in the lower plane. Then, the pumping beam was reflected upward into the vapor cell. The other laser diode module emitted a probing beam that was perpendicular to the pumping beam in the cell Fig. 3. Framework of the CW-fNIRS system. The system included 4 source boards, 6 detector boards, 3 control boards and one trigger input unit. Each source board included 2 source units, and each source unit included 3 laser diodes with different wavelengths (780 nm, 808 nm and 850 nm). Each detector board consisted 3 detector units. The electronic and optical components are indicated in black and blue, respectively. FPGA: field programmable gate array; CPLD: complex programmable logic device, DAC: digital-toanalog converter, LD: laser diode module, including single wavelength laser diode and laser drive; APD: avalanche photodiode; AMP: transimpedance amplifier; ADC: analog-to-digital converter.
plane. The probing beam was blue-detuned approximately 0.1 nm from the D2 line of Rb and reflected into the vapor cell. As a result, the OPM was sensitive along the longitudinal direction of the sensor (as shown by the arrow in Fig. 2 ). Then, the probing beam that passed through the vapor cell was reflected to a polarizing beam splitter before being detected by a pair of differential photodetectors. The output signal was pre-amplified and transmitted to customized electronics. A small current coil was attached to the vapor cell to calibrate the output voltage to the magnetic field. A flexible printed circuit board (PCB) installed at the bottom of the sensor was used for data communication and control. Several shielded and twisted-pair cables were used to transmit signals between the flexible PCB and the OPM electronics. The OPM electronics automatically optimized the sensor control parameters (including heating temperature, laser voltage and current) according to sensor feedback, completed the calibration process, and output the parameters indicating the state of the sensor, analog signals, and digital signals. The OPM control computer enabled multichannel monitoring in real time, controlled all sensor parameters and calibrated the OPM data according to the calibration information measured before the experiment.

EEG
Following the conventional EEG configuration ( Muller-Putz, 2020 ), we designed and customized 32 gel-based passive electrodes and integrated them into an EEG cap. To be compatible with our OPM-MEG, the residual magnetic field of the EEG cap needs to be less than 0.5 nT, and no additional noise should be introduced by the EEG cap. Each EEG electrode consisted of an Ag/AgCl powder sintered cylinder (3 mm diameter and 1 mm thickness) with a silver-plated copper wire (1.5 m in length). Holes were pieced in an elastic cloth cap to attach 32 electrodes with a 3D-printed base at specific locations according to the 10-10 system ( Oostenveld and Praamstra, 2001 ). The reference electrode was placed at the position between Cz and CPz on the cap (i.e., the CCPz position) and provided a reference potential for the other electrodes. The electrode wires were welded to a single medical push-pull connector, which was first connected to the impedance detection circuit to measure the impedance and subsequently to the analog input modules. EEG signals were recorded with a sampling rate of 3 kHz, a 24-bit measurement resolution and an input range of ±10 V.

fNIRS
Continuous wave (CW) fNIRS with 8 sources and 18 detectors was designed and used in this study (shown in Fig. 3 ). In each source unit, three 12-bit, 95-ks/s digital-to-analog converters (DACs) received the digital light intensity signal from the complex programmable logic device (CPLD) and output an analog voltage to three laser diode (LD) modules with different wavelengths, including a single-wavelength laser diode and a laser drive in each module. The three lasers emitted by the LD modules were transmitted through three glass fibers (length/diameter = 3000 mm/3.2 mm) and combined into a single optode (shown in Supplementary Materials Fig. 2 ) to beam near-infrared light toward the scalp. The wavelength and maximum instantaneous power of the three lasers were 780 nm and 25 mW, 808 nm and 13 mW, and 850 nm and 25 mW at the source optodes. The light signal emitted from the brain was collected by detector optodes and subsequently transmitted through a glass fiber (length/diameter = 3000 mm/2.5 mm) and an optical filter with a 770-870 nm bandpass filter before being detected by a silicon-based avalanche photodiode (APD). A transimpedance amplifier (AMP) received the analog current from the APD and output an analog voltage. Then, the analog voltage was processed through a low-pass filter and a gain circuit. An analog-to-digital converter (ADC) (24-bit, 105-ks/s) received the processed analog voltage and sent the converted digital signal to the CPLD.
A time-sharing mechanism was used to acquire light intensity signals from different channels. During one sampling cycle, different source units emitted light in a fixed order. Within each source unit, the three laser diodes sequentially emitted light with their unique wavelengths (0.5 ms duration for each light; 0.5 ms gap between two lights). This fixed emission sequence was used to identify the light source recorded by the detectors. The source and detector units were paired to form the hemoglobin concentration measurement channels. Before data acquisition, a light intensity dynamic test was conducted to determine the emitted light power. The emission power of the light source was gradually increased from 0 to the maximum power. The optimal emitted power was chosen to maintain the detector sensitivity within the range of interest based on feedback from the detector.

Data collection & analysis
All experiments involving human subjects in this study were approved by the Institutional Review Board of Peking University. All subjects who participated in this study provided written informed consent.

Multimodal imaging system testing experiments
We conducted testing experiments to validate the performance of our multimodal imaging system. For the OPM-MEG system, we evaluated the performance of the cylindrical shield and measured the noise spectrum of the OPMs for quality control. The finite element method (FEM) was used to calculate the shielding performance in the geomagnetic environment. The distribution of the residual magnetic field around the helmet (i.e., the region of interest, ROI) in the shield was measured by a fluxgate (CTM-6 W, National Institute of Metrology, Beijing, China) in 5-cm grids. The fluxgate was reversed at several grid points to correct the offset. To test the noise level of the OPM, four OPMs were placed near each other and aligned in the same direction in a small six-layer cylindrical shield with an inner diameter of 42 cm. The output analog OPM data were sampled by the ADC at a sampling rate of 3 kHz. The noise spectrum was estimated by Welch's power spectral density using MATLAB software (The MathWorks, Inc., Natick, MA).
We also investigated the electronic performance and magnetic compatibility of the EEG cap and electrodes. The electrodes were immersed in 0.9% physiological saline for 30 min before measuring the impedance. The residual magnetic field induced by the EEG cap was quantified based on the difference in the readings of the OPM sensors before and after rapidly inserting the EEG cap into the OPM helmet. The interference caused by the EEG cap was also estimated based on the change in OPM noise spectrum after the EEG cap was placed in the OPM helmet.
The source intensity drift and detector linearity were evaluated to assess the performance of the CW-fNIRS system. The source and detector optodes of the CW-fNIRS system were inserted into a polyethylene plastic phantom in a darkroom to measure the source intensity drift. To measure the detector linearity, an optical attenuator was added to the optical path to change the intensity of the emitted light, and the optical power received at the detector was measured by an optical power meter.

Evaluation of the single-modal imaging performance
We conducted an auditory evoked field (AEF) experiment to evaluate the performance of the OPM-MEG system. A pure tone auditory stimulus (1 kHz, 100 ms duration) was presented 200 times with randomized interstimulus intervals (between 0.8 s and 1.3 s) to a healthy male subject. Twenty-five OPM sensors were placed on the helmet for better source imaging quality and covered the right frontal and temporal lobes. The insertion depth of each OPM sensor was adjusted according to the shape and position of the subject's head, so the distance between the center of the vapor cell and the scalp was as low as 6 mm. The subject was instructed to keep his head still during the experiment. A 3D optical scanner (HandyScan Black Elite, Creaform Inc., Levis, Canada) was used to collect the spatial information of both OPM sensors and the subject's face for co-registration with MRI structural images before the scanning bed was slid into the shield. The structural MRI data were collected on a 3-T Siemens Prisma MRI scanner (Siemens Healthineers, Erlangen, Germany) at Peking University, Beijing, China. T 1 -weighted brain images were acquired using an MPRAGE sequence (TR/TE = 2530 ms/2.98 ms and voxel size 0 . 5 × 0 . 5 × 1 m m 3 ) with a 64-channel head coil.
Similarly, we conducted an auditory evoked potential (AEP) experiment to evaluate EEG performance. A pure tone auditory stimulus (1 kHz, 100 ms duration) was presented 200 times with a fixed interstimulus interval of 0.5 s to the same subject in the AEF experiment. The impedance between the electrodes and the scalp was reduced to less than 10 kΩ with abrasive gel.
The raw time series of the OPM-MEG and EEG were bandpass filtered between 1 and 30 Hz. Then, the AEFs and AEPs were calculated by averaging good trials from − 100 to 500 ms relative to the stimulus onset after baseline correction. Noisy trials were excluded based on predefined thresholds and additional visual inspection. T 1 -weighted images were reconstructed and segmented into scalp and brain tissues using the FreeSurfer recon-all pipeline ( https://surfer.nmr.mgh.harvard.edu/ ). The co-registration between MRI and OPM-MEG was conducted using an automatic algorithm ( Gu et al., 2021 ). Source reconstruction was conducted using the minimum-norm method implemented in the Brainstorm toolbox ( https://neuroimage.usc.edu/brainstorm/ ).
A block-design visual experiment was conducted to evaluate the fNIRS performance. Ten 20 s blocks of visual stimuli and ten 30 s resting blocks (with a uniform black screen) were alternately presented to the same subject in the AEF experiment. In each block task, a checkerboard with 8 2 subdivisions was reversed at a frequency of 7.5 Hz on the screen. Nine optodes (5 sources and 4 detectors) were distributed in a 3 × 3 grid array with a distance of 3 cm between adjacent optodes of sources and detectors. The five source optodes were placed at the four corners and center of the grid, while the four detector optodes were placed at the midpoint of each edge of the grid. The grid was centered near the O2 position on the EEG cap and mainly covered the visual cortex in the right occipital lobe.
The raw fNIRS time series were preprocessed (bandpass filtered between 0.01-0.5 Hz, motion corrected, etc.) with the Homer2 toolbox ( https://homer-fnirs.org ). The evoked changes in cerebral oxygenation and hemodynamics were obtained based on the signals averaged across task blocks from − 10 to 40 s relative to the stimulus onset. Noisy segments were excluded based on predefined thresholds and additional visual inspection.
In the single-modal imaging experiments described above, the signalto-noise ratio (SNR) was estimated by dividing the peak-to-peak change in the evoked response by the standard deviation during the prestimulus baseline ( Boto et al., 2017 ).

Multimodal imaging experiments
We conducted a block-design visual evoked experiment on the simultaneous MEG-EEG-fNIRS acquisition system. Five healthy subjects (3 males and 2 females, right-handed, aged 25 . 4 ± 6 . 1 years) were recruited. In this experiment, ten 20 s visual stimulus blocks and ten 30 s resting blocks (with a uniform black screen) were alternately presented. In each stimuli block, a checkerboard with 8 2 subdivisions was reversed at a frequency of 7.5 Hz on the screen (see Supplementary Materials Fig. 3 for the setup of the multimodal visual evoked experiment).
To capture multimodal visual evoked effects, 18 OPMs were distributed on the helmet to cover the right visual cortex. The insertion depth of the OPMs in the helmet was adjusted to minimize the distance between each sensor and the scalp. Nine fNIRS optodes were attached to the EEG cap and placed in the same manner as in the single-modal fNIRS experiment (see Supplementary Materials Fig. 2 for more details). The impedance of the EEG electrodes was reduced to less than 10 kΩ using the abrasive gel. The OPM and EEG data were acquired at a frequency of 3 kHz, while the fNIRS data were acquired at a frequency of 50 Hz. The time series of the OPM-MEG and EEG were split into ten segments that covered a time window from − 10 to 40 s relative to the stimulus onset. These segments were visually inspected and subsequently transformed into the frequency domain with a short-time Fourier transform. The power spectrums were calculated and averaged using MATLAB software. The fNIRS data were processed similarly to those in the singlemodal fNIRS experiment. . 4 shows the magnetic field inside the cylindrical shield. The residual magnetic field in the ROI was consistent with the simulation results. In the ROI, the static magnetic flux density norm and peak-topeak time variation were less than 2 nT and 0.4 nT/hour, respectively. Fig. 5a shows a picture of the compact OPM sensor. Fig. 5b shows the noise spectrum of four adjacent OPMs that simultaneously operated. The OPM sensors had a sensitivity of 15 . 0 ± 2 . 1 fT ∕

√
Hz at 10 Hz with a uniform response, a bandwidth of 90 Hz and a dynamic range of ± 0.8 nT after calibration. In summary, the shield provided a large and homogeneous ROI with acceptable residual magnetic fields for the multichannel OPM-MEG system. The OPM sensors had low background noise, high sensitivity, and a small cross-sectional area, which enables the sensors to be configured in a high-density array.
In the EEG system, the impedance between the electrodes in 0.9% physiological saline was less than 0.1 kΩ. The residual magnetic field induced by the EEG cap was less than 50 pT. There was no significant change in power spectrum density of the OPMs after the EEG cap was inserted into the helmet in the ROI of the shield. These results demonstrate that our customized EEG system is stable and compatible with MEG. In the fNIRS system, the fluctuation of the detected signal was less than 1% in one hour when the source light had a constant intensity. The fitting linearity (R 2 ) between the real light intensity received by the detector and the voltage processed by the system reached 0.999 in the light intensity range in our experiments.  The AEF elicited by presenting pure tone stimuli that was recorded by the OPM-MEG system (as shown in Fig. 6a ) was consistent with previously reported results ( Jacobson, 1994 ). The typical peaks (N1m, P2m and the sustained field) were clearly visible. The topographic distribution of the N1m component exhibited a bipolar pattern centered in the middle of the right temporal lobe, which was further sourcelocalized at the primary auditory cortex. These results show the robust performance of our multichannel OPM-MEG system and the validity of our data analysis procedure. The AEP elicited by presenting pure tone stimuli (as shown in Fig. 6b ) demonstrated the robustness of our customized EEG system. The typical peaks (N1, P2) were clearly identified and corresponded to previously reported results ( Cone-Wesson and Wunderlich, 2003 ). Fig. 6c shows the hemodynamic responses to visual stimuli recorded by the fNIRS system. These curve shapes are consistent with previous reports of responses to visual stimuli ( Villringer and Chance, 1997 ), which validates the acquisition process and data analysis procedure of our fNIRS system. The SNRs of AEFs detected by OPM-MEG and those of AEPs detected by EEG were 21 . 6 ± 2 . 9 (mean ± standard error) and 14 . 8 ± 2 . 7 , respectively (the SNR was calculated from 11 OPM sensors and 8 EEG channels in the right temporal lobe). The SNR of the oxyhemoglobin response recorded by the fNIRS channel in Fig. 6c was 15.6. Fig. 7a shows the paradigm of the multimodal block-design visual experiment. Fig. 7b shows the multimodal responses of one representative subject in two runs recorded by the MEG-EEG-fNIRS system (see Supplementary Materials Fig. 4 for the results of the other four subjects). In all subjects, the power density of the steady-state visual evoked field/potential (SSVEF/SSEVP) at 15 Hz increased immediately after the continuous visual stimulus was presented, plateaued after 4-8 s, subsequently decreased immediately after the stimulus offset, and returned to baseline levels within 1-3 s. These results were consistent with previously reported SSVEF/SSVEP responses ( Thorpe et al., 2007 ). In the hemodynamic responses recorded by the fNIRS system, the oxyhemoglobin concentration began to increase 2-5 s after the stimulus onset and plateaued at 9-14 s. The concentration began to decrease at 2-8 s after the stimulus offset and returned to baseline levels in 10-14 s. The oxyhemoglobin concentration had a similar trend to the total hemoglobin concentration and an opposite trend to the deoxyhemoglobin concentration. These results were also consistent with previously reported hemodynamic responses to visual stimuli ( Villringer and Chance, 1997 ).

Discussion
In this study, we designed and implemented our homemade OPM-MEG, EEG and fNIRS systems and integrated them into a multimodal neuroimaging system. We conducted a series of experiments and analyses to evaluate the performance of our multimodal system. The robust results demonstrate the feasibility of combining these modalities to measure brain electrophysiology and hemodynamic responses simultaneously, show that our multimodal imaging devices are compatible, and also validate our data acquisition and analysis procedures.
OPM-MEG with either a compact cylindrical shield or a MSR is a promising future direction of MEG for both scientific research and clin- Power spectrum density of four typical OPMs measured in a small six-layer cylindrical shield with an inner diameter of 42 cm. The OPMs were attached to each other and aligned in the same direction to simulate a high-density arrangement. The power spectrum was divided by a pre-measured normalized frequency response to obtain the absolute sensitivity. ical applications. At present, rubidium is a popular type of alkali metal in SERF-based OPM-MEG sensors ( Borna et al., 2017 ;Boto et al., 2017 ). To enter the SERF regime, the vapor cell of the OPM sensor needs to be heated to approximately 150 °C to achieve the desired Rb vapor density. The heat from the sensors is an issue that needs to be addressed to guarantee a comfortable experimental environment, particularly for highdensity multichannel OPM-MEG. The heat issue in our experiments (up to 25 sensors) was properly solved by pasting aerogel with a thickness of 0.5 mm onto the bottom end of each OPM sensor attached to the scalp to prevent heat conduction. All subjects reported that they felt comfortable during the experiment ( ∼ 40 min). However, for whole-head high-density OPM-MEG (e.g., 80-100 channels), additional ventilation or a cooling system is required for better heat dissipation.
Our cylindrical shield significantly reduces the cost and amount of shielding materials compared with a typical MSR while achieving the same level of shielding performance without an active compensation module. In addition, the compact cylindrical shield has fewer requirements on space and load bearing, so it is easier to implement and more compatible with different site environments. Although the magnetic field fluctuations in the ROI of our cylindrical shield were less than 0.4 nT per hour, the irregular fluctuations in the magnetic field still introduced additional environmental noise into the signals recorded by the OPM sensors. Therefore, an active compensation module based on closed-loop feedback can be incorporated with our cylindrical shield to further attenuate the magnetic field fluctuations in the ROI and enhance our system's resistance to complicated electromagnetic noise in different environments.
In this study, experiments and related analyses were conducted to evaluate the basic performance of our MEG-EEG-fNIRS system. More elaborate experimental paradigms and advanced data processing methods can be adopted to fully utilize this multimodal imaging system. Although the results of our the single-modal and multimodal imaging experiments were clear, repeatable, and consistent with previous studies, the signal-to-noise ratio can be further improved by adopting more sophisticated data processing methods for denoising or reducing artifacts.
The multimodal MEG-EEG-fNIRS system can simultaneously measure electrophysiological activity and hemodynamic responses in the brain; thus, it has potential value for studying the substrates and mechanisms underlying neurovascular coupling in both healthy individuals and those with pathological conditions (such as hypertension, Alzheimer's disease and ischemic stroke) ( Girouard and Iadecola, 2006 ). In addition, our MEG-EEG-fNIRS system can be adapted for a highperformance MSR. Combined with magnetic field compensation technology, head position tracking technology and relevant online data processing methods Holmes et al., 2018 ;Rea et al., 2021 ;Seymour et al., 2021Seymour et al., , 2022, these flexible OPM-MEG, EEG and fNIRS modules allow subjects to freely and naturally move in a magnetically shielded space. Moreover, neurovascular coupling during social interactions can be investigated through hyperscanning experiments in MSR ( Ahn et al., 2018 ). Taken together, these methods and potential improvements can greatly facilitate and expand the use of our system in research areas such as wearable neuroimaging  and hyperscanning ( Czeszumski et al., 2020 ). Moreover, OPM-MEG, EEG, and fNIRS have been used as brain-computer interfaces (BCIs) respectively in many studies ( Ahn and Jun, 2017 ;Liu et al., 2021 ;Wittevrongel et al., 2021 ). The simultaneous acquisition of complementary information from multimodal signals owns the potential to Responses of one subject (subject 3) in two runs, which were recorded by sensors of three modalities at approximately identical position (i.e., position O2 in the 10-10 system): power density of the steady-state visual evoked field at 15 Hz recorded by an OPM sensor (left column), power density of the steady-state visual evoked potential at 15 Hz recorded by an EEG electrode (middle column) and visual evoked oxyhemoglobin concentration changes recorded by an fNIRS channel (right column). The gray area indicates the presentation of visual stimuli. PD: power density; HbO: oxyhemoglobin; HbR: deoxyhemoglobin; HbT: total hemoglobin. accelerate the training process and improve the decoding accuracy in BCI applications.
In summary, our multimodal MEG-EEG-fNIRS system integrates the advantages of MEG, EEG, and fNIRS and can simultaneously measure brain electrophysiology and hemodynamics, which enables us to probe neural systems, vascular systems and their interactions underpinning brain function and dysfunction. Our multimodal imaging system can also be potentially implemented in an MSR, providing opportunities for novel applications in wearable neuroimaging, hyperscanning and braincomputer interfaces.

Data and code availability statement
All data and code for data processing are available on request from the corresponding author.

Declaration of Competing Interest
The authors report no declarations of interest.