The Kantonale Ethics Komitee (i.e., regional ethics committee) of Zurich approved involvement of the human volunteers and each participant signed a written informed consent before the experiment.

Thirty (mean age = 24.8 years, std = 1.8 years, 15 male) right-handed, non-smoking, medication-free volunteers participated in the study. Each had normal eyesight and no history of neurological disorders. The MRI or the physiological data from three participants were compromised due to computer disk error or physiological sensor failure. Because the physiological data were used for correcting time-series artifacts, henceforth this paper pertains to the remaining 27 participants (14 male, mean age = 25.75 years, std = 1.77 years).

All MRI data were acquired on a 3T Philips Achieva scanner (Philips Healthcare, Best, The Netherlands) equipped with an 8-channel receive-only head coil and single-channel body transmit coil. For both the resting-state and the task-based fMRI data sets 36 axial slices were collected at 2.6 mm thickness using single-shot double-echo GE-EPI with twofold SENSE acceleration (Pruessmann et al., 1999) in the phase encoding direction and TE1/TE2 = 17/44 ms, TR = 2.6 s, FOV 200 × 200 mm², in-plane resolution 2.5 × 2.5 mm², slice gap 0.6 mm. The total EPI readout duration for each echo was 24.792 ms. For each subject the resting-state data set contained 116 image volumes while the task-based fMRI data were acquired twice, each set containing 282 image volumes. Cardiac and respiratory signals were concurrently acquired via electrocardiogram and a breathing belt, respectively, to allow for removal of fluctuations caused by cardiac pulsation and breathing from the fMRI time-series as part of the image pre-processing pipeline.

During the task-based experiment subjects in the scanner played a competitive game against a human opponent outside of the scanner. Each trial consisted of a jittered choice and feedback epoch, in which the subject in the scanner could either win money or not depending on his/her choice and the choice of the opponent outside of the scanner. This quantity is the parametric modulator of interest, denoted as “reward.” Reward in this game is contingent on the player's ability to predict his opponent's behavior on any given round in order to be rewarded. Scanned participants played 160 trials of the inspection game in two sessions of 80 trials paired with another player seated in an adjacent room. Each trial began with a fixation-cross presented for 1.25 to 8 s, followed by a decision screen featuring two pictograms describing choice options. Participants had up to 2 s to make their decision and played simultaneously. Decisions were confirmed with a red square, which was displayed for a minimum of 100 ms, and up to the time it took for both answers to be recorded. After each player made his or her decision, a fixation cross was presented for 1.25 to 8 s, following which the feedback screen was displayed for 2.5 s. The average duration of a trial summed up to 9 s. fMRI contrast sensitivity was optimized by simulating optimal task ITIs. In total, the task lasted 25 min. It is noteworthy that the “reward”-based analysis elicits robust BOLD responses in dopamergetic pathways (Nucleus Accumbens, Substantia Nigra, Ventral Tegmental Area) and orbitofrontal cortex. We direct the interested reader toward Neurosynth (Yarkoni et al., 2011) for a meta-analytic overview for the term “Reward” across 671 studies as of October 24, 2016.

The fMRI dataset used for comparing echo combination procedures was modeled with regressors capturing latent variables from a computational model described previously (Hampton et al., 2008). Briefly, Expected Value of the Chosen option refers to the subject's prediction about his/her probability of winning the given his choice. First-order prediction error refers to the subject's estimate of observing a certain action from the opponent, given the opponent's history of choices. Influence update refers to the estimated effect of the player's choice on the opponent's model of the player.

Our design matrix included the onset of the decision period with expected value of the chosen options as parametric modulator and the onset of the feedback epoch with reward, first-order prediction error, and influence update as parametric modulators. Orthogonalization was turned off in order to capture the variance uniquely explained by each factor. In addition, we included, as nuisance regressors, all six movement parameters and physiological fluctuations related to heart rate, and breathing with the procedure described below. We constructed a boxcar epoch function whose duration corresponds to the reaction time of the decision to optimally account for variability in decision time. The feedback epoch was modeled using stick functions.

For pre-processing and statistical analysis SPM12 (Wellcome Trust Centre for Neuroimaging, UCL, UK), and Matlab (The MathWorks, MA, USA) were used.

First, physiological time-series were transformed following the RETROICOR procedure (Glover et al., 2000), as implemented in TAPAS, an open source software package (http://www.translationalneuromodeling.org/tapas) that uses Fourier expansions of various orders for the phases of cardiac pulsation (3rd order), respiration (4th order), and cardio-respiratory interaction (1st order) (Harvey et al., 2008). A cardio-respiratory interaction refers to the respiratory time-series multiplied by the cardiac time-series to form the interaction term. This is done in order to account for physiological noise not captured by the respective main effects of breathing and heart rate. Subsequently both the resting-state time-series and the two runs of the task-based datasets were realigned with Matlab such that all but the first volumes with short TE were rigid body aligned to the first image and the resulting realignment parameters were used to transform the long TE images identically. Finally, second-order polynomial de-trending was performed right after re-alignment but before calculation of the weights (see “Echo combinations of resting-state fMRI data” below).

Additionally, a high-pass filter with a cut-off frequency of 1/128 s was applied to the task-based data after echo combination, co-registration and normalization and immediately before the first-level statistical analysis.

To extract voxel values specific to the gray matter (GM) segment of the brain a binary mask was created by segmenting (Ashburner and Friston, 2005) and normalizing the T₁-weighted anatomical image of each subject in SPM. Subsequently the GM segment was smoothed by a 6-mm isotropic Gaussian kernel and thresholded to contain only voxels with a GM probability of at least at 0.35. Finally, each individual's GM mask was multiplied together to arrive at a group GM mask. The histograms and scatter plots of Figures 3–5 were produced after multiplying the t-value maps, contrast maps and total variance maps voxel-wise by this final GM mask.

The realigned and de-trended resting-state dual echo dataset of each subject was combined voxel-wise as a normalized weighted average of the two echoes to provide a single time-series, where time is represented by t

where S₁(t) and S₂(t) are the signal amplitude of a given voxel in the same spatial and temporal position acquired with TE₁ and TE₂, respectively, while w₁ and w₂ are the corresponding weights and type is one of “AVE,” “BS,” “tSNR,” or “tBS” representing the four different methods of echo combinations for calculating the weights (see below). Note, the weights, w₁ and w₂, are calculated separately for each time point in Echo combination method #2 below. For the other three echo combination methods they are identical for all time points in the time series.

Echo combination #1 (AVE): The first echo combination was a simple average of the two datasets so that the weights were identical for all voxels at all-time points (i.e., w₁ = w₂ = 1).

Echo combination #2 (BS): In the second method the weights were calculated voxel-wise and for each time point (t) as

where S₁(t) and S₂(t) are defined as in Equation (1) above. This is commonly considered BOLD sensitivity (BS) weighting (Posse et al., 1999; Deichmann et al., 2002).

Echo combination #3 (tSNR): The third method used the temporal signal-to-noise ratio (tSNR) as weights calculated separately for the two TEs by dividing the temporal mean of a voxel time-series with its temporal standard deviation

where tSNR₁ and tSNR₂ are the voxel-wise temporal signal-to-noise values for short (TE₁) and long (TE₂) echo data respectively.

Echo combination #4 (tBS): The final scheme merges Echo combinations #2 and #3 as in Poser et al. (2006) and is termed temporal BOLD sensitivity (tBS)

where tBS₁ and tBS₂ are the temporal-BOLD-contrast-to-noise ratio values, while tSNR₁ and tSNR₂ are as in Echo combination #3 above.

In each of the above echo combination methods the result is a single time-series of resting-state fMRI data. BOLD sensitivity is often defined as in Equation (2) above (Deichmann et al., 2002) and used to quantitatively compare fMRI protocols. Because the datasets have no single effective TE after dual-echo combination (e.g., AVE_data), a pseudo measure of the temporal BOLD sensitivity (ptBS) was calculated as

where type is one of “AVE,” “BS,” “tSNR,” or “tBS” to represent one of the four echo combination methods above, w₁ and w₂ are the pair of weights of the corresponding type of echo combination method, S₁ and S₂ are voxel signal intensities for the short (TE₁) and long (TE₂), respectively, and both the mean and the standard deviation are calculated across time. Note that in the case BS-based echo combination method the ptBS measure will depend on the square of the voxel signal and that in this case each time point will have a unique weight.

In order to create group-averaged results the AVE_ptBS map was co-registered to the corresponding T₁-weigthed anatomical image for each subject. Given that each of the four types of ptBS maps are calculated from the same data these maps are perfectly aligned. Therefore, the coregistration parameters from AVE_ptBS were used for the other three ptBS maps. Subsequently the T₁-weighted anatomical image was normalized in SPM12 to the Montreal Neurological Institute (MNI) template and the normalization parameters were written onto each ptBS map. No spatial smoothing was applied to the normalized images.

Similarly to the resting-state data, realigned and de-trended task-based time-series with two TEs were also combined as a weighted sum using Equation (1). Echo combinations #1 and #2 were performed by using the task-based data for calculating weights, w₁ and w₂, as well as extracting S₁ and S₂.

In Echo combinations #3 and #4 participant-specific resting-state data were used to calculate the weights, w₁ and w₂, but the task-based data were used for extracting S₁ and S₂ for Equation (1). This is because Echo combinations #3 and #4 rely on tSNR (Equations 3a,b) for calculating the weights. Hence to avoid including the BOLD response of the task-based time-series into the calculation of the temporal standard deviation, and in turn the weights, the resting-state data were used. More specifically, the tSNR₁ map of TE₁ derived from the resting-state data was realigned with the temporal mean of the realigned task-based data with TE₁. These movement parameters were then applied to the tSNR₂ map to bring it in line with the corresponding task-based data with TE₂.

The four combined task-based datasets were corrected for the different timing of the slice acquisition by temporal interpolation relative to the acquisition time of the slice in the center of the volume using the standard slice time correction method in SPM12 (Sladky et al., 2011).

All four slice-time corrected time-series were normalized to the (MNI) template detailed as follows. The time-mean of the echo averaging combination was co-registered to the T₁-weighted anatomical image and the estimated parameters were applied to all volumes of all four combined time-series. Normalization to MNI space was performed based on the anatomical image using the standard method in SPM12. All volumes were interpolated to the isotropic resolution of 3 mm during the normalization step and were subsequently smoothed with an isotropic Gaussian kernel with 6 mm FWHM.

To quantitatively assess the possible advantages of the four different echo combination methods, voxel-wise mean and voxel-wise standard deviation were computed across the subjects for each of the four normalized ptBS maps. Voxel-wise ratio of the resulting mean and SD maps were calculated for all combination methods. AVE_ptBS was considered the reference to which the other 3 echo combinations were compared. To extract group-level quantitative results, each of the four ptBS maps was multiplied by the group GM mask and the mean and standard deviation were calculated from voxels within the mask

The ROIs were generated with http://www.neurosynth.org (Yarkoni et al., 2011) using a default threshold for reverse reference. The first ROI corresponds to the keyword “Reward” (671 studies, 2291 activations). The second ROI we used keywords “Default Mode” (516 studies, 18723 activations). For the third ROI we isolated the orbitofrontal cluster from the first ROI to obtain reward-related orbitofrontal activations. Each of these ROIs was multiplied with the GM mask. The resulting intersection was then applied to the t-value difference maps before generation of the histograms.

Single subject general linear model (GLM) analysis was performed separately on each of four pre-processed combined time-series datasets (i.e., AVE_data, BS_data, tSNR_data, tBS_data) of each subject. This GLM included, 8 regressors and 24 covariates (physiology and movement) per session. Analysis was performed for one contrast defined as the sum of the two predictors representing reward in the two runs contrasted to baseline.

Random-effects group-level analyses were performed separately for each echo-weighting method using the corresponding single-subject contrast maps.

We calculated voxel-wise t-value differences between the simple echo averaging and the other three echo combination methods. Scatter plots of these t value difference maps were created from voxels inside the GM mask. Histograms were created both from within the entire GM mask as well as from within the three ROIs of interest.

To further scrutinize the effects of echo combination methods, additional analyses were performed on the t-values obtained from random-effects statistical analyses (Friston et al., 2007). By definition these t-values are affected by both the mean contrast values as well as the parameter variance. In order to disentangle the two driving factors for the different echo combinations, histograms of voxel-wise differences of contrast and variance maps were also generated for the values of

where AVE represents the echo averaging method while type represents one of the other three methods for combining the dual-echo data (BS, tSNR, tBS). Both the contrast maps and the variance maps were multiplied by the GM mask before histogram calculations.

To estimate and compare the relative contributions of inter-subject vs. intra-subject variance we followed the procedure outlined in Kirilina et al. (2016), which is based on the fact that he total variance of the random effects statistical analysis is expressible as a sum of inter-subject and intra-subject variance components.

Thus far, comparisons of group-level analysis results of statistical analyses were presented. We also performed three separate statistical analyses to compare the advanced echo combination methods with AVE. First, first-level analysis was performed for each subject by treating each the four echo-weighting methods as four different sessions. To remove any overall scaling and/or change in the variance of the resulting parameter maps between the combinations, the parameters were voxel-wise divided by their inter-subject standard deviation, separately for each echo combination. Next, from the rescaled first-level parameters three contrast maps were calculated for each subject by taking the difference of the parameters representing reward in one of the advanced echo combination method vs. that of AVE. Finally, standard group-level random-effects analysis was performed on these contrast maps separately.

According to the top two rows of Figure 1, the advanced echo combination methods have clear benefit on BOLD sensitivity of the resting-state data as assessed by the type_ptBS measure (Equation 5) for the advanced echo combination methods BS and tBS but not for tSNR. Out of the four examined echo-weighting strategies the BOLD sensitivity (BS) weighting of Echo combination #2 produces the time-series with the highest sensitivity when averaged across the entire group. The tBS method of echo combination also produces the expected improvement, which confirms the findings of Poser et al. (2006). For the mean values of the four ptBS maps from within the GM mask please see Table 1.

Note however that the center two rows of Figure 1 indicate a similar tendency (especially in the GM of the brain) in the standard deviation of the ptBS maps across the subjects, forecasting that a group-level statistical analysis may not benefit from the increased BOLD sensitivity of the advanced echo combination methods (i.e., BS, tSNR, tBS). The actual group standard deviation values within the GM mask for the four ptBS maps are listed in Table 1. Calculating the voxel-wise ratio supports this hypothesis (bottom two rows of Figure 1).

Figure 2A shows the t-value map of the group-level random effects-analysis on the time-series data in which the echoes were simply averaged (i.e., AVE_data). The group-level results were highly similar when using the time-series data resulting from the other three echo-weighting strategies (data not shown). In Figure 2B the color code in each voxel indicates the echo combination method that provides the largest t-value. Within the GM mask the actual percentages were 14.5% for AVE, 24.0% for BS, 23.3% for tSNR and 38.2% tBS. Thus, neither of the four echo combination methods comes out as the major winner, not even in the orbitofrontal cortex, where multi-echo data has been shown to provide increased sensitivity (Kirilina et al., 2016).

The scatter plots in Figure 3 show the t-value of each voxel within the GM mask from the group-level statistical analyses. In each of the three subplots the t-value from the average echo combination (i.e., using AVE_data) is plotted along the horizontal axis against that resulting from the statistical analysis of each of the other three echo combination methods (BS_data, tSNR_data and tBS_data) plotted along the vertical axis. Because all points fall in the vicinity of the identity line, neither of the echo combination methods seems superior to simple averaging. It must be noted that the spread around the identity line is broader for the statistical analysis based on BS_data.

In Figure 4 the histograms of the voxel-wise t-value differences of the random-effects group-level analyses within the GM segment are centered on zero (mean ± [STD], −0.11 ± [0.50] for BS, −0.03 ± [0.22] for tSNR, and 0.03 ± [0.15] for tBS), which is their expected value if the echo combination methods provide equivalent statistical results. However, the spreads of the histograms are not the same, showing that the spatial variance of t-values is dependent on the echo combination method, with tBS weighting showing the smallest variation (i.e., using tBS_data produces statistical results, which are most like that based on AVE_data). It is worth mentioning that the distributions of histogram differences are not symmetric. Skewness also depends on the echo combination strategy with BS weighting resulting in the highest asymmetry. The tails of these histograms represent the voxels where the group-level statistical t-values differ between the echo combination strategies. The anatomical pattern of the histograms tails is presented in the bottom of Figure 4 as maximum intensity plots. Differences between the echo combination methods are only apparent in the inferior part of the brain but the spatial arrangement of these differences precludes any possible recommendation in the optimal echo combination method.

The group-level t-values are defined by both the mean contrast (i.e., the mean effect size across subjects) and the parameter variance. Further scrutiny into the t-value histograms from above indicates that, unlike the t-values themselves, the contrast is actually dependent on the method used for echo combination (mean ± [STD], −0.03 ± [0.15] for BS, −0.06 ± [0.12] for tSNR, and 0.06 ± [0.07] for tBS). Figure 5A shows that tBS weighting results in increased contrast while tSNR weighting and BS weighting have decreased mean contrast compared to echo averaging.

The reason for comparable t-values across the different echo combination methods can be concluded from the total variance differences in histograms (Figure 5B) (mean ± [STD], −0.03 ± [0.22] for BS, −0.12 ± [0.20] for tSNR, and 0.11 ± [0.11] for tBS). Based on these it can be argued that the increased spatial variance of t-values of BS weighting strategy originates from its parameter variance, as its histogram shows an asymmetric “tail” in positive direction, according to the similar “tail” on t-value differences histogram in negative direction in Figure 4.

The ratio of inter-subject vs. intra-subject variance components in Figure 6 is spatially heterogeneous. In many large connected areas across the brain the intra-subject variance component is outweighed by inter-subject variance. In some smaller connected areas intra-subject variance is the dominant variance component.

The actual value of both the intra-subject and the inter-subject variance components depends on the echo combination methods. In Figure 7 the variance components of data derived from the simple echo averaging (i.e., AVE_data) are related to those derived from the other three echo combination methods (i.e., BS_data, tSNR_data, tBS_data) where the intra-subject (Figure 7A) and inter-subject (Figure 7B) variance components vary similarly for a given echo combination method. Furthermore, these results from the group-level random-effect analyses of the task-based data corroborate that found in the resting-state data (Figure 1) in that an increase in inter-subject variance always follows the seemingly beneficial advanced echo combination method.

Figure 8 displays the histograms of the voxel-wise t-value differences of the random-effects group-level analyses within the three ROIs (reward activation areas, default mode network, and orbitofrontal cortex). All of them are centered near zero with their mean inside the [–STD, +STD] interval, which is the result expected if the combination methods are equivalent at group level (Table 2). Similarly to the results when considering the entire GM segment (Figures 3, 4), the variance of t-values is dependent on the echo combination method, with tBS weighting showing the smallest variation, and BS weighting showing the largest variation in all ROIs.

The group-level paired t-tests between the advanced echo combinations and the echo averaging also show that there is nearly no difference between combination methods: after the FWE correction (p = 0.05) there are less than 30 significant voxels scattered around all the brain, in all combination pairs (advanced vs. AVE), for both positive and negative t-contrasts.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.