Robust high spatio-temporal line-scanning fMRI in humans at 7T using multi-echo readouts, denoising and prospective motion correction

Background: Functional magnetic resonance imaging (fMRI), typically using blood oxygenation level-dependent (BOLD) contrast weighted imaging, allows the study of brain function with millimeter spatial resolution and temporal resolution of one to a few seconds. At a mesoscopic scale, neurons in the human brain are spatially organized in structures with dimensions of hundreds of micrometers, while they communicate at the millisecond timescale. For this reason, it is important to develop an fMRI method with simultaneous high spatial and temporal resolution.


Introduction
In the human brain, neurons with similar functions cluster together in spatial structures with an extent of hundreds of micrometers, i.e. mesoscopic scale, while they communicate at the millisecond timescale (Dumoulin et al., 2018;Sabatini and Regehr, 1996). For example, in the cortex neurons are organized in columnar and layer structures (Brodmann, 1909;Hubel and Wiesel, 1968;Mountcastle, 1957). Cortical layers differ in neuronal content, but also in their connectivity to other parts of the brain. Brodmann (Brodmann, 1909) used differences in cortical layers to define distinct cortical areas. Thus, cortical layers are considered a basic anatomical and physiological unit of the cortex and for this reason it is interesting to study responses at the mesoscopic scale.
Functional magnetic resonance imaging (fMRI) plays an important role in the study of brain function (Ogawa et al., 1990), allowing detection of brain activity through the changes in blood flow and oxygenation. Blood oxygenation level-dependent (BOLD) contrast weighted imaging is one of the primarily used contrast mechanisms in both cognitive (Logothetis, 2008) and clinical (Eickhoff et al., 2020) neuroscience. BOLD imaging has two important advantages: it is non-invasive and readily available at any MRI scanner. For BOLD fMRI, high spatial and temporal (<1 mm, ~100 ms) resolutions are required to detect spatiotemporal features of the hemodynamic response function (HRF) which describes how the hemodynamic signal propagates through the cortex at the mesoscopic scale, in particular across cortical layers.
However, functional MRI is an inherently noisy acquisition method (Liu, 2016). As a result, most fMRI data are denoised in some form before statistical testing, usually by applying a spatial smoothing step (Ashburner, 2012;Jenkinson et al., 2012;Smith, 2004). For high-resolution acquisitions, smoothing is not appropriate because of the concurrent loss in spatial definition (Stelzer et al., 2014), even though noise levels increase at these higher resolutions. This means that other techniques have to be applied to increase the signal-to-noise-ratio (SNR) and decrease the noise (Caballero-Gaudes and Reynolds, 2017).
Advances in fMRI methodology have been aimed at increasing both the spatial and the temporal resolution, with the ultimate goal of recording at sub-millimeter spatial resolution and sub-second sampling rate. In the past decade, line-scanning fMRI in rodents achieved very high resolution across cortical depth (50 µm) and time (50 ms) by scanning only a single line of data, sacrificing volume coverage and resolution along the cortical surface in the process (Yu et al., 2014). Line-scanning fMRI in rodents has also been used in combination with other techniques such as fiber-based optogenetic stimulation (Albers et al., 2018) and diffusion-sensitizing gradients (Nunes et al., 2021), to investigate the fast BOLD onset times at high spatial resolution. Line-scanning fMRI in rats has been extended towards a multi-line implementation, acquiring line profiles from different cortical regions to investigate laminar-specific functional connectivity mapping under both evoked and resting-state conditions (Choi et al., 2021).
Recently, line-scanning was employed in human studies, for microstructural investigations (Balasubramanian et al., 2022(Balasubramanian et al., , 2021 and with the goal of isolating microvessel responses and characterizing the distribution of blood flow and laminar functional MRI profiles across cortical depth, at high spatiotemporal resolution (Morgan et al., 2020;Raimondo et al., 2021). Our first human line-scanning implementation (Raimondo et al., 2021) achieved resolutions of 250 µm and 200 ms and demonstrated its utility for measuring BOLD responses along cortical depth, in the visual cortex, during a visual stimulation task. The sequence was a modified 2D gradient-echo (GE) sequence employing a single-echo readout, with low bandwidth for increased SNR, but still with some deadtime within the repetition time (TR). Here, we propose several ways to improve this line-scan sequence for increased sensitivity. This increased sensitivity is necessary to make line-scanning a more robust and generalizable technique, which could be used for future neuroscientific and clinical applications. Specifically, diseases such as small vessel or sickle cell diseases (Afzali-Hashemi et al., 2021;DeBaun and Kirkham, 2016;van den Brink et al., 2022;Zwanenburg and van Osch, 2017) would benefit from high spatial and temporal resolutions to gain insights about the altered hemodynamic responses in patients across cortical depths.
First, the unused deadtime within the TR permits the acquisition of additional echo readouts without increasing the TR. The BOLD contrast is known to be maximal when the echo time (TE) is equal to the local tissue T 2 * relaxation rate (Olsrud et al., 2008). With multi-echo imaging, we can measure the T 2 * signal decay curve and reach optimal readout efficiency and BOLD sensitivity. In addition, the T 2 * signal decay curve can be used to disentangle BOLD-like (T 2 *) changes from non-BOLD signal changes (Kundu et al., 2012). Non-BOLD signals can be caused by drift, motion, physiological noise or other contaminating signals that impact the initial signal intensity (S 0 ) of the T 2 * decay curve sensitivity (Caballero-Gaudes and Reynolds, 2017;Poser et al., 2006;Posse et al., 1999).
Second, thermal noise dominates at the high spatial and temporal resolution of line-scanning. The "MR signal" can be specifically defined as an electrical current induced in the receiver coil by the precession of the net magnetization during resonance, as the manifestation of Faraday's Law of Induction (Hahn, 1953). However, the definition of noise in an fMRI time series is more complex, due to the different noise sources (Krüger and Glover, 2001). Thermal noise, classified as a zero-mean Gaussian distributed noise, is generated either from the electronics or from the sample being imaged and depends on a range of parameters including the static magnetic field strength, the electronics, the readout bandwidth and sampling scheme (Edelstein et al., 1986;Hoult and Richards, 2011). It becomes predominant when small voxel sizes are used. Other sources of signal variance include subject motion and physiological noise through respiration and heartbeat (Triantafyllou et al., 2005). For denoising, the final goal is to decrease noise without compromising any physiological aspects of the data, but, generally, some information has to be sacrificed. Many denoising techniques (Alkinani and El-Sakka, 2017;Fan et al., 2019;Kaur et al., 2018) are based on a trade-off between the removal of unwanted signal and the preservation of the data quality, such as spatial and temporal resolution, as well as the underlying biological processes. In a newly described approach from (Vizioli et al., 2021), thermal noise is selectively suppressed from high-resolution fMRI data while preserving the amplitude of the hemodynamic response, the spatial resolution and the functional point-spread function.
Finally, line-scanning is highly sensitive to subjects' motion. Generally in fMRI, participant movement leads to inconsistencies in the fMRI timecourse. These are typically corrected by coregistering the volume timepoints (Friston et al., 1994;Smith et al., 2004). Motion is even more problematic in high-resolution MRI, where the impact of smaller movements is amplified (Gallichan et al., 2016;Maclaren et al., 2012;Schulz et al., 2012). In the case of line-scanning, the 1-dimensional nature of the data only allows motion detection in the line direction, whereas rotations or displacements perpendicular to the line lead to spin-history artifacts or line acquisition outside the area of interest. These effects cannot be corrected by post-hoc motion correction in line scans; therefore, a prospective motion correction scheme (PMC) is required . This can be achieved with external hardware Schulz et al., 2012;Stucht et al., 2015) or with image-based navigators (Andersen et al., 2019;Glover et al., 2000). Here we use an image-based navigator implementation.
In this work, we aim to improve the line-scanning acquisition in three ways: first, we investigated five different multi-echo protocols varying the numbers of echoes and readout bandwidth while keeping the TRand thus the overall scan timeconstant. We compared different echo combination approaches in terms of BOLD activation (sensitivity) and temporal signal-to-noise ratio. Second, we implemented an adaptation of the NORDIC thermal noise removal for linescanning fMRI data, and finally, we tested three image-based navigators for motion correction.

Participants
Our pool of participants is composed of two groups: multi-echo section: 6 participants (3 male, 32 ± 6 years old (mean ± standard deviation)), noise removal section: 4 participants (4 male, 28 ± 4 years old) and for the motion correction section: 8 participants dataset (7 male, 31 ± 8 years old, including the individuals participating in the noise removal section). All the participants were healthy individuals who provided written informed consent as approved by the medical ethics committee of the Amsterdam University Medical Centre. The guidelines of the Helsinki Declaration were followed throughout the study, and all participants were screened for MRI compatibility prior to the experiments.
Note that the line-scanning method allows a high sampling rate within the participant, hence we can reach statistical power without the need to average over many participants. In addition, we aimed to find an effect in every single participant, which is easily achieved with the strong visual task we employed across the study. The number of participants can be seen as replication of the same effect rather than a way to measure the effect itself (Normand, 2016;Smith and Little, 2018). Ultimately, for clinical research it is essential to maximize signal and minimize noise to have information in individuals (Gratton et al., 2022).

Selection of multi-echo acquisition and echo combination version
The 7 T MRI (Philips Healthcare, Best, NL) was equipped with a 2 channel transmit and 32 channel local receive surface coil (Petridou et al., 2013), positioned close to the visual cortex.
Line-scanning acquisition used a modified 2D multi-echo gradientecho sequence where the phase-encoding in the direction perpendicular to the line, needed for conventional 2D imaging, was omitted: line res-olution= 250 µm, TR= 105 ms (108 ms for one participant), flip angle= 16 • , array size= 720, line thickness= 2.5 mm, in-plane line width= 4 mm, fat suppression using SPIR. Two saturation pulses (7.76 ms pulse duration) suppressed the signal outside the line of interest. Five different multi-echo schemes (including 3, 5, 7, 9 and 11 echoes) were compared by adapting the readout bandwidth for the different schemes. The longest echo time for all schemes was 38 ms. The order of acquisitions was randomized across participants; details are provided in Table 1.
The line was manually positioned approximately perpendicularly to the medial gray matter sheet of the occipital lobe. We acquired one run of functional data with each protocol, using a block design visual task consisting of an 8 Hz flickering checkerboard presented on the entire screen for 10 s on/off. Runs lasted 5 min and 40 s. The total run duration for the 11 echoes acquisition runs was shortened for 3 participants due to technical constraints and skipped for one other participant. Reconstruction was performed offline (MatLab, Gyrotools).
We combined the multi-channel coil data with a temporal signal-tonoise ratio (tSNR) and coil sensitivity-weighted sum of squares (SoS) weighted scheme per echo in the data reconstruction as in (Raimondo et al., 2021). The SoS was defined as: where S i is the MR signal for each receive coil channel and N is the number of channels. The resulting S(t, TE) is the MR signal as a function of time and TE. Multi-echo data were combined in 3 different ways: SoS (with the same formula as Eq. (1), but summing over echoes instead of channels), T 2 * fit and a tSNR weighted combination (wtSNR) based on Poser et al. (Poser et al., 2006). 'T 2 * fit' fits the T 2 * decay per voxel with a linear fit of Eq. (2) in a least-squares way: Functional data were analyzed in Matlab using a general linear model (GLM) approach. We used the SPM implementation of the canonical HRF (Friston et al., 1994) as block design experiments do not allow one to fit the HRF shape. T-statistic values (t-stats) were computed to identify active voxels. We also evaluated the tSNR for each voxel through: where S(t) is the mean signal over time, across the whole timecourse and σ(S(t)) is the standard deviation of the signal across time for the whole timecourse. Note that the tSNR is computed across the whole timecourse, hence including voxels with functional activation due to the visual task. However, within the line, very few voxels contain task activation and task effects on the tSNR are minimal. We compared the mean and maximum t-stats in a region of interest (ROI) for the 5 different multi-echo acquisitions, the three echocombination methods and the ROI average tSNR. The ROIs were defined as the 11 voxels, 11 × 0.25 mm = 2.75 mm, covering the gray matter (GM) ribbon (identified in a slice image centered on the line), surrounding the voxel showing the highest t-stats, for all acquisitions. For the tSNR comparison, we also defined a white matter ROI (WM ROI).

Noise removal
For the noise removal dataset, we acquired data using the 5 echoes acquisition with the same visual task as the previous section, adding 20 s baseline in the beginning.
In fMRI, it is common to perform some kind of noise reduction data processing to increase the SNR, with the final goal of maintaining signal integrity (including spatial and temporal resolutions). To achieve this purpose, for this section, we employed a thermal noise removal step in the reconstruction pipeline based on Noise reduction with Distribution Corrected (NORDIC) principal component analysis (PCA) (Vizioli et al., 2021). The denoising was applied to k-space data before the coil and the echo combination (see Fig. 1). A singular value decomposition (SVD) of the data was submitted to "hard thresholding" that eliminates all components indistinguishable from zero-mean Gaussian distributed noise (Fig. 1a). The singular value decomposition of our line-scanning k-space data matrix (for every channel and echo) was U ⋅S ⋅V H , where U and V are unitary matrices, and S is a diagonal matrix whose diagonals are the spectrum of ordered singular values. The singular values below a chosen threshold were replaced by 0, and the other singular values were unaffected. S th is a new diagonal matrix generated as a result of thresholding, and the estimate of the NORDIC-denoised data was given as U ⋅S th ⋅V H . The threshold was chosen from the elbow point of the 'scree plot' depicting the eigenvalues versus the number of components (see Fig. 1) for every channel and every echo separately. Note that, during a pilot study we ensured that noise was preferentially removed, leaving Table 1 Parameters for the 5 multi-echo line-scan acquisitions. N-echoes is the number of echoes, TE1 is the first acquired echo, echo spacing is the time difference between echoes acquisitions and readout BW is the readout bandwidth per pixel. The last echo time in the series was always 38 ms. task-driven signal variation, by verifying that t-values rose as a function of removed variance. Thermal noise removal was then followed by a tSNR and coil sensitivity-weighted SoS multi-channel data combination and SoS echo combination. To estimate the performance of the reconstructions with and without denoising, we compared the two in terms of tSNR, and tstats in an 11 voxel ROI centered around the voxel with maximum t-stats within the brain.  (Vizioli et al., 2021). Line-scanning k-space data for every channel and echo separately were decomposed through singular value decomposition (SVD); the diagonal matrices containing the eigenvalues (S) were thresholded through the elbow (red arrow) of the scree plot depicting the eigenvalues versus the components.

Prospective motion correction
For the motion correction section, we combined line-scanning with prospective motion correction (PMC) using an interleaved scanning architecture (MISS, Philips) following (Andersen et al., 2019). We inserted a navigator at every dynamic (i.e. every 440 timepoints =~46 s) (Fig. 2), as a trade-off between motion tracking and the resulting navigator-introduced signal transients in the time series (see below). Following every navigator acquisition, the navigator was reconstructed and registered to the previous one in the series in real-time.
The field of view (FOV) of both the navigator and target sequence was updated based on the estimated translation and rotation parameters (1 s waiting time). The required time gap in the line-scanning acquisition for recording and real-time processing of the PMC navigator introduces a consistent transient signal due to the temporary loss of the steady-state of the transversal magnetization. This is observed as a T 1driven return to steady-state, here dubbed 'T 1 -transient'.
Three possible navigators were compared (Table 2): a highlyaccelerated surface-coil-receive fat-navigator only covering the back of the head (surf fat-nav), a slower whole-head, transmit-coil-receive fatnavigator (vol fat-nav) and a surface-coil-receive water-excitation navigator (surf wat-nav), used to reduce the amplitude of the T 1 -transient signal. The order of acquisitions was randomized across participants.
We acquired one run of functional data (6 min 20 s) using the 5-echo acquisition and the same visual task described above. For each scan, we applied the NORDIC-denoising step in the reconstruction of the data.
We investigated three ways of managing the gaps and T 1 -transient: 1. Regressing out the T 1 -transient during the GLM analysis (regressed): 30 points around every T 1 -transient of each timecourse were averaged and fitted to an exponential decay (a − be − cx Each voxel timecourse was then submitted to the T 1 -transient regression in the GLM analysis; 2. Interpolating the points corresponding to the T 1 -transient by substituting 30 point in correspondence of the T 1 -transient with the average of the 10 points after it (imputed); 3. Analyzing the T 2 * estimates which are implicitly non-sensitive to baseline T 1 effects (T 2 *PMC fit); The time points during which the navigators were acquired and the following pause were removed from the GLM's regressors that model the visual task.
We evaluated the t-stats and tSNR values along the line, to find the optimal acquisition and analysis strategy. Specifically, for our group comparison, we plotted the tSNR evaluated from the whole timecourse following the dummy acquisition, and the mean t-stats along the line for all the participants, for all the different ways of managing the T 1 -transient and for every PMC navigator method. We compared those results using an ANOVA test, to check if the imputed data showed a significant improvement compared to the other methods and the echo combined data (raw).
The navigator method allows for motion measurements. Pilot experiments on four participants with different navigator acquisitions and length of the scans showed average frame wise displacement of ~0.2 mm. This results in motion patterns that are well within the range that can be corrected for using fat navigator-based motion correction (Andersen et al., 2019;Gallichan and Marques, 2017).

Comparison of denoised multi-echo with single-echo line-scanning
In the very last section, for one participant only, we acquired a functional scan with the same functional task as in Section 2 to demonstrate the improvements we made to our previous implementation of line-scanning (Raimondo et al., 2021). Here, we compared the single-echo gradient echo line-scanning sequence from Raimondo et al. and the 5-echo multi-echo version with SoS echo combination and NORDIC-denoising in terms of t-stats.  (3a), the placement of the saturation slabs (3b), the line signal distribution image (3c) and finally, an example of multi-echo line-scanning acquisition (for a 5-echo acquisition), for every echo separately (3d) and for the combined version (through SoS) (3e). Note the decreasing signal intensity with increasing TE in consecutive echoes.

Multi-echo acquisition and echo combination approach
We evaluated t-stats for every acquisition and echo combination method, and we averaged the maximum t-stat across participants (Fig. 4a) and the mean value of t-stats in the 11-voxel gray matter (GM) ROI centered on that maximum (Fig. 4b). As the echo combination methods used the same data, the variance was higher between acquisition types than between echo combination methods. There were no statistically significant differences for different numbers of acquired echoes, though visual inspection showed higher t-stats, both mean and peak, for the 5-echo acquisition. We found significantly higher maximum and mean t-stats for SoS and wtSNR echo combination compared to the T 2 * fit method (Student t-test, p < 0.05). Regarding the tSNR (evaluated across the whole timecourse), we averaged the values of two different ROIs: GM ROI and white matter (WM) ROI ( Fig. 4c and d, respectively). In both ROIs, the SoS echo combination gave slightly higher tSNR compared to the other two methods. Considering that the second acquisition (5-echo) led to the highest mean and maximum tstats across groups, we used this approach for the rest of the comparisons. For the reconstruction, we selected the SoS because of the slightly higher tSNR and ease of implementation.

NORDIC denoising for multi-echo line-scanning
In Fig. 5, the comparison of line-scanning data without (a) ('Standard') and with (b) ('NORDIC') denoising are shown for one exemplar participant. Note that the signal found outside the brain (i.e. mostly resulting from thermal noise) was notably reduced after NORDIC denoising. The BOLD response to the visual task was preserved after NORDIC, as demonstrated for a single voxel time course in Fig. 5c. NORDIC notably improves tSNR and the distribution of t-stats ( Fig. 5d and e, respectively). A scatter plot (Fig. 5f) of the t-stats shown in Fig. 5e, showed that the voxels with positive t-stats became more positive; they are found above the unity line on the positive x-axis (black line, Fig. 5f). There are very few voxels for which the absolute t-stat became smaller. The effect of denoising was highly consistent: in every individual participant, we saw a substantial improvement in both t-stats and tSNR by at least 200%. Table 3 showed the group average of tSNR mean, tstats mean and t-stats maximum, within an 11 voxel ROI surrounding the voxel showing the highest t-stats.

Table 2
Scan parameters of the interleaved navigators scans for the prospective motion correction during line-scanning acquisition, using the interleaved scanning architecture (MISS).  Fig. 6a, b and c showed for a single voxel the raw data timecourse for the three acquisitions; surf wat-nav, vol fat-nav, and surf wat-nav. The T 1transient signal due to the time needed for navigator coregistration and acquisition update is highlighted in Fig. 6a. The selected single voxel was at the same nominal location in the brain in an area with large taskdriven responses. Note that the faster, more undersampled, navigators (surf wat-nav and surf fat-nav) resulted in reduced T 1 -transient signal amplitude due to the shorter acquisition gap. The utilization of water excitation for the navigators reduced the T 1 -transient signal amplitude even further. The regressed, imputed and T 2 *PMC fit corrected data are shown in Fig. 6d, e and f. For this particular voxel, BOLD responses were visible in all timecourses, despite the T 1 -transient. For all acquisitions, the T 1 -transient signal was much reduced after this GLM-based signal regression and completely disappeared in T 2 *PMC fit corrected data. T 2 *PMC fit data is plotted separately as the resulting T 2 * time course has physical units in ms. Fig. 7 showed box plots for the mean tSNR across participants, evaluated on the whole timecourse, across the whole line, for the PMCinduced T 1 -transient correction approaches. With each correction approach displayed separately: (a) the water navigators acquired with surface coils (surf wat-nav), (b) the fat navigators acquired with the transmit coils (vol fat-nav) and (c) the fat navigators acquired with surface coils (surf fat-nav). Fig. 8 showed the box plots for the mean t-stats across participants across the whole line, for the PMC-induced T 1 -transient correction approaches. As seen in Fig. 6, the navigators acquired with surface coils (surf wat-nav and surf fat-nav) offered slightly higher tSNR (Fig. 7) and tstats ( Fig. 8) compared to the navigators acquired with the transmit coil (vol fat-nav). Analysis of variance (ANOVA) test was used to assess whether one of the ways of dealing with the T 1 -transient showed a significant improvement in either t-stats or tSNR compared to the other two and raw data. We found significant effects on t-stats for the surf wat-nav (F 3,15 = 7.4, p = 0.03), vol fat-nav (F 3,15 =12.8, p = 0.01), and surf fat-nav (F 3,15 = 9.3, p = 0.01). Using a pairwise post-hoc T-test, we found no significant contrasts.

Prospective motion correction using water and fat excitation navigators
Regarding the tSNR we found F 3,15 = 26.7, p = 0.003 for surf watnav, F 3,15 = 31.9, p = 0.002 for vol fat-nav and F 3,15 = 25.2, p = 0.004 for surf fat-nav. The post-hoc T-test showed that the imputed and regressed methods resulted in higher tSNR values compared to the T2* PMC fit method, as well as raw data, for every navigator acquisition strategy.

Comparison of denoised multi-echo line-scanning with previous single-echo line-scanning
As a final comparison, in Fig. 9 we showed the line-scanning data acquired with single-echo (a), NORDIC-denoised multi-echo (b) as well as a single voxel timecourse for both acquisitions (c).
We observed an increase in signal intensity when multi-echo data are acquired, compared to the single-echo acquisition, and a significant decrease of noise in the multi-echo denoised single voxel timecourse, as demonstrated from the averaged tSNR along the line, of 4.0 for the single-echo line-scanning and 36.2 for NORDIC-denoised multi-echo. Moreover, from Fig. 9d where t-stats distributions are shown, we observed an improvement in t-stats when multi-echo NORDIC-denoised data are used. On average, across the line, a mean t-stats value of 3.8 ± 8.4 (mean±standard deviation) was found for single-echo line-scanning, and 25.1 ± 20.5 for NORDIC-denoised multi-echo line-scanning.

Discussion
Line-scanning fMRI is a novel technique for high spatiotemporal resolution fMRI in humans with multiple potential applications such as cognitive neuroscience, including layer and columnar imaging, but also clinical studies on, for instance, small vessel disease. Here, we report on three improvements to our first implementation of line-scanning (Raimondo et al., 2021) to increase the sensitivity and flexibility of line-scanning and mitigate the effects of motion: 1) multi-echo acquisitions, 2) NORDIC denoising and 3) real-time motion correction using interleaved navigators.
Multi-echo fMRI is known to increase SNR and BOLD contrast-tonoise ratio (CNR) and decrease sensitivity to physiological noise (Kundu et al., 2012;Poser et al., 2006). We found that a scheme with 5 echoes showed the highest sensitivity. This is likely the result of the known interplay between SNR (and tSNR) and TE-dependent BOLD CNR as reflected by the higher t-stats from the visual task. Other multi-echo studies have also opted for 5 echoes (Hesse et al., 2009;Poser et al., 2006), suggesting that this number of echoes is a suitable balance between exploiting the power of the multiple echo acquisition (in terms of Fig. 4. (a) mean value of t-stats within the GM ROI for every multi-echo (N echoes) line-scanning acquisition, averaged across participants; Shaded areas correspond to the standard error over participants. (b) maximum value of t-stats within the ROI, averaged across participants; (c) mean tSNR within a GM ROI for every acquisition, averaged across participants; (d) mean tSNR within a WM ROI for every acquisition, averaged across participants. We found significantly higher max and mean t-stats for SoS and wtSNR echo combination compared to the T 2 * fit method. The SoS echo combination gives slightly higher tSNR compared to the other two methods. The 5 echoes acquisition led to the highest mean and maximum t-stats.
T 2 * range broadening) while maximizing the readout length without too much time spent in the risetime of the readout gradients. Further increasing the number of echoes in the same readout-time pushed the first echo to earlier echo times, resulting in the highest tSNR for 9 echoes, albeit with reduced CNR (t-stats) compared to the other multi-echo schemes, possibly due to the increased gradient ramp time.
The tSNR weighted echo combination approach (wtSNR) showed a similar behavior of tSNR profile with increasing echoes compared to SoS (Fig. 4c). The T 2 * fit echo combination approach consistently showed reduced t-stats and tSNR, reflecting the challenges of obtaining a good T 2 * fit while retaining the necessary bandwidth for our current spatiotemporal resolution. Note that the tSNR increase with respect to the increasing number of echoes appears larger in WM rather than the GM ROI (Fig. 4d). WM is known to contain less physiological noise than GM (Krüger and Glover, 2001), so the larger tSNR increase likely reflects a reduction in the contribution of thermal noise.
Regarding the echo combination methods, we investigated three echo combination strategies. Note that many other combination approaches exist, some of which might exploit the benefit of multi-echo acquisitions better (Heunis et al., 2021;Kundu et al., 2017), though most do not solve the T 1 -transient that we have to deal with in the motion-corrected data.
To further optimize line-scanning, we aimed to reduce thermal noise by adapting NORDIC denoising for multi-echo line-scanning (Vizioli et al., 2021). High spatial resolution line-scanning data is likely dominated by thermal noise, as opposed to physiological noise. We observed that most of the principal components (99.7%, on average across participants) after SVD were removed from the data, suggesting that thermal noise is indeed dominant and needs to be removed. T-stats increased after denoising, with respect to the standard reconstruction (no denoising), for every participant. The "hard-thresholding" we used in the adapted NORDIC denoising is very liberal (due to the large number of components removed from the data), hence, as for any denoising techniques, one should be wary of potential biases that can be introduced (Kay, 2022;Vizioli et al., 2021). This is particularly important when less strong stimuli (leading to smaller responses) are used. In these cases, other thresholding approaches can be adopted (such as approaches involving noise scans or g-factor maps when parallel imaging is used, or approaches involving cross-validation (Kay et al., 2013) to evaluate a threshold that allows removing specific noise sources while retaining signal components.
Importantly, besides the t-stats, the tSNR also strongly improved when NORDIC-denoising was applied, in agreement with the findings of (Vizioli et al., 2021). Note that line-scanning presents some analogies with electroencephalography (EEG) data; both have a one-dimensional nature, high noise levels, and temporal resolution in the order of ms. This suggests that denoising methods used in EEG could also be investigated for line-scanning. Specifically, some non-linear approaches denoised line-scanning fMRI data after noise removal, (c) single voxel timecourse for standard linescanning data (red line) and NORDIC-denoised data (blue line), together with the GLM model following the visual task, (d) tSNR comparison of the standard line-scanning data (red line) and the NORDIC-denoised data (blue line), (e) t-stats distributions for the standard data (red) and for the NORDIC denoised data (blue), and (f) scatter plot of standard t-stats and NORDIC-denoised t-stats for one representative participant. The black line indicates the unity line. NORDIC denoising notably improves tSNR and increases the t-stats upon a visual task. The BOLD response to the visual task is preserved after NORDIC, while NORDIC notably improves tSNR and t-stats distribution.

Table 3
Comparing NORDIC denoising to non-denoised (standard) line-scanning in terms of tSNR and t-stats upon a visual task, group results: average tSNR mean, mean t-stats and maximum t-stats in the 11 voxels ROI; mean and standard error over participants. We observe that NORDIC denoising improves tSNR and t-stats by > 200% for line-scanning.  (Aydin, 2008(Aydin, , 2009, as well as ICA-based methods, particularly on short-time Fourier transforms of EEG signals (Hyvärinen et al., 2010), appear to efficiently denoise EEG data. Line-scanning data are by definition 1D fMRI recordings in the spatial domain. The 1D nature renders line-scanning more susceptible to motion compared to standard fMRI, which is exacerbated by the fact that volume coregistration cannot be applied as a post-hoc correction method. Moreover, if motion occurs during the functional line-scanning acquisition, any drift out of the selected FOV is impossible to detect and fix. Here, we implemented a motion correction procedure using interleaved large FOV navigators to track and correct the acquisition in realtime, albeit at the temporary loss of both fMRI samples and the signal steady-state, which induces a consistent T 1 -transient signal warranting the reported correction schemes. We employed three different navigators: two surface coil navigators (surf wat/fat-nav), which have the advantage of being faster due to the possibility of strong SENSE acceleration by parallel imaging, at the cost of brain coverage and possibly reduced coregistration quality, and one transmit coil navigator (vol wat/  fat-nav), which was much slower (2-3 times slower than surface coils navigators), but whole brain movements could be tracked, favoring the coregistration.
The T 1 -transient on the timecourse can be minimized by reducing the navigator acquisition time, i.e. by employing high-undersampling afforded by dense surface receive arrays. Surface-array-recorded navigators (surf wat/fat-nav) provided large gains in t-stats and tSNR compared to a whole-head but slower navigator (vol wat/fat-nav), acquired with the transmit coil. Using water excitation rather than a fatbased excitation navigator leads to lower T 1 -transient amplitudes, as the excitation of the navigator counteracts the T 1 -driven magnetization recovery. Note that the water-based navigator images are less sparse   Raimondo et al. (2021), (b) NORDIC-denoised multi-echo line-scanning fMRI data (c) single voxel timecourse for single-echo line-scanning data (red line) and NORDIC-denoised multi-echo data (blue line), and (d) t-stats distributions for the single-echo data (red) and for the NORDIC denoised multi-echo data (blue). Substantial improvements in signal quality (tSNR single-echo = 4.0, tSNR NORDIC denoised multi-echo = 36.2; averaged across the line) and t-stats (t-stats single-echo = 3.8, t-stats NORDIC denoised multi-echo = 25.1; averaged across the line) upon a visual task is demonstrated for the NORDIC denoised multi-echo line-scanning. than fat images, resulting in potentially increased SENSE artifacts. However, our current undersampling factor of four in an array of 32 receive coils, did not prove to be problematic. In terms of signal analysis, large gaps and a T 1 -transient signal are introduced to the timecourse. From the analysis of all participants, we found that a simple interpolation of the points corresponding to the T 1 -transient was the optimal option for managing the issue, in terms of both tSNR and t-stats. Even though the interpolation method is the simplest, it completely removes the T 1 -transient distortions, (unlike our regression method, which only reduces the T 1 -transient signal), hence enhancing both the tSNR and tstats. The T 2 *PMC fit method yielded noisy timecourses and, while the T 1 -transient signal is fully eliminated, it does not fulfill the requirement of high tSNR. A more sophisticated regression might bring the results closer to the imputed data quality while retaining more of the original timecourse. However, the limited tSNR of linescan acquisitions and the sharp peak of the T 1 -transients to be removed would complicate such regressions.
We already acknowledged that the presence of the T 1 -transient is a limiting factor in the motion-corrected data, but it is a great tool when scanning specific groups, such as patients, young adults and non-trained participants (i.e. participants that are completely naive to scanning), as well as to validate the sensitivity of the method to motion. In general, considering the necessary pause for navigator acquisition and the resulting T 1 -transient effect to perform PMC, we recommend using image-based prospective motion correction only when non-trained participants are involved. Highly trained participants are capable of staying still within 250 µm for remarkably long periods of time (Zimmermann et al., 2011), if necessary supported externally, such as through head fixation.
From our last comparison, on a single participant, we observed a substantial improvement in data quality and obtained functional responses for NORDIC denoised multi-echo line-scanning compared to the original single-echo line-scanning (Raimondo et al., 2021), in terms of both t-stats and tSNR. The NORDIC-denoising proved to be a great tool for decreasing thermal noise, while the use of multi-echo data per se offers an increase in MR signal and more freedom in the processing of the data, such as different possibilities of echo combination strategies or physiological noise regression. We decided not to add the PMC in the comparison with the single-echo line-scanning version to avoid the T 1 -transient, which would bias the comparison.
The presented improvements are relatively straightforward ways to increase the data quality and make line-scanning fMRI more generalizable and open for new neuroscientific questions, as well as possible clinical research. Specifically, with a double session approach we would aim for a subject-specific line planning, in order to investigate the hemodynamic responses function across cortical depth in patients with small vessel and sickle cell diseases, compared to healthy participants (Afzali-Hashemi et al., 2021;DeBaun and Kirkham, 2016;van den Brink et al., 2022;Zwanenburg and van Osch, 2017).

Conclusion
Line-scanning is a powerful fMRI technique to detect BOLD responses at ultra-high spatial and temporal resolutions. Here, we added multiecho readouts, NORDIC denoising, and real-time motion correction. We suggest a 5-echo multi-echo acquisition with NORDIC-denoising for line-scanning fMRI in the visual cortex. For non-trained participants, we recommend using prospective motion correction and we suggest interpolating the time points corresponding to the T 1 -transient time in order to correct for it. Using multi-echo readouts and NORDIC denoising for line-scanning, we found a substantial increase in tSNR and t-stats upon a visual task compared to the original single-echo line-scanning protocol.

Funding
This work was supported by Royal Netherlands Acadamy for Arts and

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
Data will be made available on request.