Introduction to a novel T2 relaxation analysis method SAME-ECOS: Spectrum Analysis for Multiple Exponentials via Experimental Condition Oriented Simulation
Hanwen Liu, Qing-San Xiang, Roger Tam, Piotr Kozlowski, David K.B. Li, Alex L. Mackay, John K. Kramer, Cornelia Laule
FF U L L P A P E R
S u b m i t t e d t o M a g n e t i c R e s o n a n c e i n M e d i c i n e
Introduction to a novel T relaxation analysismethod SAME-ECOS: Spectrum Analysis forMultiple Exponentials via Experimental ConditionOriented Simulation Hanwen Liu ∗ | Qing-San Xiang | Roger Tam |Piotr Kozlowski | David K.B. Li | Alex L. MacKay | John K. Kramer | Cornelia Laule Physics & Astronomy, University of BritishColumbia International Collaboration on RepairDiscoveries Radiology, University of British Columbia Biomedical Engineering, University ofBritish Columbia Kinesiology, University of British Columbia Pathology, University of British Columbia
Correspondence
818 W 10th Ave, Vancouver, BC, V5Z1M9Email: [email protected]
Present address ∗
818 W 10th Ave, Vancouver, BC, V5Z1M9
Funding information
Funding support was provided by theMultiple Sclerosis Society of Canada,Natural Sciences and Engineering ResearchCouncil Discovery Grant.
We propose a novel T relaxation data analysis method whi-ch we have named spectrum analysis for multiple exponen-tials via experimental condition oriented simulation (SAME-ECOS). SAME-ECOS, which was developed based on a com-bination of information theory and machine learning neu-ral network algorithms, is tailored for different MR exper-imental conditions, decomposing multi-exponential decaydata into T spectra, which had been considered an ill-posedproblem using conventional fitting algorithms, including thecommonly used non-negative least squares (NNLS) method.Our results demonstrated that, compared with NNLS, thesimulation-derived SAME-ECOS model yields much morereliable T spectra in a dramatically shorter time, increas-ing the feasibility of multi-component T decay analysis inclinical settings. K E Y W O R D S
SAME-ECOS, T relaxation, non-negative least squares, resolutionlimit, myelin water imaging, machine learning a r X i v : . [ phy s i c s . m e d - ph ] S e p L IU ET AL ..
SAME-ECOS, T relaxation, non-negative least squares, resolutionlimit, myelin water imaging, machine learning a r X i v : . [ phy s i c s . m e d - ph ] S e p L IU ET AL .. | INTRODUCTION T relaxation in biological tissues measured with a multi-echo experiment is typically characterized by multi-exponenti-al decays because multiple water pools may exist within a single image voxel. [1] The T times of different water poolsare governed by the microenvironment of water molecules. For example, myelin water, the water trapped in myelinbilayers, exhibits a shorter T relaxation time than that of intra/extra-cellular water (IE) and free water. [2] However,accurately depicting the spectrum of constituent T components for each image voxel from MR relaxation data isnontrivial. Because the decay process of multiple water pools takes place simultaneously, the MR receiver coil can onlyrecord a signal that is the sum of multiple exponential decay components. Consequently, one has to solve a mathematicalproblem of fitting a superimposed relaxation signal into its constituent components. This mathematically complexproblem is commonly seen in many other quantitative sciences and is often considered as an ill-posed problem. [3]The analytical and numerical solutions to this problem have been comprehensively reviewed by Istratov et al. from amathematical perspective. [3]To provide a suitable solution specific to MR data, Whittall and Mackay [4] introduced the non-negative leastsquares (NNLS) [5, 6, 7] method that decomposes the multi-echo decay data into a spectrum of positive T times. TheNNLS method makes no prior assumptions about the total number of T components, which is a desirable feature formodeling complex biological tissues with heterogeneous compositions. However, without constraining the number ofcomponents, the T fitting problem becomes underdetermined with non-unique solutions, making NNLS unstable andhighly susceptible to noise even with strong regularization. In contrast, other fitting methods such as the quasi-Newtonalgorithm by Du et al. [8] and the Wald distribution by Akhondi-Asl et al. [9] usually provide more stable results andbetter noise resistance, but at the expense of modelling only two or three water pools, which limits their usage onunknown pathological tissues.Simultaneously increasing the complexity and stability of a model is a paradox, seemingly impossible to improveone without diminishing the other in the game of multi-exponential decomposition. Our quest now was to find the lineof best balance between these two factors. According to studies of information theory, the decay components canonly be resolved to a certain resolution limit at a given signal to noise ratio (SNR). [10, 11] That is to say, due to noisecontamination, there is always a limit in how closely the two neighboring components can be resolved, regardless of howsophisticated the analysis method is. This fundamental restriction on the resolution limit leads to a correlation betweenthe SNR and the maximum number of components that can be reliably resolved in a particular analysis range. [12] Theexact expression of this correlation is presented in the Methods, and plays a crucial role in our proposed approach.On the other hand, machine learning algorithms, in particular supervised neural network methods [13], have beensuccessfully implemented in many MR applications [14], especially for tasks involving parameter estimation. [15, 16]In short, a neural network can be trained to discover hidden patterns in data and to learn the mapping between twovector spaces. A trained neural network usually outperforms most conventional methods in terms of better accuracyand faster speed, particularly if the mapping is highly nonlinear. Additional background on the approximation propertiesof neural networks can be found elsewhere. [17]Based on information theory and neural network algorithm approaches, we propose a novel method which we havecalled S pectrum A nalysis for M ulti- E xponentials via E xperimental C ondition O riented S imulation ( SAME-ECOS ) forthe analysis of multi-echo T relaxation data. The general concept of SAME-ECOS ( Figure 1 ) can be briefly described bya series of calculation, simulation and model training operations: (1) determine the T range and resolution limit basedon the experimental conditions such as SNR and echo times; (2) generate sufficient examples of random T spectrawithin the T range obeying the resolution limit; (3) compute multi-echo decay data using the randomly generated T spectra; (4) train a neural network model to learn the mapping between the simulated multi-echo decay data and the IU ET AL . 3
F I G U R E 1
SAME-ECOS workflow.ground truth spectrum labels; (5) apply the trained neural network model to experimental data. In general, SAME-ECOSis a simulation-derived solver, powered by information theory and machine learning, to the problem of fitting multi-exponential decay data into a T spectrum. It is worth highlighting that SAME-ECOS has high flexibility of tailoring itselfto different experimental conditions attributed to its unique simulation workflow. The detailed SAME-ECOS algorithmis demonstrated in the Methods part by presenting an example with in-depth explanations. | METHODS
Because the SAME-ECOS analysis method is tailored to different MR experimental conditions, we use one specificexample here as a paradigm to demonstrate its simulation workflow, trained model evaluation, and experimental dataapplication. Therefore, the Methods part of this introductory paper to SAME-ECOS is presented in the followingsections as a particular experiment. | In-vivo MRI experiment ∆ TE/TR=10/10/1000ms, refocusing flip angle (FA)=180 ◦ , axialmatrix size = 232 × × × for 20 slices, reconstructed resolution = 1 × × for40 slices, acquisition time = 14.4 minutes) [18] from one healthy volunteer (male, 57 years old) was collected at a 3Tscanner (Philips Achieva) using an 8-channel head coil. | SAME-ECOS workflow (1) Choose the SNR range to be 70-300.
To initialize the simulation workflow, the SNR of the experimental data needsto be determined first, as calculations of the T range, the resolution limit, and especially the noise simulation, alldepend on the SNR. The SNR was initially estimated to be approximately 167 by examining the noise variance of the air L IU ET AL ..
To initialize the simulation workflow, the SNR of the experimental data needsto be determined first, as calculations of the T range, the resolution limit, and especially the noise simulation, alldepend on the SNR. The SNR was initially estimated to be approximately 167 by examining the noise variance of the air L IU ET AL .. voxels proximate to the skull on the 1 st echo image. However, most of today’s MRI images acquired by multi-channelcoils are reconstructed using parallel imaging [19], which complicates the estimation of SNR. Other factors, such asB1 inhomogeneity, may lead to regional variations in SNR. Therefore, the SNR of experimental data cannot be simplyassessed by a single definitive number. To accommodate these hurdles, we empirically assigned the SNR a wide range of70 to 300, instead of using a single value. This approach allows us to randomly select any SNR within the designatedrange at each simulation realization, making the resulting simulated data ‘all-inclusive’ after many realizations. (2) Define the T range to be 7-2000ms. Intuitively, the shortest detectable T decay component should haveits residual signal greater than the noise level at the first measurement. Because the first echo is used later as anormalization factor, the second echo is actually regarded as the first measurement in our analysis. In theory, a minimumof two measurement points would be needed for the purpose of T fitting. That means the residual signal of the shortestT component in our analysis should be higher than the noise level at the third echo. Thus, the lower bound of the T range can thus be obtained using equation 1. T min = − r d echo ti me ln (cid:16) SN R (cid:17) (1)Given the SNR range of 70 to 300, the lower boundary is calculated to be approximately 7ms. On the other hand, anideal experimental condition would be monitoring the decay as long as possible until the longest T component decayscompletely. [20, 21, 22] Then the upper bound for the analyzable T range can be determined using equation 2. [3] T max = − l ast echo ti me ln (cid:16) SN R (cid:17) (2)However, the ideal experimental condition rarely happens as the decay monitoring time is often compromised to achievea shorter scanning time for most in-vivo MR experiments including our own (last TE = 320ms). If equation 2 were usedwith our SNR range of 70 to 300, the longest analyzable T component would be less than 80ms, which is considered tobe too short for the analysis of in-vivo brain imaging. Therefore, the upper boundary is manually extended to 2000msbased on the literature T ranges for brain. [1, 2, 7] (3). Determine the maximum number of resolvable T components to be M = 5 and the resolution limit δ = 3.098. For a given SNR and T range, the decay components can only be resolved to a certain resolution limit. Link et al. [12]derived an expression (equation 3) relating T min , T max and SNR, to the M resolvable exponentials. M ln (cid:18) T max T min (cid:19) × sinh (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) π × M ln (cid:18) T max T min (cid:19) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) = (cid:18) SN RM (cid:19) (3)Derivations and justifications of equation 3 can be found in several publications. [3, 11, 12] Based on these previoustheoretical studies, the maximum number of resolvable T components M = range. Note that M = resolution limit δ =3.098 was determined by equation 4. [3] δ = (cid:32) T max T min (cid:33) M (4) IU ET AL . 5 (4). For each simulation realization:a. generate random integers n < M, FA ∈ [90 ◦ , 180 ◦ ], and SNR ∈ [70,300]. Note that the simulation did not includethe case of n = M because when n = M, there is only one possible configuration for the M T locations within the T range obeying the resolution limit, which would introduce an unwanted bias into the simulated dataset. Also note thatthe actual FA can deviate substantially from the prescribed FA due to B1 inhomogeneity, so we gave a 90 ◦ tolerance toaccount for the FA variations. b. generate n randomly located T components within the T range obeying the resolution limit δ . The locationsand amplitudes of the n T components were randomly assigned and normalized to one. c. synthesize 32-echo decay data. The pure decay signal without noise of the n T components (denoted as S pure )were synthesized (equation 5) using the extended phase graph (EPG) algorithm [23, 24] with T = 2000ms as a default. S pure = i = n (cid:213) i =1 Amplitude i × E P G (cid:0) T , i , selected FA (cid:1) (5)To mimic the noise profile of a real MRI image [25], the S pure was first projected into the real and imaginary axes by arandom phase factor θ ∈ [0 ◦ , 90 ◦ ], followed by adding noises on both axes, and finally producing the magnitude of thenoisy signal (equation 6). S noisy = (cid:113)(cid:0) S pure × sin θ + noi se − (cid:1) + (cid:0) S pure × cos θ + noi se − (cid:1) (6)where noi se _ and noi se _ were independently sampled from a Gaussian distribution with its mean = = 1 selected SNR × (cid:112) π / (7)such that the noisy signal S noisy would follow a Rician distribution at the selected SNR level. The synthesized noisydecay data were subsequently normalized to the 1 st echo and saved for model training. d. Embed n T components into a spectrum representation depicted by 40 basis T s . The 40 basis T s (t , t t ...t ) were equally spaced within the T range on a logarithmic scale. The 40 weighting factors of the basis T s wereused to represent the spectrum of n T components (T , T , ... T ). Explicitly, each T component was depicted as aGaussian-shaped peak by the basis T s, with the weighting factor w i of the i th basis T (t i ) being calculated as w i = n (cid:213) j =1 √ π e − (cid:16) t i − T , j (cid:17) (8)The embedded T2 basis representation was normalized to one and recorded as the ground truth spectra for modeltraining. (5) Train a neural network to map the decay data to its T spectrum. A neural network (hidden layers: × × × × , activation: SeLU [26] (hidden layers) and softmax [27] (output layer), optimizer: Adamax [28], loss:categorical cross-entropy [29]) was constructed using TensorFlow [30] to take 32-echo decay data as input and predictthe weighting factors of the 40 basis T s at the output layer. The constructed neural network was trained to map thedecay data to its T spectrum. We yielded 3,000,000 simulation realizations (step 4), 90% of which were used for theneural network training. The remaining 10% simulation realizations were used for the validation that determined the L IU ET AL ..
To initialize the simulation workflow, the SNR of the experimental data needsto be determined first, as calculations of the T range, the resolution limit, and especially the noise simulation, alldepend on the SNR. The SNR was initially estimated to be approximately 167 by examining the noise variance of the air L IU ET AL .. voxels proximate to the skull on the 1 st echo image. However, most of today’s MRI images acquired by multi-channelcoils are reconstructed using parallel imaging [19], which complicates the estimation of SNR. Other factors, such asB1 inhomogeneity, may lead to regional variations in SNR. Therefore, the SNR of experimental data cannot be simplyassessed by a single definitive number. To accommodate these hurdles, we empirically assigned the SNR a wide range of70 to 300, instead of using a single value. This approach allows us to randomly select any SNR within the designatedrange at each simulation realization, making the resulting simulated data ‘all-inclusive’ after many realizations. (2) Define the T range to be 7-2000ms. Intuitively, the shortest detectable T decay component should haveits residual signal greater than the noise level at the first measurement. Because the first echo is used later as anormalization factor, the second echo is actually regarded as the first measurement in our analysis. In theory, a minimumof two measurement points would be needed for the purpose of T fitting. That means the residual signal of the shortestT component in our analysis should be higher than the noise level at the third echo. Thus, the lower bound of the T range can thus be obtained using equation 1. T min = − r d echo ti me ln (cid:16) SN R (cid:17) (1)Given the SNR range of 70 to 300, the lower boundary is calculated to be approximately 7ms. On the other hand, anideal experimental condition would be monitoring the decay as long as possible until the longest T component decayscompletely. [20, 21, 22] Then the upper bound for the analyzable T range can be determined using equation 2. [3] T max = − l ast echo ti me ln (cid:16) SN R (cid:17) (2)However, the ideal experimental condition rarely happens as the decay monitoring time is often compromised to achievea shorter scanning time for most in-vivo MR experiments including our own (last TE = 320ms). If equation 2 were usedwith our SNR range of 70 to 300, the longest analyzable T component would be less than 80ms, which is considered tobe too short for the analysis of in-vivo brain imaging. Therefore, the upper boundary is manually extended to 2000msbased on the literature T ranges for brain. [1, 2, 7] (3). Determine the maximum number of resolvable T components to be M = 5 and the resolution limit δ = 3.098. For a given SNR and T range, the decay components can only be resolved to a certain resolution limit. Link et al. [12]derived an expression (equation 3) relating T min , T max and SNR, to the M resolvable exponentials. M ln (cid:18) T max T min (cid:19) × sinh (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) π × M ln (cid:18) T max T min (cid:19) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) = (cid:18) SN RM (cid:19) (3)Derivations and justifications of equation 3 can be found in several publications. [3, 11, 12] Based on these previoustheoretical studies, the maximum number of resolvable T components M = range. Note that M = resolution limit δ =3.098 was determined by equation 4. [3] δ = (cid:32) T max T min (cid:33) M (4) IU ET AL . 5 (4). For each simulation realization:a. generate random integers n < M, FA ∈ [90 ◦ , 180 ◦ ], and SNR ∈ [70,300]. Note that the simulation did not includethe case of n = M because when n = M, there is only one possible configuration for the M T locations within the T range obeying the resolution limit, which would introduce an unwanted bias into the simulated dataset. Also note thatthe actual FA can deviate substantially from the prescribed FA due to B1 inhomogeneity, so we gave a 90 ◦ tolerance toaccount for the FA variations. b. generate n randomly located T components within the T range obeying the resolution limit δ . The locationsand amplitudes of the n T components were randomly assigned and normalized to one. c. synthesize 32-echo decay data. The pure decay signal without noise of the n T components (denoted as S pure )were synthesized (equation 5) using the extended phase graph (EPG) algorithm [23, 24] with T = 2000ms as a default. S pure = i = n (cid:213) i =1 Amplitude i × E P G (cid:0) T , i , selected FA (cid:1) (5)To mimic the noise profile of a real MRI image [25], the S pure was first projected into the real and imaginary axes by arandom phase factor θ ∈ [0 ◦ , 90 ◦ ], followed by adding noises on both axes, and finally producing the magnitude of thenoisy signal (equation 6). S noisy = (cid:113)(cid:0) S pure × sin θ + noi se − (cid:1) + (cid:0) S pure × cos θ + noi se − (cid:1) (6)where noi se _ and noi se _ were independently sampled from a Gaussian distribution with its mean = = 1 selected SNR × (cid:112) π / (7)such that the noisy signal S noisy would follow a Rician distribution at the selected SNR level. The synthesized noisydecay data were subsequently normalized to the 1 st echo and saved for model training. d. Embed n T components into a spectrum representation depicted by 40 basis T s . The 40 basis T s (t , t t ...t ) were equally spaced within the T range on a logarithmic scale. The 40 weighting factors of the basis T s wereused to represent the spectrum of n T components (T , T , ... T ). Explicitly, each T component was depicted as aGaussian-shaped peak by the basis T s, with the weighting factor w i of the i th basis T (t i ) being calculated as w i = n (cid:213) j =1 √ π e − (cid:16) t i − T , j (cid:17) (8)The embedded T2 basis representation was normalized to one and recorded as the ground truth spectra for modeltraining. (5) Train a neural network to map the decay data to its T spectrum. A neural network (hidden layers: × × × × , activation: SeLU [26] (hidden layers) and softmax [27] (output layer), optimizer: Adamax [28], loss:categorical cross-entropy [29]) was constructed using TensorFlow [30] to take 32-echo decay data as input and predictthe weighting factors of the 40 basis T s at the output layer. The constructed neural network was trained to map thedecay data to its T spectrum. We yielded 3,000,000 simulation realizations (step 4), 90% of which were used for theneural network training. The remaining 10% simulation realizations were used for the validation that determined the L IU ET AL .. Pre-defined spectra T locations (ms) T amplitudes (normalized) Spectrum 1 100 1Spectrum 2 25, 120 0.3, 0.7Spectrum 3 15, 80, 50 0.3, 0.5, 0.2Spectrum 4 10, 60, 300, 1200 0.2, 0.4, 0.3, 0.1
TA B L E 1
Four pre-defined ground truth spectra.
The T locations of each spectrum are selected obeying theresolution limit. Amplitudes are normalized to one.stopping criterion for the training process. The training was stopped when the accuracy on the validation set did notimprove further. This particular trained neural network is denoted as the SAME-ECOS model hereafter and can beapplied to new 32-echo decay data to obtain T spectrum. | SAME-ECOS model performance evaluation
The performance of the SAME-ECOS model was evaluated using three designed tests and compared respectivelywith the results determined by a regularized NNLS solver equipped with stimulated echo correction [24] (analysisprogram can be requested here: https://mriresearch.med.ubc.ca/news-projects/myelin-water-fraction/ ).The kernel matrix for the NNLS analysis was adjusted accordingly to match our experimental and simulation parameters.The regularization parameter was chosen to be the largest value that allows a misfit of less than 1.02 times of theminimum misfit, which is commonly used in many studies. [31, 32, 33, 34, 35]
Test 1: spectra and their noisy 32-echo decay data were randomly generated followingthe workflow described in section 2.2. The decay data were analyzed by the SAME-ECOS model and NNLS respectivelyto produce the T spectra. The processing time was recorded. The ‘goodness’ of each spectrum fitting by both methodswas quantitatively assessed using cosine similarity scores [36], which report values between 0 to 1, with 0 being theleast similar and 1 being the most similar to the ground truth labels. The cosine similarity score is a commonly usedmetric that measures the similarities between two vectors, especially when the vectors are high dimensional. It is asuitable metric for our task since each T spectrum can be treated as a vector of 40 dimensions. The calculation of thecosine similarity score is defined explicitly in the following formula cosi ne si mi l ar i t y scor e = X • Y (cid:107) X (cid:107) × (cid:107) Y (cid:107) (9)Where X and Y are the vector representations of the predicted and the ground truth spectra; (cid:107) X (cid:107) and (cid:107) Y (cid:107) are theirEuclidean norms respectively. Paired t-test was performed to determine whether there was a significant difference(P<0.05) in the cosine similarity scores calculated by SAME-ECOS and NNLS. Test 2:
To examine the model robustness to noise, simulated decay data of 4 pre-defined ground truth T spectra( Table 1 ) at SNR = 100 and FA = 180 ◦ , each with 100 different noise realizations, were passed to the SAME-ECOS modeland NNLS for spectrum predictions. The location and amplitude of each ground truth spectrum were manually chosento provide a visual-friendly data presentation. The spectrum analysis results of both SAME-ECOS and NNLS werenormalized (sum to unity) prior to comparison. The similarity between each predicted and ground truth spectrum wasassessed by cosine similarity scores defined above. The mean and standard deviation of the cosine similarity scoreswere also calculated. IU ET AL . 7
Test 3: spectra and their 32-echo decay data with noise realizations were randomlygenerated according to steps described in section 2.2. The decay data were analyzed by the SAME-ECOS model andNNLS respectively to produce the T spectra. The myelin water fraction (MWF, a fraction of signal with T s < 40ms) was extracted from each predicted spectrum and compared with the ground truth MWF. Mean absolute error (MAE) inthe MWF estimation was computed. The correlation between the errors of MWF estimation and the FA was evaluatedusing Pearson correlation analysis. | Apply SAME-ECOS to experimental data
The SAME-ECOS model was applied to the experimental in-vivo GRASE data, which were pre-processed by normalizingto the first echo image. The processing time for the whole brain data analysis was recorded. From the resulting T spectra, the MWF (T s < 40ms) was extracted for each voxel. The masks of regions of interest (ROI) including wholebrain, whole white matter, corpus callosum, corticospinal tract, forceps major, and forceps minor were produced usingthe first echo image via the FSL segmentation tool. [37] The T spectra and MWF map was also produced using NNLS inthe T range of 7-2000ms as a reference for comparison. | RESULTS3.1 | Performance evaluation via simulation tests3.1.1 | Test 1: general performance
The processing time to make 300,000 spectra predictions was 22 seconds for the SAME-ECOS model and 1,620 secondsfor NNLS (CPU: Intel(R) Core(TM) i7-5930K @ 3.5 GHz, 32 GB RAM). The mean cosine similarity score of the 300,000spectra predicted by the SAME-ECOS model (0.838 ± ± | Test 2: robustness to noise
The resulting SAME-ECOS and NNLS spectra from each individual noise realization (grey) were plotted in
Figure 2 tocompare against the ground truth spectrum (green). The average spectrum from 100 different noise realizations wasalso calculated for all SAME-ECOS and NNLS scenarios (blue: SAME-ECOS, red: NNLS). SAME-ECOS produced visuallybetter results than NNLS for all scenarios. For spectra consisting of one or two T components, the SAME-ECOS modelwas able to make almost perfect predictions (cosine similarity score: 0.998 ± ± ± ± components, the performances of bothmethods started to degrade. However, SAME-ECOS (cosine similarity score: 0.681 ± ± ± ± | Test 3: MWF prediction accuracy
MWF values were extracted using both SAME-ECOS and NNLS methods and plotted against 10,000 ground truth MWFvalues in
Figure 3 . SAME-ECOS MWF (blue, MAE=0.050) demonstrated slightly better agreement with the ground L IU ET AL ..
Figure 3 . SAME-ECOS MWF (blue, MAE=0.050) demonstrated slightly better agreement with the ground L IU ET AL .. F I G U R E 2 T spectra produced by SAME-ECOS (blue) and NNLS (red) respectively are compared with theground truth spectra (green). Simulated decay data, which are generated from 4 pre-defined ground truth spectra(
Table 1 ) each with 100 different noise realizations, are fed into the trained SAME-ECOS model and NNLS algorithm togenerate the T spectra. Faded gray lines indicate the produced spectra for each noise realization. The mean andstandard deviation of cosine similarity scores of 100 realizations are shown for each sub-figure. IU ET AL . 9
F I G U R E 3
Myelin water fraction (MWF) produced by SAME-ECOS (blue) and NNLS (red) are compared withthe ground truth MWF (green). s < 40ms) was extracted from each predicted spectrum. IU ET AL ..
Myelin water fraction (MWF) produced by SAME-ECOS (blue) and NNLS (red) are compared withthe ground truth MWF (green). s < 40ms) was extracted from each predicted spectrum. IU ET AL .. F I G U R E 4
Flip angle (FA) dependence of myelin water fraction (MWF) prediction errors. s <40ms) was calculated and compared with FA. The mean error of each FA was also presented for a visual check for biases.truth (green) than NNLS MWF (red, MAE=0.058). The MWF prediction errors of both methods are plotted against theFA in Figure 4 , where the mean error of each FA was also presented for a visual check for biases. A small but noticeablepositive bias (0.019, overestimation of MWF) was observed for the NNLS method, whereas SAME-ECOS did not showany obvious bias (0.007). | In-vivo experimental data: MWF maps
In-vivo GRASE data were analyzed by the SAME-ECOS model and NNLS. The processing times of the whole braindata were 3 minutes for the SAME-ECOS model and 86 minutes for NNLS. Six representative slices of the GRASEfirst echo, the resulting MWF maps produced by both methods (T s < 40ms), and the voxel-wise MWF difference map(SAME-ECOS MWF − NNLS MWF), are presented in
Figure 5 . The SAME-ECOS MWF map is visually very similar tothe NNLS MWF map, but subtle differences are still visible. Quantitatively, the SAME-ECOS mean MWF of whole whitematter (0.131 ± ± Figure 6 shows the voxel spectra, themean spectra, and the mean MWF within the ROIs of the corpus callosum, corticospinal tract, forceps major and forcepsminor. The SAME-ECOS mean MWF was lower than the NNLS mean MWF for most ROIs (except for forceps major). Itis observed that for all ROIs, the SAME-ECOS mean spectra are able to resolve a short T peak (attributed to myelinwater) in addition to a more dominate middle peak (attributed to IE water), whereas the NNLS mean spectra can onlyresolve a single broad middle peak. IU ET AL . 11
F I G U R E 5
SAME-ECOS and NNLS derived in-vivo myelin water fraction maps.
Six representative slices of theGRASE first echo, the resulting MWF maps produced by both methods (T s < 40ms), and the voxel-wise MWFdifference map (SAME-ECOS MWF − NNLS MWF), are shown.
IU ET AL ..
IU ET AL .. F I G U R E 6
Voxel spectra (gray), the mean spectra (blue: SAME-ECOS; red: NNLS), and the mean MWF (T s <40ms) within the four ROIs of the experimental in-vivo data. ROIs include the corpus callosum, corticospinal tract,forceps major, and forceps minor. The voxel spectra are only plotted for a fraction of the total number of voxels for avisual-friendly presentation as indicated in each subfigure. The mean spectrum and mean MWF are calculated from allvoxels within each ROI.
IU ET AL . 13 | DISCUSSION4.1 | SAME-ECOS vs. NNLS
From the results of all three simulation tests, the SAME-ECOS model largely outperformed NNLS. Using NNLS asthe baseline, the SAME-ECOS model achieved 13.1% higher overall cosine similarity scores (
Test 1 ), and 16.0% lowerMAE of MWF predictions (
Test 3 ), as well as demonstrated better robustness to noise (
Test 2 ). Specifically, from visualinspection of the results presented in
Figure 2 , both methods were able to produce accurate T spectra when thenumber of T components n ≤
2. However, when n > 2, the NNLS predictions became extremely unstable, resulting inover-smoothed mean spectra This observation illustrates that, unlike NNLS being highly susceptible to noise, the SAME-ECOS model has a desirable feature of being relatively noise-inert. It is also noticed that both methods were incapableof producing reliable spectrum predictions particularly for the long T components (e.g. T s > 500ms), which is likely dueto the relatively short measurement time (the last TE=320ms) so there was not enough information acquired to resolvethe long T component accurately. In addition, as illustrated in Figure 4 , NNLS was systematically overestimating MWFthroughout the FA range under our investigation. In contrast, SAME-ECOS produced more accurate MWF valueswithout any obvious biases. Overall, the current SAME-ECOS model demonstrated a better performance than NNLS inthe interrogations using simulation tests.Due to the lack of ground truth, it is difficult to judge which method is more accurate when it comes to in-vivoexperimental data. However, in terms of data processing speed, SAME-ECOS is approximately 30 times faster thanNNLS, achieving a whole-brain analysis in 3 minutes, which is more feasible in clinical settings. From a visual inspectionof the MWF maps shown in
Figure 5 , these two methods produced similar results, but the SAME-ECOS MWF mapappeared to be a little less noisy, making its white and gray matter distinction slightly more prominent. Interestingly,compared with SAME-ECOS, the NNLS mean MWF values in white matter regions were higher, which coincided withthe NNLS MWF overestimation issue presented in
Figure 4 .A thorough analysis of the NNLS spectra in a few white matter ROIs (
Figure 6 ) revealed that the NNLS MWF signal(T cutoff at 40ms) largely originated from the broad middle peak, which is commonly designated as the IE water pool.There was not much contribution from a separate myelin water pool (T s < 40ms) to the NNLS MWF, and the myelinwater peak did not exist on the NNLS mean spectra for all ROIs. In contrast, the myelin water peak and its separationwith the IE water peak were easily seen on the SAME-ECOS mean spectra for all ROIs ( Figure 6 ), aligning with theintuitive interpretation of MWF arising entirely from a separate myelin water pool. Thus, this observation raises aconcern that reporting MWF values alone may not fully unveil the underlying information, especially for those myelinwater imaging studies that used NNLS. Nevertheless, future studies using histological validation [38, 39] should beconducted to better compare these two analysis methods in dealing with experimental data. | Advantages of SAME-ECOS
The biggest advantage of SAME-ECOS lies in the fact that it is a simulation-derived solver, which takes a fundamentallydifferent route to yield the solutions compared with conventional solvers such as NNLS. Its unique workflow makesSAME-ECOS highly tailorable to different experimental conditions and the needs of analysis. Most parameters used inSAME-ECOS are tunable. For example, the T parameter was set to be a constant for all T components in our simulationbecause the T weighting is minimized by using GRASE. But if a sequence sensitive to both T and T were used, onecould simply turn the T parameter into another variable for the production of training data and make the training labelsa 2D map that resolves distinct components with different T , T times. [40] In contrast, it is presumably more difficult IU ET AL ..
The biggest advantage of SAME-ECOS lies in the fact that it is a simulation-derived solver, which takes a fundamentallydifferent route to yield the solutions compared with conventional solvers such as NNLS. Its unique workflow makesSAME-ECOS highly tailorable to different experimental conditions and the needs of analysis. Most parameters used inSAME-ECOS are tunable. For example, the T parameter was set to be a constant for all T components in our simulationbecause the T weighting is minimized by using GRASE. But if a sequence sensitive to both T and T were used, onecould simply turn the T parameter into another variable for the production of training data and make the training labelsa 2D map that resolves distinct components with different T , T times. [40] In contrast, it is presumably more difficult IU ET AL .. to include either T or any other quantities as an additional variable in the NNLS analysis. Another example would behow easily the influence of B1 inhomogeneity is handled. SAME-ECOS simply treats the refocusing FA as a variablewhen producing the simulated decay data to account for the FA variations caused by B1 inhomogeneity. However, thesame problem was not solved using NNLS for many years until Prasloski et al. [24] proposed to integrate the stimulatedecho correction into the original NNLS algorithm, which was a tremendous effort for such integration. Furthermore, ifone wants to subject the data analysis to a simpler 3-pool model for instance, SAME-ECOS can easily be simplified to a3-pool solver by fixing the number of T components to be n=3 in the simulation, although it is not encouraged to do sofor reasons discussed later. Conclusively, SAME-ECOS is extremely tunable to accommodate either a simpler or morecomplex model.Another advantage of SAME-ECOS is to utilize the concept of resolution limit from information theory to preventoverfitting of noise. Due to unavoidable noise contamination, T components can only be recovered to a certain extent,and any violation to the law of resolution limit should be theoretically prohibited. The essential part of SAME-ECOSworkflow is to obey the resolution limit to prevent any solutions with adjacent T components being too close together,a feature which is not guaranteed in the solutions of NNLS or any other methods to our knowledge. One reasonableconcern that is raised here is that nature does not know about resolution limits, so what happens when adjacent T components really are too close together? Unfortunately, the answer is that we might never resolve these two peaksin a solution where they are bound to degenerate into one peak due to the limited SNR, as Istratov et al. [3] pointedout about this particular problem: “Many physicists have discovered after much wasted effort that it is essential tounderstand the ill-conditioned nature of the problem before attempting to compute solutions”. On the other hand, thecommonly used NNLS has incorporated regularization techniques [41] to mitigate the noise overfitting, but the choiceof regularization parameters is poorly justified in most of the MR literature. From a technical perspective, the strengthof regularization should be selected based on the local SNR rather than being universally fixed. Unfortunately, thishas never been practiced in the NNLS analysis due to technical challenges. In addition, the noise of a magnitude MRIimage originates from Gaussian distributed noise on both real and imaginary channels [25], which is not accounted forin the NNLS method [42] but is correctly modelled in the SAME-ECOS approach. Particularly in the last few echoes,when there is no residuals of MR signal (only pure noise), the noise profile would follow a Rayleigh distribution on themagnitude image that always has positive values, such that NNLS could misinterpret the Rayleigh noise as a long T component that has not been completely decayed yet. This phenomenon is observed in Figure 6 , where NNLS meanspectra (red) for all ROIs always give a little rise on the far side of long T s (T = 2000ms), but SAME-ECOS mean spectra(blue) is free of this problem.SAME-ECOS also takes advantage of the strong predictive power of a fine-tuned neural network. For regressionproblems such as our T -fitting task, modern machine learning methods like neural networks usually outperformconventional statistical methods in terms of better prediction accuracy and faster data processing speed [17], which isexactly what we have observed in our simulation tests. Out of various machine learning methods, the neural networkapproach is favored for use in the current SAME-ECOS workflow because it has been successfully implemented insimilar tasks by different research groups. [15, 16, 43] In theory, other machine learning methods may also achieve asimilar predictive power. Comparisons between different machine learning methods are beyond the scope of this paper,but could be an area of future investigation.Finally, it is worth highlighting that the simulated training dataset of SAME-ECOS is not informed by any priorknowledge (e.g. a 2-pool model is informed by the T times of myelin water and IE water, or a typical range of theMWF values obtained from previous studies). Instead, SAME-ECOS is completely driven by a large number of randomsimulations, which are only regulated by the experimental conditions such as SNR. This approach should be valued andfavored because it is absolutely immune to (1) the potential errors of previous findings where the prior information IU ET AL . 15 is acquired from; (2) the biases that prior information may introduce into the analysis results. Simply speaking, aninformed model is more likely to perform the analysis with a bias naturally towards the prior information. Similar toNNLS, no prior information being needed is a desirable and specially designed feature of SAME-ECOS. | Disadvantages of SAME-ECOS
A major concern of SAME-ECOS is that its fitting results are not easily verified mathematically due to the use of aneural network. SAME-ECOS is capable of producing reliable results only on data that are similarly distributed as thetraining data. SAME-ECOS may become unpredictable and yield uninterpretable results when applied to unfamiliardata. Although we have generated a large training dataset by the simulations to account for many types of variations,there is no doubt that real experimental data have far more complexity in them. Factors such as artifacts are detrimentalto any analysis methods, but SAME-ECOS is potentially less predictable. Unlike SAME-ECOS which empirically analyzesthe data, conventional methods process the data either analytically or numerically, making them somewhat moremathematically explainable and predictable when encountering these problems.Another limitation of SAME-ECOS is related to standardization. The variety of parameters such as the T range, thenoise profiles, the number of neural network hidden layers etc. can be manually selected at the user’s will, which offersflexibility and customizability, but it is difficult to propose a standard SAME-ECOS model that works universally. Theeffects of changing these parameters warrants further investigation. Although SAME-ECOS demonstrated excellentperformance in our simulation tests and in-vivo example application, users should further validate SAME-ECOS on theirown experimental data before replacing any conventional analysis methods. | Other quantitative MRI techniques and beyond
Within MRI, the usage of SAME-ECOS is not just limited to multi-echo relaxation sequences. As long as the spinsare trackable by simulations, then the SAME-ECOS methodology should also be applicable to other quantitative MRItechniques, such as multi-component driven equilibrium single pulse observation of T and T (mcDESPOT) [44] andneurite orientation dispersion and density imaging (NODDI) [45], by modifying the current workflow accordingly.Beyond MRI, it is also possible to apply SAME-ECOS to any quantitative sciences that involve multi-exponential decays.In general, we believe a simulation-derived solver like SAME-ECOS is an alternative way to produce at least comparableresults to conventional methods. It may deliver better performance, especially when analytical and numerical solutionsstart to fail due to factors such as a limited amount of data points and low SNR. Nevertheless, the SAME-ECOSmethodology seems to have the potential to be generalized for the analysis of decay data within and beyond the MRIfield. | CONCLUSION
We have introduced a novel method SAME-ECOS, which can decompose multi-exponential MR relaxation data into a T spectrum. SAME-ECOS is highly tailorable to different experimental conditions and various analysis models. Comparedwith the commonly used method NNLS, our results have demonstrated that SAME-ECOS can yield much more reliableT spectra and MWF values in a dramatically shorter processing time, by utilizing information theory and machinelearning simultaneously. IU ET AL ..
We have introduced a novel method SAME-ECOS, which can decompose multi-exponential MR relaxation data into a T spectrum. SAME-ECOS is highly tailorable to different experimental conditions and various analysis models. Comparedwith the commonly used method NNLS, our results have demonstrated that SAME-ECOS can yield much more reliableT spectra and MWF values in a dramatically shorter processing time, by utilizing information theory and machinelearning simultaneously. IU ET AL .. A C K N O W L E D G E M E N T S
We thank the study participants and the excellent MRI technologists at UBC MRI Research Center. Funding support wasprovided by the Multiple Sclerosis Society of Canada, Natural Sciences and Engineering Research Council DiscoveryGrant.
R E F E R E N C E S [1] Whittall KP, MacKay AL, Li DK. Are mono-exponential fits to a few echoes sufficient to determine T2 relaxation forin vivo human brain? Magnetic Resonance in Medicine: An Official Journal of the International Society for MagneticResonance in Medicine 1999;41(6):1255–1257.[2] Mackay A, Whittall K, Adler J, Li D, Paty D, Graeb D. In vivo visualization of myelin water in brain by magnetic resonance.Magnetic resonance in medicine 1994;31(6):673–677.[3] Istratov AA, Vyvenko OF. Exponential analysis in physical phenomena. Review of Scientific Instruments1999;70(2):1233–1257.[4] Whittall KP, MacKay AL. Quantitative interpretation of NMR relaxation data. Journal of Magnetic Resonance (1969)1989;84(1):134–152.[5] Lawson CL, Hanson RJ. Solving least squares problems. SIAM; 1995.[6] Provencher SW. CONTIN: a general purpose constrained regularization program for inverting noisy linear algebraic andintegral equations. Computer Physics Communications 1982;27(3):229–242.[7] Kroeker RM, Henkelman RM. Analysis of biological NMR relaxation data with continuous distributions of relaxationtimes. Journal of Magnetic Resonance (1969) 1986;69(2):218–235.[8] Du YP, Chu R, Hwang D, Brown MS, Kleinschmidt-DeMasters BK, Singel D, et al. Fast multislice mapping of the myelinwater fraction using multicompartment analysis of T decay at 3T: A preliminary postmortem study. Magnetic Resonancein Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine 2007;58(5):865–870.[9] Akhondi-Asl A, Afacan O, Mulkern RV, Warfield SK. T2-Relaxometry for Myelin Water Fraction Extraction Using WaldDistribution and Extended Phase Graph. In: International Conference on Medical Image Computing and Computer-Assisted Intervention Springer; 2014. p. 145–152.[10] Bertero M, Boccacci P, Pike ER. On the recovery and resolution of exponential relaxation rates from experimental data:a singular-value analysis of the Laplace transform inversion in the presence of noise. Proceedings of the Royal Society ofLondon A Mathematical and Physical Sciences 1982;383(1784):15–29.[11] McWhirter J, Pike ER. On the numerical inversion of the Laplace transform and similar Fredholm integral equations ofthe first kind. Journal of Physics A: Mathematical and General 1978;11(9):1729.[12] Link N, Bauer S, Ploss B. Analysis of signals from superposed relaxation processes. Journal of applied physics1991;69(5):2759–2767.[13] Schmidhuber J. Deep learning in neural networks: An overview. Neural networks 2015;61:85–117.[14] Lundervold AS, Lundervold A. An overview of deep learning in medical imaging focusing on MRI. Zeitschrift für Medi-zinische Physik 2019;29(2):102–127.[15] Cohen O, Zhu B, Rosen MS. MR fingerprinting deep reconstruction network (DRONE). Magnetic resonance in medicine2018;80(3):885–894.
IU ET AL . 17 [16] Liu H, Xiang QS, Tam R, Dvorak AV, MacKay AL, Kolind SH, et al. Myelin water imaging data analysis in less than oneminute. NeuroImage 2020;210:116551.[17] Hornik K, Stinchcombe M, White H, et al. Multilayer feedforward networks are universal approximators. Neural net-works 1989;2(5):359–366.[18] Prasloski T, Rauscher A, MacKay AL, Hodgson M, Vavasour IM, Laule C, et al. Rapid whole cerebrum myelin water imagingusing a 3D GRASE sequence. Neuroimage 2012;63(1):533–539.[19] Deshmane A, Gulani V, Griswold MA, Seiberlich N. Parallel MR imaging. Journal of Magnetic Resonance Imaging2012;36(1):55–72.[20] Shapiro FR, Senturia SD, Adler D. The use of linear predictive modeling for the analysis of transients from experimentson semiconductor defects. Journal of applied physics 1984;55(10):3453–3459.[21] Dobaczewski L, Kaczor P, Hawkins I, Peaker A. Laplace transform deep-level transient spectroscopic studies of defectsin semiconductors. Journal of applied physics 1994;76(1):194–198.[22] Thomasson W, Clark Jr J. Analysis of exponential decay curves: A three-step scheme for computing exponents. Mathe-matical Biosciences 1974;22:179–195.[23] Hennig J. Multiecho imaging sequences with low refocusing flip angles. Journal of Magnetic Resonance (1969)1988;78(3):397–407.[24] Prasloski T, Mädler B, Xiang QS, MacKay A, Jones C. Applications of stimulated echo correction to multicomponent T2analysis. Magnetic resonance in medicine 2012;67(6):1803–1814.[25] Cárdenas-Blanco A, Tejos C, Irarrazaval P, Cameron I. Noise in magnitude magnetic resonance images. Concepts inMagnetic Resonance Part A: An Educational Journal 2008;32(6):409–416.[26] Klambauer G, Unterthiner T, Mayr A, Hochreiter S. Self-normalizing neural networks. In: Advances in neural informationprocessing systems; 2017. p. 971–980.[27] Bishop CM. Pattern recognition and machine learning. springer; 2006.[28] Kingma DP, Ba J. Adam: A method for stochastic optimization. arXiv preprint arXiv:14126980 2014;.[29] LeCun Y, Bengio Y, Hinton G. Deep learning. nature 2015;521(7553):436–444.[30] Abadi M, Agarwal A, Barham P, Brevdo E, Chen Z, Citro C, et al. Tensorflow: Large-scale machine learning on heteroge-neous distributed systems. arXiv preprint arXiv:160304467 2016;.[31] Dvorak AV, Ljungberg E, Vavasour IM, Liu H, Johnson P, Rauscher A, et al. Rapid myelin water imaging for the assessmentof cervical spinal cord myelin damage. NeuroImage: Clinical 2019;23:101896.[32] Liu H, Rubino C, Dvorak AV, Jarrett M, Ljungberg E, Vavasour IM, et al. Myelin water atlas: a template for myelin distri-bution in the brain. Journal of Neuroimaging 2019;29(6):699–706.[33] Liu H, Ljungberg E, Dvorak AV, Lee LE, Yik JT, MacMillan EL, et al. Myelin water fraction and intra/extracellular watergeometric mean T2 normative atlases for the cervical spinal cord from 3T MRI. Journal of Neuroimaging 2020;30(1):50–57.[34] Lee LE, Ljungberg E, Shin D, Figley CR, Vavasour IM, Rauscher A, et al. Inter-vendor reproducibility of myelin waterimaging using a 3D gradient and spin echo sequence. Frontiers in Neuroscience 2018;12:854.[35] Ljungberg E, Vavasour I, Tam R, Yoo Y, Rauscher A, Li DK, et al. Rapid myelin water imaging in human cervical spinal cord.Magnetic resonance in medicine 2017;78(4):1482–1487.
IU ET AL ..