Modelling Temporal Stability of EPI Time Series Using Magnitude Images Acquired with Multi-Channel Receiver Coils

In 2001, Krueger and Glover introduced a model describing the temporal SNR (tSNR) of an EPI time series as a function of image SNR (SNR0). This model has been used to study physiological noise in fMRI, to optimize fMRI acquisition parameters, and to estimate maximum attainable tSNR for a given set of MR image acquisition and processing parameters. In its current form, this noise model requires the accurate estimation of image SNR. For multi-channel receiver coils, this is not straightforward because it requires export and reconstruction of large amounts of k-space raw data and detailed, custom-made image reconstruction methods. Here we present a simple extension to the model that allows characterization of the temporal noise properties of EPI time series acquired with multi-channel receiver coils, and reconstructed with standard root-sum-of-squares combination, without the need for raw data or custom-made image reconstruction. The proposed extended model includes an additional parameter κ which reflects the impact of noise correlations between receiver channels on the data and scales an apparent image SNR (SNR′0) measured directly from root-sum-of-squares reconstructed magnitude images so that κ = SNR′0/SNR0 (under the condition of SNR0>50 and number of channels ≤32). Using Monte Carlo simulations we show that the extended model parameters can be estimated with high accuracy. The estimation of the parameter κ was validated using an independent measure of the actual SNR0 for non-accelerated phantom data acquired at 3T with a 32-channel receiver coil. We also demonstrate that compared to the original model the extended model results in an improved fit to human task-free non-accelerated fMRI data acquired at 7T with a 24-channel receiver coil. In particular, the extended model improves the prediction of low to medium tSNR values and so can play an important role in the optimization of high-resolution fMRI experiments at lower SNR levels.


Introduction
The optimization of EPI acquisition parameters is important for maximizing the sensitivity to relatively small BOLD signal changes in fMRI studies. Factors such as spatial resolution, echo time (TE), flip angle, and parallel imaging strategies can be optimized for a given magnetic field strength and receiver coil arrangement. In particular, these factors can influence the contribution of physiological noise, for example from cardiac and respiratory functions as well as motion, to the temporal noise, thereby impacting on BOLD sensitivity. The noise model introduced by Krueger and Glover [1] describes the temporal signal to noise ratio (tSNR) of an EPI time series as a function of image SNR and has been used to demonstrate that the ratio of physiological noise to thermal noise increases with image SNR [1][2][3]. The model allows the time series noise to be separated into thermal and physiological components and can be used to estimate the maximum attainable tSNR for a given set of MR acquisition parameters [2] and processing strategies [4].
In its current form, the application of this noise model requires an accurate estimate of the actual image SNR (SNR 0 ). The estimation procedure is particularly elaborate for images acquired with multi-channel receiver coils, involving detailed steps to correctly scale the image SNR [3,5]. Further steps are necessary for images reconstructed using partial parallel imaging methods [3,5]. Some of the main factors that must be taken into account are noise correlation between receiver channels, change of statistical distributions in the case of magnitude images, geometry factors in under-sampled parallel imaging and temporal autocorrelation in the k-space raw data [5]. Although a principled way to reconstruct and output images in SNR units is proposed in [5], this is not part of the standard image reconstruction available on clinical MRI scanners. Thus, SNR 0 can only be correctly estimated by custom-made image reconstruction based on k-space raw data. This approach is inefficient and in many cases even impossible, since it requires export, storage and computation of large amounts of raw data (e.g. ,70 GB for ,10 minutes of highresolution fMRI acquired using a 32-channel receiver coil).
A common approach to combining information from multiple receiver channels is to calculate the root-sum-of-squares (RSS) across channels, after reconstructing the image separately for each channel. This method is regularly used for non-accelerated image acquisition [6] and parallel imaging employing GRAPPA [7]. It has been demonstrated [8] that the image SNR can be estimated from the background noise and image signal in RSS combined non-accelerated images using a simple correction for the change in statistical distributions due to rectification of multiple channel data, i.e. the noncentral chi distribution. However, this approach is only valid if no noise correlations between channels are present, which is generally not the case for RSS composite images acquired with state-of-the-art multi-channel RF head coils [3,9]. In fact, as pointed out in [5], correctly estimating SNR has caused significant confusion and discussion in the community. For example, the SNR correction methods proposed in [10][11][12] account only for the RSS combination and magnitude operation. This correction does not take into account noise correlations in the data from the different receiver channels as explained in [8]. Typical levels of noise correlation (i.e. ,0.2 to 0.3 [9]) cause a spatial variation of the noise in the composite images [6,13,14] and lead to errors in SNR measurement when noise correlations are not taken into account. To overcome this problem, the methods described in [5] and [3] use the noise covariance matrix to account for correlations between different channels when estimating the SNR and are in line with the formulation of [6].
In this study, we extend the model proposed in [1] to allow the tSNR of an EPI time course acquired without acceleration using multi-channel receiver coils to be characterized, without the necessity for complex custom-made image reconstruction of absolute image SNR estimates. Instead, the straightforward approach described in [8] is used for estimating an apparent image SNR9 0 and the impact of noise covariance between receiver coils is accounted for by introducing a constant scaling factor k = SNR9 0 /SNR 0 , which modulates the actual image SNR 0 . To our knowledge, no other published methods can estimate the parameters of the model proposed in [1] from RSS combined magnitude data acquired using multi-channel receiver coils. We derive our proposed extended model from basic physical principles and use Monte Carlo simulations to optimize and demonstrate the accuracy and precision of the resulting fit. We validate the extended model fit using phantom data acquired at 3T and an independent estimate of the actual SNR 0 . Finally we demonstrate the improved fit of the extended model in human task-free fMRI data acquired at 7T.

Theory
Following the theory outlined in [1] the total variance or noise in an MR image time series s 2 can be described by a combination of independent raw noise s 0 2 and physiological noise s p 2 : The noise variance s 0 2 is assumed to comprise thermal noise from the subject and scanner electronics. The physiological noise s p 2 is assumed to arise from cardiac and respiratory functions, which lead to oscillatory signal fluctuations in the vascular system as well as small pulsatile movements of the brain and modulation of the magnetic field. It has been shown that s 0 increases with field strength [15] but is independent of the MR signal strength, whereas physiological noise has been shown to be signaldependent [1]. Note that in [1] the physiological noise is treated as two separate components. One describes T2* related fluctuations, which are TE-dependent and also lead to the BOLD effect and the other arises from image-to-image fluctuations such as pulsatile effects and scanner imperfections with no TE-dependencies. For the purpose of the theory developed here, we do not make a distinction between these components of physiological noise and only consider the case of physiological noise represented by s 2 p . Therefore, following on from [1] and the formulation used in [2], the temporal SNR (tSNR) can be defined as where I is the mean image signal intensity. If the actual image SNR is defined as SNR 0 = I/s 0 and the physiological noise is defined as a signal dependent function s p = lI, then the relationship between tSNR and SNR 0 is given by [2]: The parameter 1/l provides a measure of the maximum attainable tSNR for a given set of acquisition parameters. It can be estimated by measuring tSNR for different values of the actual image SNR 0 (e.g. by varying flip angle, voxel size or echo time). From equation 3 it can be seen that, in the absence of physiological noise, i.e. when l = 0, tSNR = SNR 0 . Equation 3 requires an accurate, unbiased estimate of SNR 0 as discussed above.
In the following we extend the model so that instead of the actual SNR 0 , the apparent SNR9 0 measured from RSS combined magnitude images using the approach proposed by [8] can be used. For this method the SNR9 0 is the ratio of the signal intensity in the measured object to the noise level in the background and equals the actual SNR 0 in the absence of noise correlations. For a receiver system with n coils a measure of background noise s9 0 in RSS combined images can be estimated using (equation 7 in [8]): It is important to note that although this noise estimate includes the corrections for the RSS combination and magnitude operation (i.e. the noncentral chi distribution), it does not take into account noise correlations and is therefore only accurate in their absence. Equation 4 is therefore correct for data acquired using a single-coil (which of course can not suffer from noise correlations) but not for multi-channel receiver coil data with noise correlations between receiver channels. We therefore use equation 4 to define the apparent image SNR as, SNR9 0 < I n /s9 0 , where I n is the mean signal in the RSS combined images. For low SNR9 0 further corrections must be applied to account for changes in the noise statistics. We will neglect this correction in the following since they are less than 2% for image SNR exceeding 50 and number of receiver channels not exceeding 32 [3], (i.e. the currently maximum number of channels in commercially available coils). These additional corrections for the noise statistics are only required at even lower SNR values for data acquired with less than 32 channels (Figure 1b, [3]).
In general SNR9 0 only equals the actual image SNR 0 for data acquired with a single-channel receiver coil. When n $2 receiver coils are used to acquire array coil images and if noise correlations exist between the different channels of the receiver system, the actual variance s In 2 in a voxel with composite intensity I n will deviate from the background noise variance estimate s9 0 2 [13,14]. We therefore define a constant scaling factor k to account for correlated noise, which relates SNR 0 and SNR9 0 in a voxel by: Replacing SNR 0 by SNR9 0 in equation 3 yields our proposed extended model: where k and l can be estimated from equation 6 by measuring tSNR for different values of SNR9 0 (i.e. in a similar way as for equation 3). Note that because correlated noise is spatially varying, k will also vary spatially over an image. Furthermore, we assume that k will be dependent to some extent on the loading of the coil and so will vary for phantoms and human subjects. In this study we will therefore demonstrate the validity of equations 5 and 6 in phantom and human data.

Accuracy and Precision of Extended Model Fit Estimated from Monte Carlo Simulations
The estimation of l and k from the measured SNR9 0 and tSNR data poses a non-linear problem, which was solved by a multidimensional unconstrained nonlinear minimization (Nelder-Mead) method in Matlab (2009a, The MathWorks, Natick, MA).
To address potential problems with model estimability, Monte Carlo simulations were performed to study the accuracy and precision of the estimation of k and l, and to determine how a set of five SNR9 0 values should be best distributed to obtain fit parameters with highest accuracy.
For the simulation, synthetic tSNR values were generated using the extended model (equation 6) with a realistic set of parameters: l = 1/90, k = 1.4 or 1.8 and five evenly distributed SNR9 0 values in the range of 50 to 600 or 50 to 300. These SNR9 0 values are typical for 7T and 3T data respectively; SNR9 0 values ,50 were excluded, since they would require further correction for changes in the noise distribution [8], if 32 or more channels were used. Gaussian distributed noise with zero mean and a standard deviation of 5 (typical for group studies [2,4]) was added to the tSNR values to simulate measurement noise. The model was then fitted to this set of synthetic data, yielding model parameters l and k. This process was repeated 500 times to get a reliable estimate of the mean and standard deviation (SD) of the fitted model parameters l and k.
To determine how a set of five SNR9 0 values should be best distributed to obtain unbiased and low noise estimates of l and k, the above simulation was repeated 5000 times, each with a different set of five arbitrary SNR9 0 values. Each repetition resulted in an estimate of l and k calculated using five unique SNR9 0 values, which were selected randomly from the given range based on a uniform probability distribution. From the 5000 sets of five SNR9 0 values, 250 were selected (i.e. five percent), for which the estimate of l and k had the lowest percent bias. For these 250 estimates of l and k values, the accuracy was calculated from their mean absolute bias and the precision from the SD. For visual exploration, histograms of the best 250 sets were generated to show the distribution of each of the five SNR9 0 values in increasing order ( Figure 1). The precision of the estimates based on the best sets was estimated from the mean SD over the sets. For comparison, the sets with the highest precision out of the 5000 sets and their SD were also determined, providing an approximate upper limit of the highest precision achievable (regardless of accuracy).
An additional simulation explored the model behaviour for phantom data. Temporal SNR measurements on a gel phantom are usually not significantly affected by signal dependent noise if standard acquisition parameters are used due to the high temporal stability of modern MRI scanners. Thus, the tSNR is small compared to 1/l for the phantom and the tSNR is not saturated.
To determine whether k can be estimated robustly under these conditions, the same Monte Carlo simulations described previously were repeated with the following parameters: l = 1/1800; SNR 0 = (SNR9 0 /k) = 60 or 120 or 180 for k = 1.0 to 2.0 (increased in steps of 0.1). The parameters were typical for phantom measurements as determined from independently estimated SNR 0 maps (see below) and quality assurance experiments [16,17].
The third type of Monte Carlo simulation addressed whether l and k can be estimated robustly without relevant correlation between the parameters. If the non-linear minimization problem is ill posed, the estimated parameters may be spuriously correlated, skewing the results. The correlation between parameters were estimated for a measurement approximating the in-vivo dataset acquired at 7T with SNR 0 = (SNR9 0 /k) = 80, 180, 270, 350, 450. In one simulation k was varied from 1 to 2 (in steps of 0.0001) and the effect on the estimation of l was studied. In the second simulation 1/l was varied from 80 to 140 (in steps of 0.005) and the effect on the estimation of k was studied. As in the previous simulations Gaussian distributed noise with zero mean and a standard deviation of 5 was added to the tSNR values for each simulated 1/l or k set. Robust regression analysis (iteratively reweighted least squares with the bisquare weighting function as implemented in Matlab) was used to assess whether there was a relevant correlation between the parameters.

Validation of Extended Model Fit using Phantom Data and Independent Estimate of Actual Image SNR
We used EPI time series acquired from a phantom to confirm that the scaling factor k estimated from the extended model (equation 6) equals the ratio SNR9 0 /SNR 0 (i.e. to validate equation 6) where SNR9 0 is estimated from RSS combined images and SNR 0 is an independent estimate of actual image SNR, according to [3,5]. An agar gel phantom (construction based on the Stanford Agar Phantom Recipe [16]) was scanned using a 3T whole-body MRI scanner (Magnetom TIM Trio, Siemens Healthcare, Erlangen, Germany) operated with a 32-channel RF head receive and RF body transmit coil. Three EPI runs were acquired with different RF excitation flip angles to manipulate the image SNR. Each run comprised of 205 volumes and the flip angles were 17u, 37u and 90u. A thermal noise measurement of 20 The results for different parameters l = 1/90, k = 1.4/1.8, SNR9 0 ranges = 50-600/50-300 are presented as histograms of the top 5% best sets of SNR9 0 in (A-D). Acq # orders the 5 acquisitions according to their SNR9 0 . It can be seen that these best sets sample low and high SNR' values more densely. For the small SNR9 0 range (C, D), high SNR9 0 values are more strongly represented, in order to better estimate 1/l. doi:10.1371/journal.pone.0052075.g001 EPI volumes with no RF excitation (i.e. 0u flip angle) was also acquired. These flip angles resulted in images with equally spaced image SNR levels from 0u up to a maximum at 90u, which is approximately the Ernst angle for gray matter at 7 T for the used TR [20].
The EPI data were collected with the following parameters [18]: matrix = 64674, in-plane resolution = 3 mm63 mm, 49 sequentially acquired slices, slice thickness = 2 mm, interslice gap = 1 mm, TE = 30 ms, volume TR = 3.43 s, echo-spacing = 500 ms, BW = 2298 Hz/Px. At the beginning of the experiment the manufacturer's automatic adjustment procedure Figure 2. Validation of extended model fit for phantom data acquired at 3T using a 32-channel RF head receive coil. Results for a single central slice through phantom data show maps of A) tSNR, B) SNR9 0 C) k, D) SNR 0 E) SNR9 0 /k F) SNR9 0 /(k ?SNR 0 ) . The tSNR and SNR maps are shown for the EPI data acquired with the largest flip angle (90u). The k map was estimated by fitting the extended model to each voxel of the tSNR and SNR9 0 maps calculated from EPI time series acquired at 3 different flip angles (17u, 37u and 90u). In G, a plot of the mean and SD of tSNR against SNR9 0 (red squares), SNR 0 (green triangles) and SNR9 0 /k (blue circles) within a 20620 voxels ROI at the centre of the phantom is shown for all flip angle time series and the red solid line shows the fit of the extended model to the SNR9 0 values. The black dashed line is the line of identity. doi:10.1371/journal.pone.0052075.g002 was performed to correct for first and second order distortions in the static magnetic field. The k-space raw data from all time series were reconstructed using an algebraic trajectory-based reconstruction approach designed to minimize ghosting [19] followed by RSS calculation to combine the images from the 32 receiver channels. The resulting magnitude images were then processed using routines implemented in Matlab.
A measure of noise (s9 0 ) was estimated from a background ROI from the RSS images acquired without RF excitation (i.e. with 0u flip angle) using equation 4. According to [8], this particular noise measure is not affected by noise correlations. However, since in general this assumption can be violated, we use this noise measure to calculate apparent image SNR using SNR9 0 = I 32 /s9 0 where I 32 is the mean signal at each voxel in the RSS combined images. The tSNR was calculated using tSNR = I 32 /s t I32 where s t I32 is the temporal standard deviation at each voxel. Voxelwise maps of image SNR and tSNR were calculated for each flip angle time series after discarding the first 5 volumes and correcting for low frequency temporal drifts by fitting and removing a linear and quadratic function of image number. The extended model defined in equation 6 was fitted to the tSNR and SNR9 0 values at each voxel in the central slice of the phantom resulting in maps of l and k values.
The independent estimate of SNR 0 was based on an established method to estimate actual image SNR for RSS combined images using the individual channels of k-space raw data [3,5,6]. The SNR 0 at each voxel in the RSS combined image is given by.
where S is the vector of complex image values across all coils for a given voxel, S H is the Hermitian transpose of S, Y is the noise covariance matrix and b is a correction factor to account for the algebraic reconstruction method used to transform the k-space raw data to images and autocorrelation in k-space data as well as other effects, e.g. see [3,5]. The complex valued trajectory-based reconstructed data from each channel in the first non-discarded image of each flip angle EPI time series formed the vector S. Y was calculated from the covariance between all pairs of channels of k-space raw data acquired as the thermal noise measurement (i.e. acquired with no RF excitation). The correction factor b was estimated from the square root of the ratio between the integrated power in the k-space raw noise data and the trajectory-based reconstructed noise data. Note that this estimate for SNR 0 takes into account noise correlations between receiver channels. Resulting maps of SNR 0 were compared with maps of tSNR, SNR9 0 and k values estimated from the RSS combined EPI time series. The extended model was also fitted to the tSNR and SNR 0 values to estimate l and k.

Fit of Extended Model to Human task-free fMRI Data
A 7T whole-body MRI scanner (Siemens Healthcare) using a 24-channel receive head coil with dedicated CP coil for RF transmission (Nova Medical, Inc., Wilmington, MA) was used to acquire EPI time series from 5 subjects. Written informed consent was obtained from each participant for the study approved by the ethics committee of the University of Magdeburg. The data were acquired as part of a comprehensive study to investigate the impact of physiological noise correction at 7T (further details of the full study are reported in [4]).  Parameters are given for model fits to tSNR versus SNR9 0 (data shown in Figure 3a) and tSNR versus SNR 0 (data shown in Figure 3b). doi:10.1371/journal.pone.0052075.t001 For each subject, 5 EPI runs were acquired while subjects were presented with a blank screen and instructed to rest with their eyes open. Each run comprised 150 volumes and was acquired with one of the following flip angles: 8u, 16u, 26u, 38u and 70u, selected in a randomized order. These flip angles were selected to produce images with an equally spaced range of signal levels from 0u up to a maximum at 70u, the Ernst angle for grey matter at 7 T at a TR = 2 s, (assuming T1 = 1.9 s [20]). At higher field strengths it is important to account for large B1 inhomogeneities [21] when setting the RF excitation flip angle to maximize the image SNR in a particular region of interest. Here the scanner adjustment procedures resulted in an actual RF excitation flip angle, which was reasonably close (within ,10%) of the nominal flip angle in the studied region of interest, i.e., visual cortex. A thermal (raw) noise measurement of 20 EPI volumes with no RF excitation (i.e. 0u flip angle) was also acquired.
At the beginning of the experiment the manufacturer's automatic adjustment procedure, optimized for 7T, was performed to correct for first and second order distortions in the static magnetic field and to set the RF transmitter voltage. The EPI data were collected with the following parameters [18]: matrix = 64664, in-plane resolution = 3 mm63 mm, 40 sequentially acquired slices, slice thickness = 2 mm, interslice gap = 1 mm, TE = 25 ms, volume TR = 2 s, echo-spacing = 500 ms, BW = 2298 Hz/Px. The slice block was axial-to-coronal single oblique and aligned and centered manually with the calcarine fissure. For each subject, an axial, dual echo, gradient echo field map and a T1-weighted anatomical image (MPRAGE, resolution = (1 mm) 3 , TE = 3.72 ms, TR = 2000 ms, flip angle = 5u, TI = 1050 ms) were also acquired.
The k-space raw complex data from all time series were reconstructed using a trajectory-based reconstruction [19] followed by RSS calculation to combine the images from the 24 receiver channels. The resulting magnitude images were then processed using SPM8 (http://www.fil.ion.ucl.ac.uk/spm/, [22]) with additional routines implemented in Matlab. After discarding the first five volumes of each run, the EPI data were spatially coregistered to the first volume of the first run. The gradient echo field map was processed to create a voxel displacement map and used to correct the realigned images for geometric distortions [23]. Each subject's anatomical image was registered to their corre-sponding undistorted, realigned EPI data and segmented into grey and white matter tissue probability maps using the unified segmentation procedure in SPM8 [24]. The resulting inverse spatial normalization parameters were used to match a brain atlas (AAL toolbox, [25]) defined in MNI space [26] to the anatomy of each subject and hence define a visual cortex (VC) ROI. Finally the ROIs were restricted to grey matter by masking with the grey matter tissue probability maps (from the segmentation step) after thresholding at a probability value greater than 0.01.
For each subject the noise measure (s9 0 ) was estimated from the thermal noise measurement (i.e. with no RF excitation). For the time series acquired with the different (non-zero) flip angles, maps of SNR9 0 and tSNR were calculated as described previously for the phantom data. The mean SNR9 0 and tSNR were calculated for each subject's VC ROI and the mean and SD were estimated over subjects for each flip angle. The original noise model defined in equation 3 and the extended model defined in equation 6 were fitted to the mean tSNR and SNR9 0 values. Values for l and model fit sum-of-squares errors (SSE) were estimated for the two models and k was estimated for the extended model only.
An independent estimate of actual SNR 0 was also calculated for each subject and each flip angle time series using the complex trajectory-based reconstructed data from each channel in the first non-discarded image in the EPI time series and the noise covariance as described previously for the phantom data. The mean SNR 0 was calculated for each subject's VC ROI and the mean and SD were estimated over subjects for each flip angle. As for the SNR9 0 values, the original and extended noise models were fitted to the mean tSNR and SNR 0 values. Values for l and model fit sum-of-squares errors (SSE) were estimated for the two models and k was estimated for the extended model only.

Accuracy and Precision of Model Fit Estimated from Monte Carlo Simulations
The Monte Carlo simulations of the accuracy and precision indicated that the non-linear model fit is an unbiased estimator, provided that the set of five SNR9 0 values (over which the fitting is performed) are well distributed. Figure 1 shows the top 5% best sets of five SNR9 0 values for achieving model fits with a minimal bias for different SNR9 0 ranges (50-600 (a,b), 50-300 (c,d)) and k values (1.4 (a,c) 1.8 (b,d)). Approximately speaking, unbiased fits were observed for SNR9 0 sets which cover rather densely the low and high SNR9 0 range and less densely the mid SNR9 0 range. The bias in 1/l and k was always smaller than 1.2% when the best SNR9 0 sets were chosen, indicating a high accuracy.
For the set of best SNR9 0 values yielding minimal bias, the precision in estimating 1/l was high with SD always smaller than 7.0 for the large SNR9 0 range (50-600) and smaller than 11.3 for the small SNR9 0 range (50-300). These SD values compared reasonably with the overall lowest SD of 2.7 and 4.5 achieved in all 5000 simulations of different SNR9 0 , i.e., the highest achievable precision. The precision in estimating k was comparably lower, since the SD was always lower than 0.45 for the large SNR9 0 range and 0.27 for the small SNR9 0 range (compared to the minimal SD ranging between 0.12-0.15 for all simulations, i.e., the highest precision achievable).
For the gel phantom, Monte Carlo simulations demonstrated that it is possible to determine k robustly with high precision and accuracy even when only relatively low values for SNR9 0 /k were measured compared to 1/l, i.e., a maximal SNR9 0 /k = 180 compared to 1/l = 1800. The maximal bias in k was smaller than 2.3% and the SD was smaller than 0.085 for all simulated k values from 1.0 to 2.0. Note that 1/l values could not be reliably estimated as expected from the fact that only low SNR measurements informed the fit.

Validation of Extended Model Fit using Phantom Data and Independent Estimate of Actual Image SNR
The results of the validation using phantom data and an independent estimate of actual image SNR are shown in Figure 2. Maps of tSNR, SNR9 0 , and SNR 0 for the data acquired from the largest flip angle (90u) are shown in Figures 2a, b and d respectively and the map of k values estimated from the fit of the extended model to the SNR9 0 maps from all flip angles is shown in Figure 2c. In terms of magnitude and spatial distribution, the similarity between the tSNR and SNR 0 maps and the difference between those and the SNR9 0 map are apparent. However, voxelwise scaling of the SNR9 0 map by the k map (Figure 2c) resulted in a map (Figure 2e) that was comparable to both the tSNR and SNR 0 maps. This result validates the proposed linear relationship between SNR9 0 and SNR 0 (equation 5). The map of the ratio SNR9 0 /(k?SNR 0 ) shown in Figure 2f confirmed that the k value estimated from the phantom data satisfied the equation k = SNR9 0 /SNR 0 to within 4% averaged over the signal in the central slice of the phantom (mean 6 SD of the ratio SNR9 0 /(k?SNR 0 ) = 1.0460.1). The map of k values ranged from a minimum of 0.4 up to 1.7 in the central slice of the phantom. Since SNR9 0 is affected by noise correlations whereas SNR 0 accounts for them, this result demonstrates the spatial variation of the impact of correlated noise. Figure 2g shows a plot of tSNR against SNR9 0 , SNR 0 and SNR9 0 /k (mean 6 SD) within a 20620 voxels ROI at the centre of the phantom for all flip angle time series. For this ROI, the fit of the extended model to the tSNR and SNR9 0 values resulted in k = 1.560.1 (mean 6 SD) and the fit to the tSNR and SNR 0 values resulted in k = 1.060.04 (mean 6 SD), since noise correlations were accounted for in SNR 0 . As expected from the results of the Monte Carlo simulation, the fit of the l value was imprecise and biased with l = 1/935.861/833.3 for the first fit and l = 1/913.661/813.3 for the latter fit (mean 6 SD). Figure 3a shows plots of mean tSNR versus mean SNR9 0 measured in the VC ROI of 5 subjects for the task-free fMRI data acquired at 7T. The parameters estimated from the fits of the two models are given in Table 1. The extended model using equation 6 (red solid curve) gave a much improved fit compared to the original model using equation 3 (blue dashed curve), with SSE = 17.8 and 191.5 respectively. This was particularly apparent at low to medium SNR9 0 values and was reflected by the value estimated for the parameter k = 1.560.2 (mean 6 SD). Figure 4 shows the subject specific EPI images acquired at the 70u flip angle together with the maps of k estimated for the fit of the extended model to the tSNR versus SNR9 0 values at each voxel. The value estimated for 1/l was slightly higher for the extended model compared to the original model (95.269.2 and 89.169.4 respectively, (mean 6 SD)). Note that the difference between these values fell within the precision in estimating 1/l for this SNR9 0 range as predicted by the simulation results (i.e. SD ,7.0). Figure 3b shows plots of tSNR against the actual image SNR (SNR 0 ) for the same data as in Figure 3a. The fits of the extended model (solid green line) and the original model (blue dashed line) are shown for these data and the parameters are given in Table 1. For these results, the fit of the extended model was still improved compared to the original model (SSE = 17.6 and 68.8 respectively), even though one may have expected the original model to describe the data as well as the extended model since actual image SNR values were used. The plots show that again the fit of the original model was more problematic at low to medium SNR 0 values. This was also highlighted by the fact that the k value estimated from the extended model fit was not 1.0 as expected, but 1.260.1 (mean 6 SD). However, it should also be noted that the difference in estimates of 1/l for the two models was not greater than their SD values.

Discussion
The model introduced by Krueger and Glover [1] describes the tSNR of an EPI time series as a function of image SNR. It can be used to estimate the maximum attainable tSNR (1/l) or BOLD sensitivity for a given set of MR acquisition parameters [2] and image processing strategies [4]. This model requires that the image SNR be estimated accurately, which is not straightforward for data acquired using multi-channel coils. First of all, access to the kspace raw data is necessary which can be cumbersome and often practically impossible for large numbers of receiver coils (e.g. up to 30 GB per EPI time series in this study). Secondly, the steps required to correctly scale the image SNR measurements are complicated, and may require corrections for the RSS combination and magnitude operation, as well as for noise correlations between receiver coils and the algebraic reconstruction. The complexity and practical challenge are reflected by extensive discussion in the literature [3,5,6,8,10,11,17]. In this work we have proposed an extension to this model that allows the tSNR of an EPI time series acquired using multi-channel receiver coils to be characterized by estimating an apparent image SNR value (SNR9 0 ) simply and directly from non-accelerated RSS combined images. The extended model includes an additional parameter k, which scales the apparent image SNR to be equal to the actual SNR (SNR 0 ). Since the apparent SNR9 0 does not account for noise correlations between receiver coils whereas the actual SNR 0 accounts for them, the parameter k reflects the impact of correlation on the data.

Robustness of the Extended Model Studied with Monte Carlo Simulations
The Monte Carlo simulations showed that the extended model parameters l and k can be accurately determined without significant bias from a series of five SNR9 0 and tSNR measurements when the set of SNR9 0 values is appropriately chosen. In other words, tSNR should be measured at high and low SNR9 0 values with less emphasis on the medium SNR9 0 range (Figure 1). Based on only 5 measurements, a reasonable precision for the estimated 1/l is achieved with a coefficient of variation of ,5-10%. However, the precision of k estimates is somewhat low with coefficients of variation between 20-30%. This is less critical when the extended model is used for studies on the impact of physiological noise whose primary outcome measure is 1/l. A higher precision in k could be achieved by increasing the sampling of SNR9 0 values. Furthermore, a robust estimation of k (coefficient of variation ,5%) can be achieved even when the maximal SNR9 0 is low compared to 1/l (i.e. when l approaches 0 as it can for phantom data when the MRI scanner is stable). In this case, estimates of 1/l are not reliable. For typical SNR 0 values observed in vivo at 7T, the interdependence between estimates of 1/l and k was low, corroborating a good estimability of the model.
We would like to note that the original model proposed by Krueger and Glover [1] shows similar characteristics in estimating 1/l and also relies on a sufficient range of SNR9 0 values to avoid bias in its estimation.

Validation of Extended Model using Phantom Data
The fit of the extended model to the phantom data resulted in a map of k values that accurately (within 4%) scaled the map of apparent SNR values (SNR9 0 ) to be equal to an independent measure of the actual image SNR (SNR 0 ) calculated using an established method [3,5,6]. The results showed that k estimated for SNR9 0 was spatially highly variable (.50% change) making the extended model essential for bias free estimates. The spatial variability of k can be well explained by the spatial variation of noise in composite images as a result of noise correlations between the different channels of the receiver system [13,14]. When the extended model was fitted to the estimate of actual image SNR, the estimate for k was as expected 1.060.01 (mean 6 SD). Note that noise correlations between the different channels of the receiver system were taken into account for the estimate of SNR 0 via the noise covariance between all pairs of receiver channels (i.e. equation 8). The large variability of values estimated for 1/l for the two data sets could be explained by the results of the Monte Carlo simulations, which demonstrated that 1/l values could not be reliably estimated when all measured SNR values are low compared to 1/l.

Fit of Extended Model to Human Task-free fMRI Data
The results for the human data showed a much improved fit of the extended model to tSNR and SNR9 0 measured from nonaccelerated RSS combined images, compared to the original model by Krueger and Glover [1]. Although the estimate of 1/l was very similar for the two models (within the SD values), the fit of the extended model was particularly improved for low to medium values of SNR9 0. This suggests that at low to medium SNR9 0 values, the original model may lead to incorrect estimates of achievable tSNR if SNR values are not measured accurately. This is particularly important for optimizing EPI parameters for this lower SNR range where thermal noise components contribute significantly, e.g. by varying spatial resolution [2,27] or flip angle [28]. Using the extended model to characterize the tSNR as a function of apparent image SNR (SNR9 0 ) estimated from RSS combined images avoids these potential problems.
Notably, when actual SNR values (SNR 0 ) were calculated using the established method, the extended model still resulted in a better fit than the original model and the estimate for k was slightly higher than expected (1.260.1, mean 6 SD versus 1.0). Possible reasons for this discrepancy are discussed below under methodological considerations. Interestingly, the variation in estimates of 1/l for both SNR measures and both models was within the given precision. This highlights the relatively high accuracy of estimates of 1/l for both the original and extended models, suggesting that in general these models are rather robust for estimating the maximum attainable tSNR for a given set of MR acquisition parameters.

Methodological Considerations
The validation of the extended model using phantom data clearly demonstrated that the estimate of k satisfied the equation k = SNR9 0 /SNR 0 when SNR 0 was calculated using an established independent measure. This was less clear-cut for the human data. A discrepancy of between 10% and 30% was observed between SNR9 0 /k and the independently measured SNR 0 , possibly due to the limited precision of the k estimates. The discrepancy might be reduced by sampling more SNR9 0 values. Furthermore, k was estimated for tSNR and SNR9 0 values averaged over an ROI for each subject, which could further decrease its precision due to the presence of outliers or effects of non-linearity (i.e. averaging before the non-linear fit). Examples of these are apparent in the estimated maps of k shown in Figure 4.
The discrepancy may also have been related to limitations in both the original and the extended noise models. For example, the original model assumes that tSNR < SNR 0 at lower SNR 0 values, (i.e. where noise is dominated by thermal noise). However, the estimation of tSNR in human data requires motion correction, which may alter the relationship between tSNR and SNR 0 . Furthermore, correlations between the physiological noise from different channels may have an effect on the tSNR, e.g. by spatial resampling, which is not included in the original model presented here and is also not accounted for by the image reconstruction. These two latter effects may influence the estimate of k from human EPI time series data in the extended model. As shown in Figure 4, estimates for k are much higher around the edges of the brain, and in the ventricles and other fluid spaces, where pulsatile motion from fluid and head movement are typically most apparent.
In this work, the parameter k has been interpreted as a factor which accounts for the effect of noise correlations by scaling the apparent image SNR (SNR9 0 ) so that it equals the actual SNR (SNR 0 ). In theory, k could be predicted from the receiver coils' flux lines and noise correlation coefficients. Practically this is complicated by large numbers of receiver coils, their geometric arrangement, electrical coupling properties and differing sensitivity profiles [8,13,14].
Many methods have been proposed to estimate the noise and hence the SNR in MR images (for example see [5] for a detailed review of methods). The noise measurement method used in this work (i.e. using equation 4) was selected because it can be estimated from a short thermal noise reference scan and has been shown to be independent of noise correlations [8]. Furthermore, it does not require any pulse sequence or image reconstruction modifications. It should be pointed out that this noise measure and the proposed extended model is applicable to all situations where there is a constant scale factor between the measured apparent SNR9 0 and the actual SNR 0 [5], although the spatial variation of k is likely to be different for different situations. For example, the extended model could be used to characterize the tSNR and physiological noise properties of GRAPPA parallel imaging acquisitions [7] and reconstructions because g-factors and undersampling factors introduce a constant scaling factor between apparent and actual SNR. The extended model in its current form is not valid for low SNR values (,50), if 32 or more channels are used since correction factors must be applied to account for changes in the noise distribution [8]. However, in principle a correction for low SNR values could be added to the model [3,8], although an iterative estimation approach may be necessary.

Conclusion
We have presented a simple extension to the noise model of [1] which allows the characterization of tSNR and SNR of nonaccelerated RSS combined images from multi-channel receiver coils without needing to handle raw data or to perform arduous custom-made image reconstruction. To our knowledge, no other published methods can estimate the model proposed in [1] from RSS combined magnitude data from multi-channel receiver coils.
We derive the proposed extended model from basic physical principles and show that our simple SNR correction is valid for non-accelerated data with an SNR .50 and for multi-channel receiver coils with #32 channels. Given these constraints, the model can be applied to all situations where a constant scaling factor exists between the actual SNR 0 and the apparent SNR9 0 measured in the reconstructed and combined images (e.g. nonaccelerated, RSS combined images). We successfully validated the model by Monte Carlo simulations as well as application to phantom data at 3T and task-free human fMRI data at 7T. In particular, the extended model improves the prediction of low to medium tSNR values, which is especially important for optimizing high-resolution fMRI experiments, where thermal noise components contribute significantly. In summary, the extended model offers a simple and robust method to investigate the fundamentals of physiological noise, to assess the impact of physiological noise correction and image reconstruction and to optimize EPI acquisition parameters.