Improved brain stimulation targeting by optimising image acquisition parameters

Functional connectivity analysis from rs-fMRI data has been used for determining cortical targets in therapeutic applications of non-invasive brain stimulation using transcranial magnetic stimulation (TMS). Reliable connectivity measures are therefore essential for every rs-fMRI-based TMS targeting approach. Here, we examine the eﬀect of echo time (TE) on the reproducibility and spatial variability of resting-state connectivity measures. We acquired multiple runs of single-echo fMRI data with either short (TE = 30 ms) or long (TE = 38 ms) echo time to investigate inter-run spatial reproducibility of a clinically relevant functional connectivity map, i.e., originating from the sgACC. We ﬁnd that connectivity maps obtained from TE = 38 ms rs-fMRI data are signiﬁcantly more reliable than those obtained from TE = 30 ms data sets. Our results clearly show that optimizing sequence parameters can be beneﬁcial for ensuring high-reliability resting-state acquisition protocols to be used for TMS targeting. The diﬀerences between reliability in connectivity measures for diﬀerent TEs could inform future clinical research in optimising MR sequences.


Introduction
Resting-state functional MRI (rs-fMRI) is extensively used in neuroimaging to assess functional connectivity within brain networks ( Biswal et al., 2010 ) and low-frequency fluctuations in general ( Woletz et al., 2019 ;Yu-Feng et al., 2007 ). It is based on recording the blood oxygenation level-dependant (BOLD) signal changes related to spontaneous neural activity ( Biswal et al., 1995 ). Functionally connected areas are identified by determining spatially distinct brain regions that demonstrate synchronous BOLD fluctuations at rest. These so-called resting-state networks (RSNs) have provided insight into typical cerebral functioning, as well as changes in neurological and psychiatric disorders ( Canario et al., 2021 ). For instance, in Major Depressive Disorder (MDD) multiple RSNs show altered connectivity patterns and disrupted correlation with executive function ( Connolly et al., 2013 ;P. Liu et al., 2020 ).
Recently, functional connectivity analysis from rs-fMRI data has been used for defining the optimal cortical target in the therapeutic application of non-invasive brain stimulation using transcranial magnetic stimulation (TMS) ( Cash et al., 2019 ;Cole et al., 2020 ;Fox et al., 2012 ;Weigand et al., 2018 ) and to better understand its mechanism of action ( Tik et al., 2017 ). TMS is an FDA-approved treatment for pharmacotherapy-resistant MDD and other neurological and psychiatric disorders ( George et al., 1997 ;O'Reardon et al., 2007 ) and is based on short pulses of strong, localised magnetic fields that are used to stimulate brain areas. For MDD treatment, the typical target region is the left dorsolateral prefrontal cortex (DLPFC), an area involved in mood regulation and executive functioning. While standard approaches for defining the cortical target within the DLPFC use scalp-based measures like distance to the motor cortex hand area ( Pascual-Leone et al., 1996 ), recent studies have indicated that methods identifying TMS stimulation targets using personalised brain connectivity patterns derived from rs-fMRI yield increased response and remission rates (for a review see Cash et al., 2021 ).
Of particular importance is the connectivity between the subgenual anterior cingulate cortex (sgACC) and the DLPFC, as sgACC-DLPFC connectivity has been shown to be associated with treatment response ( Cash et al., 2019 ;Fox et al., 2012 ;Weigand et al., 2018 ) More specifically, treatment outcome is improved for DLPFC stimulation targets showing anti-correlated functional connectivity with the sgACC. Typically, connectivity between sgACC and DLPFC is assessed in a seed-based approach where the correlation between the rf-fMRI time-series in a seed region and the rest of the brain is calculated ( Biswal et al., 1995 ) The resulting maps, however, are strongly affected by noise in the seed and target region time courses, clearly showing the need for adequate noise reduction procedures ( Bright and Murphy, 2013 ;Pruim et al., 2015 ;Weissenbacher et al., 2009 ).
Reliability of RSN results is a sine qua non condition for every rs-fMRI-based TMS targeting approach. In a recent meta-analysis, Noble et al. (2019) found only modest reliability of RSNs (mean ICC = 0.29) across 25 RSN reliability studies examined. They reported higher reliability with shorter inter-session intervals and increased data size. With regards to TMS targeting, a recent study showed that the spatial variability in the target location is lower when a combination of two scans from the same day is used for calculating RSNs ( Ning et al., 2019 ). While this confirms the overall principle of longer rs-fMRI runs yielding increased reliability, the effects of other imaging sequence factors on RSN reliability are less evident.
The quality of the functional images acquired during rs-fMRI critically depends on several sequence protocol parameters. Images used for functional MRI are typically acquired with gradient-echo echo-planar or spiral sequences yielding T2 * -weighted data sets. For these images, signal intensities critically depend on the selection of the echo time (TE), i.e., the time between excitation and signal reception. As the MRI signal decays after excitation, shorter TEs will yield higher signal intensities and thus higher signal-to-noise ratios (SNRs) in the images acquired. On the other hand, BOLD contrast increases with longer TE ( Havlicek et al., 2017 ). While accurate modelling the change of BOLD signal over TE needs separate treatment of extravascular and intravascular BOLD effects ( Uluda ğ et al., 2009 ), for longer TEs it can be approximated by a mono-exponential decay function: With S 0 representing the initial signal amplitude. Assuming this monoexponential decay, the maximum signal change Δ will occur at TE = T2 * . While typical T2 * values in adult grey and white matter at 3 Tesla range from around 40 to over 65 ms ( Fera et al., 2004 ;Peters et al., 2007 ;Wansapura et al., 1999 ;Windischberger and Moser, 2000 ), echo times of 30 ms or less are typically used in adult resting state fMRI studies at 3T ( Greicius et al., 2003 ;D. Liu et al., 2013 ;Mouraux et al., 2011 ;Wei et al., 2018 ;Zou et al., 2015 ). Such short-TE acquisition yields images with higher SNR, but BOLD contrast will be reduced and RSN estimation might be suboptimal. It is important to note that in contrast to T2, T2 * not only depends on the tissue type but also on the static, local magnetic field inhomogeneities. In areas of strong field changes like the sgACC, MR signals decay much faster compared to homogeneous fields due to intravoxel dephasing effects ( Merboldt et al., 2001 ). A straightforward approach to reduce field inhomogeneities within a voxel is to reduce voxel sizes, i.e., increase spatial resolution. Several studies have shown that such increased spatial resolution acquisition schemes reduce intravoxel dephasing effects and effectively regain signal in areas affected ( Geissberger et al., 2020 ;Robinson et al., 2004 ;Sladky et al., 2013 ).
Multi-echo acquisition represents an elegant alternative approach as images are recorded with multiple different TEs within the same repetition time ( Boyacio ğlu et al., 2015 ;Dipasquale et al., 2017 ;Lynch et al., 2021 ;Preibisch et al., 2015 ). However, if high spatial resolution is required, the use of multi-echo acquisition is limited as k-space trajectory durations increase.
In the current study, we examine the effect of TE on the reproducibility and spatial variability of resting-state connectivity measures. We acquired multiple runs of single-echo fMRI data with either short (TE = 30 ms) or long (TE = 38 ms) echo time to investigate inter-run spatial reproducibility of functional connectivity maps originating from a sgACC seed. For comparison purposes, we also examined the motor functional connectivity network as an example for an area without strong magnetic field inhomogeneities.

Subjects
Fifteen healthy right-handed volunteers (9 female, 6 male, mean age: 26 years; std: 7.3) participated in this study. None of the participants reported any history of neurological or psychiatric disorders. Participants were recruited from the general public and were financially reimbursed for their participation. The study protocol was approved by the ethics board of the Medical University of Vienna and conducted in accordance with the declaration of Helsinki. All subjects gave written informed consent.

Data acquisition
Data was acquired on a Siemens Prisma 3T whole-body MR scanner using a 64-channel head coil with the CMRR multiband EPI sequence ( Moeller et al., 2010 ). We acquired resting-state data with two different TEs, i.e., TE = 30 ms and TE = 38 ms, while keeping all other acquisition parameters the same (TR = 2000 ms, flip angle 77°, MB = 2, slice thickness 2 mm, 170 vol, 84 slices, voxel size = 1.8 mm 3 ). We acquired three runs per TE in a counterbalanced order, resulting in a total of six runs for each subject. Subjects were instructed to let their minds wander, keep their eyes open and fixate on a white cross presented on a black background. The subject's head was fixated with additional MRcompatible pillows in order to prevent motion during scanning. All runs were acquired within one session. Each resting-state run lasted 6 min 16 s. An additional anatomical T1-weighted image was acquired for registration purposes (MPR; TR/TE = 2100 ms/3.67 ms, flip angle 8°, voxel size = 1 mm 3 ).

Functional connectivity analysis
For each subject, functional connectivity maps of the sgACC and the motor network were calculated separately for each run. A seed-based approach was used ( Biswal et al., 1995 ) and Pearson's correlation coefficients were calculated for each voxel and converted to z-scores with Fisher's r-to-z transform. To demonstrate robustness, spatial smoothing of the resulting functional connectivity maps with a Gaussian kernel was varied from 4 mm to 12 mm in steps of 2 mm. The upper bound of 12 mm was selected as it has previously been shown to result in more stable and reliable target estimates ( Ning et al., 2019 ).

Group analysis
Resulting single-subject connectivity maps were used for group analysis. Separate flexible factorial designs (as implemented in SPM12) were computed for the sgACC ROI and the M1 ROI, each representing one factor with two conditions (TE = 30 ms and TE = 38 ms). Linear regression was performed at each voxel, using generalised least squares. Contrasts were calculated for TE = 30 ms and TE = 38 ms as well as the contrasts between the TEs, i.e., TE 30 > TE 38 and TE 38 > TE 30. Threshold of the t-statistics was set to p < 0.05 and cluster-wise FWE correction was performed (initial cluster defining threshold p < 0.001). The resulting FWE corrected t-maps from the M1 seed for the contrast of TE = 30 ms and TE = 38 ms were added and binarized to create a motor network masque.

Intraclass correlation of connectivity maps
Intraclass correlation (ICC, Shrout and Fleiss, 1979 ) was used to quantify the reliability of the RSNs for TE = 38 ms and TE = 30 ms. ICC measures the extent of consistency, agreement, or reliability of an effect (e.g., BOLD response) across two or more measures. Whole-brain voxel-wise ICC values were calculated using the 3dICC program as included in AFNI . The resulting ICC maps were assessed at whole-brain level and within specific masks. For the sgACC, we were especially interested in the ICCs within the DLPFC masque, whereas for the motor network we used the motor network masque created based on the group analysis. The analysis was performed on the unthresholded individual functional connectivity maps in normalised MNI space. In accordance with Fleiss (1986) , we denote ICC values below 0.4 as poor, from 0.4 to 0.59 as fair, 0.60 to 0.74 as good and above 0.75 as excellent.

Signal-to-noise ratio
In addition to ICCs, voxel-wise temporal signal-to-noise ratio (tSNR) was calculated as the mean of time series divided by its standard deviation . The equation for a time series is defined by: Subject-wise tSNR maps were computed separately for each run and TE. Whole brain tSNR maps were first averaged across runs and across subjects to obtain mean tSNR maps for both TEs. In addition, the standard deviation (SD) across subjects was calculated.

Results
None of the subjects reported difficulties in undergoing the six resting-state fMRI runs. After each resting-state run, subjects were briefly asked for their status. All subjects confirmed that they stayed awake during the entire duration of the scanning session. In addition, none of the scans were excluded based on subject motion.

Temporal SNR
The tSNR calculated from the preprocessed unsmoothed images showed that average whole-brain tSNR was higher for the TE = 30 ms data set ( M = 45.76, SD = 15.04) compared to TE = 38 ms data set ( M = 39.99, SD = 12.98). When comparing the mean tSNR values for each subject using a paired t -test, the values for TE = 30 ms were found to be significantly higher compared to TE = 38 ms ( t = 19.20, p < 0.001). The average tSNR for voxels within the sgACC seed was 17.29 (SD = 10.44) for TE = 30 ms and 14.70 (SD = 8.68) for TE = 38 ms. Compared to sgACC regions, tSNR values were higher in left M1 seeds. The average tSNR within M1 seed was 59.47 (SD = 4.99) TE = 30 ms and 55.15 (SD = 5.76) for TE = 38 ms.

Resting-state networks
Group-level connectivity maps for the sgACC and M1 seed are shown in Figs. 1 and 2 , respectively. Results reported here are based on functional connectivity maps smoothed with a 12 mm Gaussian kernel. Figures depicting connectivity maps for the remaining smoothing kernels (4 mm, 6 mm, 8 mm, and 10 mm) can be found in the Supplementary Materials. The brain areas included in the resting-state connectivity network for the sgACC network were identical for the short and long TE. No statistically significant differences were found between the sgACC connectivity maps calculated from TE = 38 ms versus TE = 30 ms. The sgACC showed negative connectivity to regions including caudate, thalamus, middle cingulate gyrus, supramarginal gyrus, angular gyrus and importantly the bilateral medial and dorsolateral prefrontal cortex. Regions with positive connectivity to the sgACC included medial prefrontal cortex, precuneus and posterior cingulate cortex.
Also for the motor network, no statistical differences were found for connectivity patterns obtained from TE = 30 ms and TE = 38 ms rs-fMRI data. Positive connectivity to the left-M1 seed region was found in right primary motor cortex, as well as supplementary motor area (SMA) and bilateral superior temporal gyrus. Negative connectivity was found in subcortical brain areas including the thalamus ( Fig. 2 ).

Reliability
Reproducibility of the RSNs derived from the sgACC and the motor cortex were assessed by calculating ICC values over the runs and TEs. A comparison of whole-brain ICC values between the two TEs shows higher reliability for TE = 38 ms (mean/std = 0.31/0.20) compared to TE = 30 ms data (mean/std = 0.24/0.19). As TMS targeting approaches typically aim for a target within the greater DLPFC area, we specifically compared the ICC values within the predefined DLPFC masque for each TE. ICC values of the voxels within the DLPFC were significantly higher ( p < 0.001) for TE = 38 ms (mean/std = 0.35/0.17) compared to TE = 30 ms (mean/std = 0.23/ 0.18). Voxel-wise ICC DLPFC maps are displayed in Fig. 3 .
It has been previously suggested that a higher smoothing kernel may be a better strategy when using resting-state fMRI for the application of identifying cortical neurostimulation targets ( Ning et al., 2019 ). We thus also assessed the robustness of the TE-effect on functional connectivity reliability across different spatial smoothing kernels ranging from 4 mm to 12 mm with steps of 2 mm. The results of the ICC analyses from different smoothing kernels presented here are consistent with this suggestion ( Fig. 4 ).
To analyse if the higher tSNR within the M1 seed area affects the reliability, we calculated ICC values for the motor network similar to the sgACC analyses. Average whole-brain ICC values are similar at both echo times (TE = 30 ms: mean/std = 0.23/0.17; TE = 38 ms: 0.23/0.18). Within the motor network, however, reliability is slightly higher for TE = 38 ms (mean/std = 0.45/0.16) compared to TE = 30 ms (mean/std = 0.44/0.15). ICC motor network maps are displayed in Fig. 5 .
In Table 1 the distributions of voxel wise ICC values as represented by the percentage of voxels within the RSNs are summarised.

Discussion
This study investigated the effect of echo time on the reproducibility of resting-state connectivity measures. In a group of healthy subjects, three rs-fMRI scans were acquired per TE within one scanning session, and the reliability of connectivity measures was assessed. Considering the lower SNR in the sgACC region compared to more superficial cortical regions, we also tested the reliability of the motor network using a control seed region in the left primary motor cortex. Fig. 1. Group level seed-voxel based connectivity maps (12 mm smoothing) for the sgACC seed (MNI 6, 16, − 10). No statistical difference was found when comparing the connectivity patterns for long and short TE ( t = 4.74, p > 0.05). The sgACC was positively connected to the posterior cingulate cortex, the precuneus and the medial prefrontal cortex. Negative connectivity was found between the sgACC and the thalamus, angular gyrus and the DLPFC. Threshold for the connectivity maps was set to p < 0.05, FWE cluster-level correction for multiple comparison. Axial slices starting at Z = − 16, in steps of 4 to Z = 66.

Fig. 2.
Group level seed-voxel based connectivity maps (12 mm smoothing) for the motor cortex. For the seed in M1 (MNI 35, − 22, 55) we found positive connectivity to the bilateral primary motor cortex and temporal lobes. Negative connectivity was found to the precuneus and thalamus. Threshold was set to p < 0.05, FWE clusterlevel correction for multiple comparison. Axial slices starting at Z = − 10, in steps of 6 to Z = 62.  Our results clearly demonstrate that sgACC-DLPFC connectivity maps obtained from TE = 38 ms rs-fMRI data are significantly more reliable than those obtained from TE = 30 ms data sets. This increase in reliability was found despite lower SNR values at TE = 38 ms compared to TE = 30 ms.
A direct statistical comparison of connectivity maps at the two TEs did not show any significant difference. This is also a remarkable result, as the lower SNR at longer TEs would suggest a reduction of sensitivity to detect RSNs. Regarding inter-run reproducibility, we calculated voxel-wise ICC values and showed that for both the sgACC RSN and the motor RSN average reliability measures were higher for the longer TE compared to the shorter TE.
Following the ICC classification from Fleiss (1986) , the average reliability for both TEs was poor (ICC = 0.35 for TE = 38 ms and ICC = 0.23 for TE = 30 ms). These results are in line with Noble et al. (2019) who found modest reliability (ICC = 0.29) for several RSNs. However, for the area of interest within the sgACC RSN, i.e., the greater DLPFC area, 31.6% of the voxels had fair reliability (ICC = 0.4-0.6) for TE = 38 ms compared to only 16.3% for TE = 30 ms. Additionally, 61.5% of the voxels within the DLPFC were classified as having poor (ICC = 0-0.4) reliability for TE = 38 ms, compared to 80.1% for TE = 30 ms. Lastly, 6.8% of the voxels within the DLPFC had good reliability for TE = 38 ms compared to 3.6% for TE = 30 ms ( Table 1 ).
In addition, a possible effect of spatial smoothing on the reliability measures was examined. Recent work has recommended larger smoothing kernels to yield higher intraindividual reproducibility ( Ning et al., 2019 ). Our results corroborate these previous findings. In general, reliability improves with a wider smoothing kernel, regardless of TE choice. In addition, functional connectivity maps calculated from images acquired with TE = 38 ms consistently show higher reliability compared to TE = 30 ms. The difference in reliability between the TEs is largest for functional connectivity maps smoothed with a FWHM of 12 mm. However, even with minimal smoothing (FWHM = 4 mm), differences in reliability between TEs are apparent ( Fig. 4 ). Based on these results, we conclude that a slightly longer TE is advised when acquiring restingstate fMRI data for the purpose of identifying brain stimulation targets based on functional connectivity. Further, spatial smoothing of the data improves the reliability of the functional connectivity pattern.
In the current study we tested a TE of 30 ms and 38 ms. However, previous literature on the reliability of functional connectivity measures suggests that the optimal TE for subcortical structures is much shorter, i.e., below 20 ms . Lynch et al. (2021) adopted a multi-echo (ME) sequence in four subjects and show improved test-retest reliability of functional connectivity measures compared to resting-state fMRI data sets where images were acquired with a single echo sequence. As an explanation for these results the authors propose the biophysical signal mechanisms of the variability in T2 * decay of different brain tissues. However, it is important to note that T2 * not only depends on the tissue type but also on the local magnetic field inhomogeneities also affect the MR signal. Therefore, it is important to minimise field inhomogeneities and their effects. In areas of the ventral brain, such as the sgACC, strong field inhomogeneities exist due to the different susceptibilities of adjacent tissues, causing signal loss through intra-voxel signal dephasing. In the image acquisition approach presented in this manuscript, spatial resolution is increased to reduce field inhomogeneity per voxel to reduce dephasing of the MR signal ( Robinson et al., 2004 ). In accordance with previous studies ( Geissberger et al., 2020 ;Sladky et al., 2015Sladky et al., , 2018Sladky et al., , 2022 we show that small voxel sizes effectively recover susceptibility-related signal loss. Utilizing this approach, Fig. 4. ICC values within the DLPFC component of the sgACC resting-state network for unsmoothed data and different Gaussian smoothing kernels (FWHM of 0 mm, 4 mm, 6 mm, 8 mm, 10 mm and 12 mm). The ICC values are consistently higher for TE = 38 ms compared to TE = 30 ms.
we show that using a longer TE in a resting-state scan lasting just under 7 min, we derive functional connectivity measures with higher reliability compared to a short TE.
Regarding TMS targeting strategies, it is important to note that assessing the stability of potential TMS targets was not within the aims of the current study. In clinical practice as well as in clinical research, there are several methods for defining personalised connectivity-guided TMS targets. For instance, ( Cash et al., 2021 ) describe a method in which contiguous clusters of DLPFC voxels at which connectivity was most anticorrelated with the sgACC were identified and the centre of gravity of the largest cluster was defined as the target coordinate. Another method described by Fox et al. (2013) computes the connectivity for half-spheres positioned at the surface of the DLPFC and selects the target based the highest anti-correlation to the sgACC. In the recent Stanford Accelerated Intelligent Neuromodulation Therapy (SAINT) protocol, two separate algorithms are adopted to identify personalised connectivityguided TMS targets, for details see ( Cole et al., 2020 ). Moreover, recently the positive peaks within the DLPFC area of the sgACC network have been targeted to understand the mechanisms of action behind TMS. Oathes et al. (2021) showed that the neural responses corresponding to the negatively correlated sites were similar to the positively correlated sites. This suggests that positive sites within the DLPFC might also have therapeutic effects. In sum, there are multiple methods to determine personalised functional connectivity-guided targets within the DLPFC. All these different methods have one important concern in common: the need for reliable connectivity estimates within the greater DLPFC area. As demonstrated in this study, acquiring rs-fMRI data sets with longer TE can yield higher reliability.
A limitation of the current study is the inclusion of just two TEs. To make a recommendation for an optimal TE length, a larger number of TEs should be sampled. It is well known that the choice for a single TE to acquire the whole brain is often based on either a value that is best on average across all brain regions or an optimised value for a specific region of interest ( Peters et al., 2007 ;Stöcker et al., 2006 ). For this reason, the recommendation on optimal TE length is limited to the comparison between 30 ms and 38 ms. Nonetheless, the pairing of these two TEs is of high importance since on the one side most studies use a TE of 30 ms, while on the other side longer TEs have been suggested for adult grey and white matter at 3 Tesla ( Fera et al., 2004 ;Peters et al., 2007 ;Wansapura et al., 1999 ;Windischberger and Moser, 2000 ). Our results suggest that the selection of TE = 30 ms is at least suboptimal compared to a slightly longer TE of 38 ms when considering reliability of connectivity estimates based on seed-voxel correlations. Along the same line, we have shown improvement in one clinically relevant RSN and a stable control network, further research is needed to assess the impact of a longer echo time on other RSNs derived from seeds throughout the brain.

Conclusion
The choice of TE in resting-state acquisitions affects the reliability of the resulting functional connectivity maps. This can influence TMS targeting approaches that base their stimulation target on its functional connectivity to other cortical areas. Our results clearly show that optimizing sequence parameters can be beneficial for ensuring highreliability resting-state acquisition protocols to be used for TMS targeting. Taken together, our study provides important insights into the differences between reliability in connectivity measures for different TEs, which could inform future clinical research in optimising MR sequences.

Declaration of Competing Interest
None of the authors have any interests to declare.

Data availability
All anonymised data and analysis codes are available upon request in accordance with the requirements of the institute, the funding body, and the institutional ethics board.