==== Front bioRxiv BIORXIV bioRxiv Cold Spring Harbor Laboratory 37398157 10.1101/2023.05.26.542419 preprint 1 Article No Replication of Direct Neuronal Activity-related (DIANA) fMRI in Anesthetized Mice Choi Sang-Han 1 Im Geun Ho 1 Choi Sangcheon 2 Yu Xin 2 Bandettini Peter A. 3 Menon Ravi S. 4 Kim Seong-Gi 15 1 Center for Neuroscience Imaging Research, Institute for Basic Science, Suwon, Republic of Korea 2 Athinoula A. Martinos Center for Biomedical Imaging, Department of Radiology, Harvard Medical School, Massachusetts General Hospital, Charlestown, Massachusetts, USA 3 Section on Functional Imaging Methods and Functional MRI Facility, NIMH, NIH, Bethesda, MD, USA 4 Centre for Functional and Metabolic Mapping, Robarts Research Institute, Western University, London, Ontario N6A 5B7, Canada 5 Department of Biomedical Engineering, Sungkyunkwan University, Suwon, Republic of Korea Author contributions S.-G.K. and X.Y. initiated the study, S.C. and X.Y. provided a line scanning pulse sequence, G.H.I. and S.-H.C. performed the experiments, S.-H.C. analyzed data and simulated noise, P.A.B. and R.S.M. provided insights of human fMRI, S.-G.K. wrote an initial draft, and all authors edited and approved the manuscript. * Correspondence to seonggikim@skku.edu 29 5 2023 2023.05.26.542419https://creativecommons.org/licenses/by-nd/4.0/ This work is licensed under a Creative Commons Attribution-NoDerivatives 4.0 International License, which allows reusers to copy and distribute the material in any medium or format in unadapted form only, and only so long as attribution is given to the creator. The license allows for commercial use. nihpp-2023.05.26.542419.pdf Toi et al. (Science, 378, 160–168, 2022) reported direct imaging of neuronal activity (DIANA) by fMRI in anesthetized mice at 9.4 T, which could be a revolutionary approach for advancing systems neuroscience research. There have been no independent replications of this observation to date. We performed fMRI experiments in anesthetized mice at an ultrahigh field of 15.2 T using the identical protocol as in their paper. The BOLD response to whisker stimulation was reliably detected in the primary barrel cortex before and after DIANA experiments; however, no direct neuronal activity-like fMRI peak was observed in individual animals’ data with the 50–300 trials used in the DIANA publication. Extensively averaged data involving 1,050 trials in 6 mice (1,050×54 = 56,700 stimulus events) and having a temporal signal-to-noise ratio of 7,370, showed a flat baseline and no detectable neuronal activity-like fMRI peak. Thus we were unable to replicate the previously reported results using the same methods, despite a much higher number of trials, a much higher temporal signal-to-noise ratio, and a much higher magnetic field strength. We were able to demonstrate spurious, non-replicable peaks when using a small number of trials. It was only when performing the inappropriate approach of excluding outliers not conforming to the expected temporal characteristics of the response did we see a clear signal change; however, these signals were not observed when such a outlier elimination approach was not used. direct neuronal activity fMRI sample size noise distribution ==== Body pmc1. Introduction Blood oxygenation level-dependent (BOLD) functional magnetic resonance imaging (fMRI) has revolutionized how neuroscientists investigate human brain functions and networks (Kwong et al., 1992; Ogawa et al., 1992; Bandettini et al., 1992). However, BOLD fMRI measures hemodynamic responses as a surrogate of neuronal activity, thus spatial and temporal resolution is closely dependent on BOLD sensitivity and neurovascular coupling. To further understand brain functions, it would be ideal to measure causality and sequences of neuronal events by fMRI. Initial approaches have been to utilize differences of BOLD fMRI onset times between layers or regions (Menon and Kim, 1999; Silva and Koretsky, 2002; Bellgowan et al., 2003; Yu et al., 2014; Jung et al., 2021). The rationale for our ultra-high resolution, capillary-specific studies was that since the contribution of capillaries to BOLD fMRI increases with magnetic field strength, and these capillaries are in close proximity to neurons, early hemodynamic responses at ultrahigh fields could reflect changes in neuronal activity. However, while this approach is feasible in ultrahigh fields with extensive averaging, it may be beyond the capabilities of current human MRI systems. Recently, Toi et al. reported that direct imaging of neuronal activity (DIANA) can be achieved by high temporal resolution fMRI in the anesthetized mouse model at 9.4 T (Toi et al., 2022). When the whisker pad electric stimulus was applied, the peak intensity of ~0.15% was detected at ~12 ms in the thalamus and ~25 ms in the primary somatosensory cortex following the stimulus pulse. When an optogenetic pulse was applied in the primary sensory cortex (S1), the peak intensity of ~0.3% was observed at ~15 ms in the S1 and ~25 ms at the thalamus after the stimulus onset. These DIANA data were shown to have low variability among animals and were highly consistent with electrophysiology latency data. Although low temporal signal-to-noise ratio (tSNR) is expected due to high spatial (220×220×1000 μm3) and temporal resolution (5 ms) in conjunction with an intrinsically noisy acquisition method (FLASH), 10.8-s long trials were averaged only 40 times for whisker stimulation and even less (13–20 times) for optogenetic stimulation. This breakthrough fMRI approach has an unusual combination of high spatial resolution, high temporal resolution, high detectability, low inter-animal variation, and most importantly, direct neuronal detection. Many laboratories around the world have been attempting to reproduce DIANA signals in various experimental conditions. In preliminary human studies (Hodona et al., 2023), DIANA signals were not observed, possibly due to species difference (human vs. mouse) or subject status (awake vs. anesthetized). Although data in the original publication (Toi et al., 2022) are highly convincing in multiple experimental conditions, it is critical to independently reproduce these DIANA findings. Thus, we repeated the reported whisker stimulation experiments in virtually the same mouse model used in Toi et al. (Toi et al., 2022). Compared to this study, two improvements were made; 1) continuous infusion of anesthetics (You et al., 2021) rather than intermittent bolus injection of anesthetics (Shim et al., 2018) for maintaining stable animal physiology during fMRI studies, and 2) the use of 15.2 T rather than 9.4 T for enhancing SNR (Jung et al., 2019). All other parameters were the same as in Toi et al. (2022). 2. Results Six anesthetized mice were used for fMRI studies of whisker pad stimulation at 15.2 T (Fig. 1A). For anesthesia, a mixture of ketamine and xylazine was initially injected intraperitoneally (IP), and was continuously infused intravenously (You et al., 2021; Shim et al., 2022) for maintaining stable animal physiology. Note that Toi et al. (2022) used the bolus IP delivery of supplementary anesthetics when needed (Shim et al., 2018; Jung et al., 2019; Shim et al., 2020; Dinh et al., 2021; Jung et al., 2021), which modulates anesthetic depth during experiments (Jaber et al., 2014; Shim et al., 2018). For whisker stimulation, electrodes were placed on the mouse’s right whisker pad (You et al., 2021). Stimulation parameters were at a pulse width of 0.5 ms and current intensity of 0.5 mA. 2.1. BOLD fMRI responding to whisker stimulation at 15.2T All fMRI data were acquired at 15.2 T for enhancing sensitivity. One 1-mm thick slice covering the primary somatosensory barrel field (S1BF) was chosen using scout BOLD fMRI studies to map S1BF. To ensure that anesthetized animal’s condition enabled reliable detection of the stimulus over the duration of the imaging session, BOLD fMRI studies were performed using standard gradient-echo (GE) echo planar imaging (EPI) with repetition time (TR) of 1 s and echo time (TE) of 11.5 ms before and after DIANA experiments (Figure 1B). These control experiments to confirm neuronal activity during the imaging session utilized 200-s runs with two blocks of 20-s whisker stimulation. One animal’s fMRI data are presented in Figure 1C–D. The BOLD map responding to whisker pad stimulation at 4 Hz (Fig. 1C for one mouse) shows reliable activation in the contralateral S1BF area and contralateral thalamus. Two 5×5-voxel regions of interest (ROI) were chosen, an active ROI at the contralateral S1BF and an inactive ROI in the ipsilateral thalamus for further time course analyses. The post-DIANA BOLD fMRI response in the active ROI was higher than the earlier pre-DIANA BOLD response (Fig. 1D for one animal’s time courses, and Fig. 1E for individual data), which is consistent with our previous time-dependent BOLD fMRI studies (You et al., 2021). This indicates that BOLD fMRI responses are intact pre- and post-DIANA, confirming the presence of neuronal activity, and thus the animals’ physiological state should allow for the detection of a stimulus in the DIANA experiments. To ensure that the DIANA acquisition method, a shuffled line scanning k-t pulse sequence originally developed by Silva and Koretsky (Silva and Koretsky, 2002), was working correctly, we performed BOLD fMRI with TR/TE of 100/11.7 ms, flip angle = 17° in the S1BF (Ernst angle), and in-plane resolution of 0.22 × 0.22 mm2. Acquisition of 160 frames was repeated for 54 phase-encoding steps, resulting in experimental run time of 16 s × 54 = 864 s. The BOLD response to 1-s whisker stimulation was reliably detected - even from one run of 54 stimulus events (Fig. 1F–G). An average signal change in the contralateral S1BF area was 0.49 ± 0.27% (SD, n = 5), while small or negligible response was observed in the ipsilateral thalamus ROI (Fig. 1G–H). Our data indicate that the shuffled k-t pulse sequence was implemented correctly. We do note that the k-t approach is intrinsically more noisy despite 1 s × 54 phase-encoding steps = 54 s of total stimulation for a single BOLD-LS run versus 20 s × 2 runs = 40 s in BOLD-EPI, leading to lower correlation and p values for the k-t approach. 2.2. DIANA fMRI with high temporal resolution at 15.2T: Sensitivity and number of averages To reproduce DIANA findings in anesthetized mice with the same imaging parameters (Toi et al., 2022), the shuffled line scanning k-t pulse sequence was used with TR/TE of 5/2 ms and spatial resolution of 0.22 × 0.22 × 1 mm3. Acquisition of 40 (200 ms) or 200 frames (1000 ms) was repeated at 54 phase-encoding steps, lasting 10.8 or 54 s for each run. In each mouse, at least 50 runs were repeated for each experimental condition (see Table 1 for experimental design). Since an average of only 40 10.8-s trials at 9.4 T showed a 0.17% peak response in the S1BF ROI at 25 ms after the stimulus in the original publication (see Fig. 1D in Toi et al., 2022), we expected to easily detect the DIANA peak with 50 averages at 15.2 T given that tSNR of thermal-noise-dominant high-resolution fMRI increases with magnetic field strength as (15.2/9.4)1.65 = 2.2 (Pohmann et al., 2016; Uludağ and Blinder, 2018). To measure tSNR, data with 1000-ms duration (200 frames) and 50 trials were analyzed in detail (Mouse #3). tSNR was calculated by standard deviation divided by the signal mean over the 200 frames. Within the S1BF ROI, average voxel-wise tSNR was 32 for a single trial (ranging between 28 and 42 in 6 animals), and 230 for 50 trials-averaged data (Fig. 2A). The activation map obtained using a cross-correlation analysis with the DIANA response function (Toi et al., 2022) did not show any visible activation cluster in the S1BF (Fig. 2B). This may be explained by the low sensitivity of direct neuronal-related fMRI response on a single voxel level. To increase tSNR and detectability, time courses were obtained from 5×5-voxel ROIs using a moving average of 3 frames (Fig. 2C). Assuming that all data points are independent, spatial averaging and temporal smoothing increases tSNR by (number of voxels × number of temporal averages)1/2. Trial-wise tSNR across 200 frames was 279.1±23.5 (SD, n = 50 trials), and tSNR of the averaged time series was 1585 for the entire 200 frames and 2628 for the initial 20 frames. If stimulus-related responses were added to the variance, tSNR obtained from the entire frames was underestimated. At first glance, the averaged time course showed a direct neuronal activity-like peak of 0.2% in the S1BF ROI right after somatosensory stimulation and another peak of −0.15% at 340 ms. When pre-stimulus data points (20 frames × 50 trials) were compared with data in the peak (50 trials), both peaks were statistically significant (p = 8.7e-05 and 0.0025, two-sample t-test). However, a time to positive peak of 10 ms (Fig. 2D inset time course) did not match with the neuronal time to peak of 25 ms (Toi et al., 2022), thus our level of certainty that this peak reflects neuronal activity is low. The statistically significant peak may be a consequence of noise in a limited sample size. When all time course data points (50 trials × 200 frames) were considered, a Gaussian noise distribution was observed (Fig. 2E). However, when only 50 trials of the peak frame (10 ms post-stimulus) were considered, a skewed distribution was observed in a histogram (Fig. 2F), indicating that a noise peak can appear statistically significant due to an insufficient, biased sample. If an observed peak is genuine, then it should be reproducible across animals and be differentiated better from the noise level when more averaging is performed. Thus, an increased number of repeated trials is needed to determine whether a genuine neuronal activity-related peak can be detected. 2.3. Extensively averaged DIANA fMRI data: No observation of direct neuronal activity-related peak All repeated runs in each mouse were stimulation-pulse-locked, leading to 50 – 300 trials in each animal (Table 1). There was no difference of baseline variations between runs with and without radio frequency (RF) spoiling. In each subject, time courses of all trials (50 – 300) were obtained from the contralateral S1BF (Fig. 3A) and inactive ipsilateral thalamic ROIs (Fig. 3C). Then, an averaged time course was obtained by averaging all repeated trials in each subject (Fig. 3B and D). ROI-wise tSNR for individual animals with 50 – 300 averages ranged between 1950 and 4583 in the S1BF ROI and between 1500 and 2731 in the inactive ipsilateral thalamus ROI. Even in averaged time series with high tSNR values, no DIANA-like activity was observed in any of the animals. Note that the statistically significant positive peak observed in Figure 2 proved to be insignificant when 100 trials were averaged (3rd row in Fig. 3B). To enhance the sensitivity further, time courses were averaged from all 1050 trials in 6 mice (1050 × 54 = 56,700 pulse events) (Fig. 3E and G) and from individual animal’s averaged time courses (Fig. 3F and H). When 1050 trials were averaged, tSNR was 7370 for the S1BF ROI. Even when extensive averaging was achieved, no direct neuronal activity-like fMRI peak was detected. 2.4. Artifactual DIANA peak can be generated by exclusion of presumably bad trials. Given our inability to detect a DIANA signal in data with extremely high tSNR and extensive averaging, we investigated how one might observe a DIANA peak when dealing with more limited sampled data. In anesthetized animal fMRI studies with the intermittent IP bolus injection of anesthetics as Toi et al. (2022) adopted, anesthesia depth is modulated (Jaber et al. 2014; Shim et al., 2018). At a deep anesthesia condition (such as right after the bolus injection), neuronal activity and corresponding fMRI responses are suppressed (see Shim et al., 2018). In such an intermittent anesthesia protocol, the observation of no BOLD response was regarded as an indication of deep anesthesia, and these runs are often appropriately excluded for further data analysis. Throwing away bad trials is often practiced, as these trials presumably occur due to poor anesthetized animal conditions that lead to a null, independently measured BOLD response. However, it would be incorrect to use temporal similarity to the hypothesized DIANA response as a metric to select “good” from “bad” trials, as spurious noise could then average together to produce the hypothesized result. To demonstrate the impact of a biased exclusion of “presumably bad trials” on DIANA studies, we re-processed our DIANA data with 50 – 300 trials (Fig. 3). Initially, the correlation between each time series and a negative response function (Fig. 4B) was calculated for identifying presumably outlier trials systematically. Trials were ranked based on their cross-correlation values with this function, and the top 6%, 10%, and 20% ranks were used as an outlier exclusion threshold. Then, the averaged time course of the included (94%, 90%, and 80% of all trials) or excluded trials was obtained for each subject in each exclusion threshold. The original data (Fig. 4A and Supplementary Fig. 1A for individual animals; Fig. 4C and Supplementary Fig. 1C for an average of 6 animals) show no neuronal activity-like peak. When some trials were excluded, the averaged time course of included trials for each individual animal showed noisy positive responses (Fig. 4D for active ROI and Supplementary Fig. 1D for inactive ROI), while the averaged time course of excluded trials showed a negative peak at 25 ms (Fig. 4E for active ROI and Supplementary Fig. 1E for inactive ROI). When only 10% of total trials (e.g., 5 out of 50) were excluded, noisy time courses having an observable peak were observed (Fig. 4F for active ROI and Supplementary Fig. 1F for inactive ROI). Peak intensity was also closely dependent on the exclusion threshold and tSNR; a larger spurious peak intensity was observed at a higher exclusion threshold (Fig. 4F and Supplementary Fig. 1F), not surprisingly as one is effectively rejecting more outlier responses and selectively averaging positive noise peaks that occur when the hypothesized signal is occurring. When data exclusion is performed based on any pre-selection model of the expected response implicitly (visual inspection) or explicitly (cross-correlation analysis), then statistical circularity occurs, invalidating results, as seen in Fig. 4 and Supplementary Fig. 1. 3. Discussion In our high temporal resolution fMRI studies with 50–300 averages at 15.2 T, ROI-wise tSNR was ~2,000 to ~4,000 in the somatosensory barrel cortex, values that are very high by normal BOLD fMRI standards. A tSNR of 3000 to 1 corresponds to a standard deviation of 0.033%, and thus our data have more than sufficient statistical power to detect the reported ~0.15% DIANA peak in individual animals. However, we did not find any direct neuronal-related signal in our extensively averaged fMRI studies (Fig. 3), suggesting that the lack of observation of the DIANA peak cannot be due to insufficient tSNR. Since BOLD responses were reliably observed in all the animals, poor control of anesthetic depth is also not a cause of no DIANA observations. Our data suggest that the reported DIANA peak in anesthetized mice (Toi et al., 2022) is not likely neuronal-related fMRI peak, but artifacts due to the noise characteristics of the insufficiently averaged data. This is apparent in Fig. 2B, where spatially distributed random voxels can be found that match the cross-correlation template when using only 50 trials. This random distribution shows apparent activation in S1BF, the thalamus and many other regions of the brain. This random pattern is in contrast to the BOLD activation maps in Fig. 1C and Fig. 1F. When time points and averages (i.e., sample size) are limited, then it is possible that a noise peak can be mistakenly identified as a genuine peak as seen in Fig. 2. This potential problem increases with decreasing tSNR. To separate between genuine peaks and artifacts, extensive averaging should be performed until a peak amplitude is much greater than baseline signal fluctuations (see Fig. 3). Alternatively, if each run has sufficient temporal frames with repeated stimuli, then we can evaluate whether the assumed peak is reproducible across different stimulation pulses. Thus, one should be extremely careful when trying to identify and interpret an unknown peak from limited sampled data. Our data indicate that a spurious, neuronal activity-like peak can be generated from noisy data with limited time points and trials when some trials are excluded by comparison to the hypothesized signal change. We were able to reproduce DIANA signal characteristics observed in Toi et al. (2022) by the exclusion of outlier-like trials from noisy time courses. For most practical fMRI studies, a few tens of repetitions are all that are acquired due to the limited experimental time in vivo. Here we show that the exclusion of even a few trials based on selective filtering may produce artifactually positive findings which are not seen when no outliers are rejected (Fig. 4). This demonstrates that positive results can be produced from selective exclusion of spurious time series. Extreme care should be taken to ensure that trial exclusion is based on unbiased metrics such as separately measured BOLD response magnitude, so that circularity is avoided. Objective, statistically justifiable inclusion/exclusion criteria, having no bias based on any temporal features of the hypothesized result, should be used for data screening and for ensuring that overall findings are preserved, regardless of the exclusion criteria. While the biased metric approach was the only way that we were able to reproduce the findings of Toi et al. (2002), it should be noted that they did not report using this approach in their paper. Therefore, it remains unknown as to why the results were not replicated even with more extensive averaging, a more precise and steady anesthetic protocol, and higher field strength. Other possibilities include subtle differences in pulse sequence (i.e., RF spoiling approach) or stimulus quality. 4. Methods 4.1. Animal preparation and stimulation Six adult male C57BL/6 mice (23–27 g, 10–12 weeks old; Orient Bio, Korea) were used with approval by the Institutional Animal Care and Use Committee (IACUC) of Sungkyunkwan University in accordance with the standards for humane animal care. All MRI experiments were performed under anesthesia in accordance with the guidelines of the Animal Welfare Act and the National Institutes of Health Guide for the Care and Use of Laboratory Animals. For anesthesia, a mixture of ketamine and xylazine (100 mg/kg and 10 mg/kg, respectively) was initially injected intraperitoneally (IP), and a dose of ketamine 45 mg/kg/h and xylazine 2.25 mg/kg/h was continuously infused intravenously (IV) (You et al., 2021). Note that Toi et al. (2022) used the supplementary bolus IP delivery of 25 mg/kg ketamine and 1.25 mg/kg xylazine when needed (Shim et al., 2018; Jung et al., 2019; Shim et al., 2020; Dinh et al., 2021; Jung et al., 2021), which induces anesthetic depth deep right after the bolus injection and a slow recovery to wakefulness (Jaber et al., 2014; Shim et al., 2018). The animals were self-breathing under continuous supply of oxygen and air gases (1:4 ratio) through a nose cone at a rate of 1 liter/min (Shim et al., 2020). To reduce head motions during image scanning, the mouse was carefully positioned on a customized cradle with two ear plugs, a bite bar and nose cone. Body temperature was maintained at 37 ± 0.5°C with a warm-water heating system via rectal temperature probe. For whisker electric pad stimulation, two anodes with 2 mm apart and a cathode center (You et al., 2021) were placed on the mouse’s right whisker-pad. Pulse parameters were at a pulse width of 0.5 ms, a frequency of 4 Hz (for BOLD studies), and current intensity of 0.5 mA, controlled by a pulse generator (Master 9; World Precision Instruments, Sarasota, FL, USA). 4.2. fMRI data collection Data were acquired on a 15.2 T (Bruker BioSpec MRI, Billerica, MA, USA) equipped with a 11 cm horizontal bore magnet and actively shielded 6-cm gradient. A 15 mm ID surface coil was used for both transmission and reception. Note that Toi et al used a 10 mm diameter surface coil at 9.4 T. Mouse brain was placed at the isocenter of magnet and field inhomogeneity was minimized via MAPSHIM protocol in Paravision 6.0.1 software. Detailed experimental procedures were described in our previous publications (Jung et al., 2019; Jung et al., 2021). Three fMRI approaches were used as follows. 1) BOLD-EPI: Standard gradient-echo echo planar imaging (EPI) was used to obtain BOLD fMRI with the following imaging parameters; image matrix = 96 × 48, FOV = 16 × 12 mm2 (0.17 × 0.25 mm2 in-plane resolution), 1.0 mm-thick slices, TR/TE = 1000/11.5 ms, bandwidth = 300,000 Hz, and 50° flip angle. Then, a single 1-mm slice was chosen based on scout BOLD fMRI (approximately −1.75 mm from bregma) for subsequent fMRI studies. Single-slice BOLD fMRI was performed for ensuring reliable BOLD activity. The stimulation paradigm consisted of 40 s control, 20 s stimulation, 60 s control, 20 s stimulation, and 60 s recovery, lasting 200 s. Two runs of BOLD fMRI were repeated before and after DIANA scans to ensure their stable responsiveness during entire fMRI studies. 2) BOLD line scanning (BOLD-LS): The 2-D shuffled line scanning originally proposed by Silva & Koretsky (2002) was adopted for single-slice BOLD fMRI with TR/TE = 100/11.7 ms, FOV =16 × 12 mm2, image matrix =72 × 54, bandwidth = 50,000 Hz, slice thickness = 1 mm, and flip angle = 17° in the S1BF (Ernst flip angle with T1 of 2.2 s). The stimulation paradigm consisted of 2 s control, 1 s stimulation, 6s control, 1 s stimulation, and 6 s recovery. Each trial lasted 16 s × 54 phase encoding steps = 864 s. 3) DIANA: DIANA experiments were performed with the 2-D shuffled line scanning approach with the same imaging parameters used in Toi et al. (2022), TR/TE = 5/2 ms, 4° flip angle in the S1BF, FOV = 16 × 12 mm2, image matrix =72 × 54, slice thickness = 1 mm, and bandwidth = 50,000 Hz. In some studies, RF spoiling and steady-state dummy scans were used (see Table 1). Two different numbers of frames for each phase-encoding step were acquired, 40 (200 ms) or 200 frames (1000 ms). For DIANA studies with 40 time frames, 50 ms pre-stimulation (10 frames) and 150 ms post-stimulation images (30 frames) were acquired, lasting 200 ms × 54 phase encoding steps = 10.8 s/trial, which is the same as Toi et al. (2022) used. Alternatively, 20 pre-stimulus and 180 post-stimulus images, or 20 pre-stimulus and 2 times 90 post-stimulus images were obtained (see Table 1). At least 50 trials were repeated. 4.3. fMRI data processing All the acquired data were analyzed with AFNI software (https://afni.nimh.nih.gov/, Cox, 1996), FMRIB Software (FSL, https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/, Smith et al., 2004), and Matlab codes (Mathworks). Generating fMRI maps: Each individual animal’s fMRI maps were generated using preprocessing and a cross-correlation analysis with a hemodynamic response function (HRF) (Jung et al., 2019) or an expected DIANA response function (Toi et al., 2022). The DIANA response function was modeled as a Gaussian window function in Matlab, ‘gausswin’, with a peak position at 25 ms after the stimulus and the full width at half maximum (FWHM) of ~25 ms. The following preprocessing steps were performed to improve the detection of signal activation: linear detrending for signal drift removal and normalization by the average of the pre-stimulus baseline volumes. All repeated fMRI trials on each animal were averaged and spatial smoothing was applied using a Gaussian kernel with FWHM of one voxel. ROI analysis: Active and inactive ROI with 5×5 voxels was chosen in the contralateral S1BF and the ipsilateral thalamic area from BOLD-EPI functional maps, respectively. Since EPI and line scanning used a slightly different spatial resolution, the ROIs determined from BOLD-EPI were similarly positioned in the line scanning images. In each run (BOLD-EPI, BOLD-LS, DIANA), signals in the selected 5×5-voxel ROI were averaged from detrended, normalized image frames (without any spatial or temporal smoothing), then were temporally smoothed with 3-point moving averages, which was used in Toi et al. (2022). In each subject, repeated DIANA runs with 10 pre-stimulus and 30 post-stimulus frames were temporally aligned, based on the stimulus time. The 10th image (right before the stimulus pulse) was set to 0 ms for being consistency with Toi et al. (2022). Since TE of 2 ms was used, the exact time of the 11th image (the image right after the stimulus pulse) was 2 ms, but assigned to 5 ms. The percent signal change time course was calculated by 100×ΔS/Sbase, where ΔS is a difference between signals at the reference and at the baseline period and Sbase is the averaged baseline pre-stimulus signal (1st to 10th frame). In each subject, 50–300 trials were averaged. Finally, averaged subject-wise time courses were averaged. Since different animals had different number of averages, a grand total average was also performed from 1050 trials of the six subjects. BOLD percent change and tSNR calculation: In each subject, ROI-wise BOLD-EPI responses to whisker stimulation were calculated from 0 – 35 s time frames after the stimulus onset, while ROI-wise BOLD-LS changes were calculated from 1 – 4.5 s time frames after the stimulus onset. To determine the detectability of DIANA measurements, tSNR values were calculated by standard deviations divided by its means from time series data of individual voxels or ROIs. Although tSNR calculation generally uses only pre-stimulus data points, we used the entire time series due to no obvious observation of stimulation-induced responses. Histogram analysis: Histograms of individual DIANA time series were obtained by counting the number of data points every 0.1 bins. A Gaussian fitting to histogram was performed and its mean and standard deviation were determined. Selection of trials: To evaluate the impact of trial exclusion, each trial-wise ROI time course with 3-point moving averaging was cross-correlated with a hypothetical outlier Gaussian function with a negative peak time of 25 ms and FWHM of ~25 ms. Then, trials were ranked, based on cross-correlation values. Various thresholds of top 6%, 10%, and 20% were used for the separation between included and excluded trials. Then, the included and excluded trials were separately averaged for each animal. Statistical analysis: All statistical analyses between the two BOLD groups (Fig. 1) were conducted with paired t-test. To determine a statistical significance in Fig. 2, all pre-stimulus baseline data (20 time point × 50 trials) were compared with peak data points (1 time point × 50 trials) with two-group t-test. Quantitative values were presented as mean ± standard deviation and plots were presented as mean ± standard error of the mean (SEM). Supplementary Material Supplement 1 Acknowledgments This research was supported by the Institute of Basic Science (IBS-R015-D1 to S.-G.K.), and by NIH (RF1NS113278, R01NS124778, R01NS122904 to X.Y and S.C). Figure 1. BOLD fMRI of somatosensory stimulation in anesthetized mice at 15.2 T. (A-B) Experimental schematics and stimulation paradigms of BOLD fMRI with EPI and line scanning (LS). (C) BOLD fMRI map overlaid on original EPI image of one animal (Mouse #1). Black box: active 5×5-voxel S1BF ROI; red box: inactive 5×5-voxel ROI. (D) BOLD time courses of the active S1BF ROI before and after DIANA experiments in Fig. 1C. Vertical shaded bars: stimulus durations. (E) Average percent changes of the S1BF ROI before and after DIANA experiments (n = 6 mice). Each data point: each mouse. (F) Conventional BOLD fMRI map obtained from only one shuffled line scanning fMRI run. (G) Time courses of active and inactive ROIs in Fig. 1F. Vertical shaded bars: stimulus durations. (H) Averaged percent changes of active and inactive ROIs (n = 5). Each data point: each mouse; N.A: non active; error bars: SEM. Figure 2. Systematic analysis of 50 DIANA trials in one animal. One DIANA trial composed of total 20 pre-stimulus and 180 post-stimulus frames (i.e., 1000-ms inter-stimulus interval) (Mouse #3, see Table 1). (A) tSNR map calculated from the average of 50 trials. (B) Cross-correlation map obtained with an expected neuronal response curve from Toi et al. (Toi et al., 2022). (C) Trial-wise time courses of the active S1BF ROI. To visualize individual trials, 5 trials were plotted per row. Vertical bar: stimulus. (D) The averaged time course of 50 repeated trials with an expanded view in inset (red). An expected DIANA response was also plotted (blue). A statistically significant positive peak was detected at 10 ms after the stimulus. Shaded area: SEM. (E) A histogram of Figure 2C data points (50 trials × 200 time points). A Gaussian noise distribution was observed with a standard deviation of 0.377%. (F) A histogram of 50 individual trials at a time (10 ms post-stimulus) of the statistically significant positive peak. Figure 3. No observable DIANA responses in all six subjects. (A-D) Individual (A, C) and averaged time courses (B, D) from the active contralateral S1BF (A-B) and inactive ipsilateral ROIs (C-D) in each individual animal (row). ROIs were shown with the total number of trials in each animal (images between B and C). 50–300 trials were repeated in each animal. (E-H) Averages of all 1050 trials in 6 animals (E, G) and 6 averaged animal’s time courses (F, H) from the active S1BF (E, F) and inactive ROI (G, H). No noticeable neuronal activity-like peak was observed. Vertical bar: stimulus; shaded area: SEM. Figure 4. Effect of data exclusion on noisy time courses. (A) Averaged time courses of six individual animals in the active S1BF ROI (same as Fig. 3B). (B) A Gaussian reference outlier function with a full width of half minimum of ~25 ms peaked at 25 ms. (C) An average of all animals’ time courses without excluding trials. This is the same as Fig. 3F. In six animals, an individual trial’s time course was correlated with a reference function shown in B, and ranked among all trials, based on its cross-correlation value. Then, in each animal, trials were separated into included and excluded categories, based on an outlier threshold of top 6%, 10%, and 20% cross-correlation values. (D-E) Averaged time courses of the included (D) and excluded trials (E) in individual animals. Averaged time courses of the included trials showed noisy positive responses, while those of the excluded trials showed a negative peak around the expected peak time. Each color time course: each animal. (F) Means of subject-wise included trials time courses. Clearly, an exclusion process leads to an erroneous peak from noisy data. It is fundamentally important to note that these peaks are erroneous because such a preselection process is statistically circular, summing noise in a biased manner to produce precisely the peak that was preselected. Error bars: SEM. Table 1. Detailed experimental information Mouse ID Paradigm1 (ms) ISI2 (s) Pulse # Average Trial #3 RF-spoil4 Steady State 1 50–150 0.2 1 50 50 x x 2 50–150 0.2 1 100 100 x x 3 50–150 0.2 1 50 100 x x 100–900 1.0 1 50 x x 4 50–150 0.2 1 100 200 x x 100–450–450 1.0 2 50 x x 5 50–150 0.2 1 100 300 x o 50–150 0.2 1 100 o o 100–450–450 1.0 2 50 o x 6 50–150 0.2 1 100 300 x o 50–150 0.2 1 100 o o 100–450–450 1.0 2 50 o x Total 50–150 1 1050 1 Stimulation paradigm. “−” indicates a 0.5-ms single whisker pad electric stimulus, and image TR is 5 ms. 2 Inter-scan interval. 3 Number of trials time-locked with a stimulus pulse. 4 Radio frequency (RF) spoiling. ==== Refs References Bandettini PA , Wong EC , Hinks RS , Tikofsky RS , Hyde JS . Time course EPI of human brain function during task activation. Magn Reson Med. 1992 Jun;25 (2 ):390–7. doi: 10.1002/mrm.1910250220.1614324 Bellgowan PS , Saad ZS , Bandettini PA . Understanding neural system dynamics through task modulation and measurement of functional MRI amplitude, latency, and width. Proc Natl Acad Sci U S A. 2003 Feb 4;100 (3 ):1415–9. doi: 10.1073/pnas.0337747100. Epub 2003 Jan 27. 12552093 Cox RW . AFNI: software for analysis and visualization of functional magnetic resonance neuroimages. Comput Biomed Res. 1996 Jun;29 (3 ):162–73. doi: 10.1006/cbmr.1996.0014.8812068 Dinh TNA , Jung WB , Shim HJ , Kim SG . Characteristics of fMRI responses to visual stimulation in anesthetized vs. awake mice. Neuroimage. 2021 Feb 1;226 :117542. doi: 10.1016/j.neuroimage.2020.117542. Epub 2020 Nov 10. 33186719 Hodono S , Rideaux R , van Kerkoerle T , Cloos MA . Initial experiences with Direct Imaging of Neuronal Activity (DIANA) in humans. 2023. https://arxiv.org/abs/2303.00161 Jaber SM , Hankenson FC , Heng K , McKinstry-Wu A , Kelz MB , Marx JO . Dose regimens, variability, and complications associated with using repeat-bolus dosing to extend a surgical plane of anesthesia in laboratory mice. J Am Assoc Lab Anim Sci. 2014 Nov;53 (6 ):684–91.25650976 Mouse BOLD fMRI at ultrahigh field detects somatosensory networks including thalamic nuclei. Jung WB , Shim HJ , Kim SG . Neuroimage. 2019 Jul 15;195 :203–214. doi: 10.1016/j.neuroimage.2019.03.063. Epub 2019 Apr 1. 30946950 Jung WB , Im GH , Jiang H , Kim SG . Early fMRI responses to somatosensory and optogenetic stimulation reflect neuronal information flow. Proc Natl Acad Sci U S A. 2021 Mar 16;118 (11 ):e2023265118. doi: 10.1073/pnas.2023265118.33836602 Kwong KK , Belliveau JW , Chesler DA , Goldberg IE , Weisskoff RM , Poncelet BP , Kennedy DN , Hoppel BE , Cohen MS , Turner R , Dynamic magnetic resonance imaging of human brain activity during primary sensory stimulation. Proc Natl Acad Sci U S A. 1992 Jun 15;89 (12 ):5675–9. doi: 10.1073/pnas.89.12.5675.1608978 Menon RS , Kim SG . Spatial and temporal limits in cognitive neuroimaging with fMRI. Trends Cogn Sci. 1999 Jun;3 (6 ):207–216. doi: 10.1016/s1364-6613(99)01329-7.10354573 Ogawa S , Tank DW , Menon R , Ellermann JM , Kim SG , Merkle H , Ugurbil K . Intrinsic signal changes accompanying sensory stimulation: functional brain mapping with magnetic resonance imaging. Proc Natl Acad Sci U S A. 1992 Jul 1;89 (13 ):5951–5. doi: 10.1073/pnas.89.13.5951.1631079 Pohmann R , Speck O , Scheffler K . Signal-to-noise ratio and MR tissue parameters in human brain imaging at 3, 7, and 9.4 tesla using current receive coil arrays. Magn Reson Med. 2016 Feb;75 (2 ):801–9. doi: 10.1002/mrm.25677. Epub 2015 Mar 29. 25820458 Shim HJ , Jung WB , Schlegel F , Lee J , Kim S , Lee J , Kim SG . Mouse fMRI under ketamine and xylazine anesthesia: Robust contralateral somatosensory cortex activation in response to forepaw stimulation. Neuroimage. 2018 Aug 15;177 :30–44. doi: 10.1016/j.neuroimage.2018.04.062. Epub 2018 May 4. 29730495 Shim HJ , Lee J , Kim SG . BOLD fMRI and hemodynamic responses to somatosensory stimulation in anesthetized mice: spontaneous breathing vs. mechanical ventilation. NMR Biomed. 2020 Jul;33 (7 ):e4311. doi: 10.1002/nbm.4311. Epub 2020 Apr 15. 32297409 Shim HJ , Im GH , Jung WB , Moon HS , Dinh TNA , Lee JY , Kim SG . Protocol for mouse optogenetic fMRI at ultrahigh magnetic fields. STAR Protoc. 2022 Dec 16;3 (4 ):101846. doi: 10.1016/j.xpro.2022.101846. Epub 2022 Dec 12. 36595930 Silva AC , Koretsky AP . Laminar specificity of functional MRI onset times during somatosensory stimulation in rat. Proc Natl Acad Sci U S A. 2002 Nov 12;99 (23 ):15182–7. doi: 10.1073/pnas.222561899. Epub 2002 Oct 29. 12407177 Smith SM , Jenkinson M , Woolrich MW , Beckmann CF , Behrens TE , Johansen-Berg H , Bannister PR , De Luca M , Drobnjak I , Flitney DE , Niazy RK , Saunders J , Vickers J , Zhang Y , De Stefano N , Brady JM , Matthews PM . Advances in functional and structural MR image analysis and implementation as FSL. Neuroimage. 2004;23 Suppl 1 :S208–19. doi: 10.1016/j.neuroimage.2004.07.051.15501092 Toi PT , Jang HJ , Min K , Kim SP , Lee SK , Lee J , Kwag J , Park JY . In vivo direct imaging of neuronal activity at high temporospatial resolution. Science. 2022 Oct 14;378 (6616 ):160–168. doi: 10.1126/science.abh4340. Epub 2022 Oct 13. 36227975 Uludağ K , Blinder P . Linking brain vascular physiology to hemodynamic response in ultra-high field MRI. Neuroimage. 2018 Mar;168 :279–295. doi: 10.1016/j.neuroimage.2017.02.063. Epub 2017 Feb 22. 28254456 You T , Im GH , Kim SG . Characterization of brain-wide somatosensory BOLD fMRI in mice under dexmedetomidine/isoflurane and ketamine/xylazine. Sci Rep. 2021 Jun 23;11 (1 ):13110. doi: 10.1038/s41598-021-92582-5.34162952 Yu X , Qian C , Chen DY , Dodd SJ , Koretsky AP . Deciphering laminar-specific neural inputs with line-scanning fMRI. Nat Methods. 2014 Jan;11 (1 ):55–8. doi: 10.1038/nmeth.2730. Epub 2013 Nov 17. 24240320