- Email: [email protected]
Optimization of q-space sampling for mean apparent propagator MRI metrics using a genetic algorithm
NeuroImage 199 (2019) 237–244
Contents lists available at ScienceDirect
NeuroImage journal homepage: www.elsevier.com/locate/neuroimage
Optimization of q-space sampling for mean apparent propagator MRI metrics using a genetic algorithm Daniel V. Olson a, e, *, Volkan E. Arpinar b, d, L. Tugan Muftuler c, d a
Department of Biophysics, Medical College of Wisconsin, Milwaukee, WI, USA Department of Radiology, Medical College of Wisconsin, Milwaukee, WI, USA Department of Neurosurgery, Medical College of Wisconsin, Milwaukee, WI, USA d Center for Imaging Research, Medical College of Wisconsin, Milwaukee, WI, USA e Magnetic Resonance Imaging, GE Healthcare, Waukesha, WI, USA b c
A R T I C L E I N F O
A B S T R A C T
Keywords: Mean apparent propagator (MAP) Optimization q-space sampling Genetic algorithm
Mean Apparent Propagator (MAP) MRI is a recently introduced technique to estimate the diffusion probability density function (PDF) robustly. Using the estimated PDF, MAP MRI then calculates zero-displacement and nonGaussianity metrics, which might better characterize tissue microstructure compared to diffusion tensor imaging or diffusion kurtosis imaging. However, intensive q-space sampling required for MAP MRI limits its widespread adoption. A reduced q-space sampling scheme that maintains the accuracy of the derived metrics would make it more practical. A heuristic approach for acquiring MAP MRI with fewer q-space samples has been introduced earlier with scan duration of less than 10 minutes. However, the sampling scheme was not optimized systematically to preserve the accuracy of the model metrics. In this work, a genetic algorithm is implemented to determine optimal q-space subsampling schemes for MAP MRI that will keep total scan time under 10 min. Results show that the metrics derived from the optimized schemes more closely match those computed from the full set, especially in dense fiber tracts such as the corpus callosum.
1. Introduction Diffusion MRI (dMRI) is a powerful technique capable of probing tissue microstructure in vivo at a resolution that is much higher than the imaging resolution. Diffusion MRI achieves this by indirectly estimating the net displacement of water molecules in the biological environment. Effective diffusion modeling necessitates a series of diffusion-weighted images (DWI). However, acquiring a DWI set with sufficient diffusionencoding vectors (b-vectors) for advanced modeling is time-intensive. This makes advanced dMRI techniques impractical when working with certain patient populations and unfeasible for clinical use. The distribution of acquired b-vectors affects robustness to sample orientation and noise. At the same time, model fitting can be overdetermined such that a subset of DWI might yield similar results as the full set (Dubois et al., 2006). An ideal q-space sampling scheme balances robustness with minimizing redundancy in the measurements. For the most part, q-space sampling schemes have been designed to approximate uniform angular coverage on the unit sphere. Jones et al. formulated the optimization as a numerical minimization of electrostatic
charge distributions on a sphere (Jones et al., 1999). For 30 directions, this electrostatic solution approximates the icosahedral group which is rotationally invariant (Batchelor et al., 2003; Hasan and Narayana, 2003). An analytical model for distributing antipodally symmetric points on a sphere has been developed, which closely matches the electrostatic model (Koay, 2011). Extending the optimization to multiple shells is required for advanced diffusion models, as single-shell acquisitions are innately limited in the radial direction in q-space. Standard multi-shell acquisitions have uniformly distributed sample points on each shell. Recent works have also sought to maximize the angular incoherence between q-space samples taken across shells (Koay et al., 2012; Caruyer and Deriche, 2012; Caruyer et al., 2013). Generally, more samples are acquired on the high b-value shells (Wu and Alexander, 2007; Assaf and Basser, 2005; € Ozarslan et al., 2013). This results in more uniform sampling but also denser sampling where the diffusion signals are highly attenuated, which has raised concerns that the rectified noise floor will influence modeling (Andre et al., 2014; Hutchinson et al., 2017). Optimizing q-space sampling for a specific diffusion model can
* Corresponding author. Department of Biophysics, Medical College of Wisconsin, Milwaukee, WI, USA. E-mail address: [email protected] (D.V. Olson). https://doi.org/10.1016/j.neuroimage.2019.05.078 Received 4 March 2019; Received in revised form 28 April 2019; Accepted 29 May 2019 Available online 1 June 2019 1053-8119/© 2019 Elsevier Inc. All rights reserved.
D.V. Olson et al.
NeuroImage 199 (2019) 237–244
produce more accurate parameter estimates than general optimization. An example is the optimization for Diffusion Kurtosis Imaging (DKI) performed by Poot et al. (2010), which used simulated annealing as its search method. The optimal q-space sampling did not conform to a multi-shell distribution. Instead, the b-values had a varied radial distribution with a substantial portion of the gradients at the maximum allowed b-value of 2800 s/mm2. The b-vectors were approximately isotropically distributed on the unit sphere. A separate multi-shell scheme from that work performed best in an independent study (Sprenger et al., 2016). This suggests that a departure from multi-shell acquisitions to q-space sampling schemes of higher radial resolution could better characterize the diffusion probability density function. A recently proposed diffusion model called Mean Apparent Propagator (MAP) MRI is a mathematical framework connecting the diffusion€ weighted imaging q-space with the molecular displacements (Ozarslan et al., 2013). The basis functions are eigenfunctions of the quantum-mechanical analog of the simple harmonic oscillator Hamiltonian. The MAP MRI model yields novel metrics to characterize the diffusion properties. Derivation of the metrics from the harmonic co€ efficients is demonstrated in Ozarslan et al. (2013). While MAP MRI may offer better characterization of diffusion in complex tissue structures compared to DTI or DKI, an intensive acquisition impedes its widespread adoption, both clinically and in research. € The original MAP acquisition scheme proposed in Ozarslan et al. (2013) consisted of seven b-values of b ¼ {200, 800, 1800, 3200, 5000, 7200, 9800} s/mm2 along 5, 14, 32, 56, 87, 125, and 170 directions, respectively. An excised marmoset brain was scanned, so there were no concerns about long acquisition time or motion. A human adaptation proposed by Avram et al. (2016) reduced the acquisition time to 72 min. The sampling scheme consisted of six b-values b ¼ {1000, 2000, 3000, 4000, 5000, 6000} s/mm2 along 23, 43, 83, 131, 148, and 256 directions, respectively. A heuristically-derived q-space subsampling scheme was also proposed in that paper for clinical feasibility, which further reduced the number of directions per shell to 4, 7, 11, 17, 23, and 31 directions for an acquisition of about 10 min. In this work, we describe an optimization framework implemented for MAP MRI. A genetic algorithm is employed, which determined clinically feasible subsets of q-space sampling schemes that minimize the mean-squared error (MSE) relative to the MAP metrics produced by the full q-space. The optimized q-space sampling schemes were validated using both a digital phantom and also in vivo data. In the next section, genetic algorithms are reviewed first to provide a context for our implementation.
Fig. 1. Overview of genetic algorithm.
should refer to the citations provided. The GA consisted of four parts: initialization, fitness evaluation, selection, and breeding. Initialization is the construction of the first population. For the initial population, individuals were assigned random subsamples of b-vectors from the full set without prohibition of duplicates. Future generations are determined through iterations of evaluation, selection, and breeding. Each individual was then evaluated with MAP fitting, and the MSE was calculated relative to the true value of the metric from the full data set, which was similar to € the 489 directions scheme used in Ozarslan et al. (2013). The selection method was Stochastic Universal Sampling (SUS), a fitness-proportionate method that minimizes deviation of the selections from their reproduction expectation values (Baker, 1987; Mitchell, 1998). Fitness values were adaptively scaled with sigma scaling (Forrest and Mitchell, 1993). In populations with high fitness variance, sigma scaling limits the impact of individuals with outlying fitness values, as would occur in early generations. Conversely, low fitness variance amplifies the differences between individuals in later generations to maintain a relatively constant selection pressure over the duration of the algorithm (Back, 1994). A final selection concept is elitism, where a small number of the best individuals are guaranteed survival (De Jong, 1975). GA with elitism have been shown to outperform those without it in almost every trial (Baluja and Caruana, 1995). Breeding is comprised of the genetic operators of crossover and mutation. Ninety-eight percent of the offspring were generated through crossover. In uniform crossover, each b-vector is inherited independently from any other for maximal disruption without loss of b-vectors (Syswerda, 1989). A disruptive crossover operator is desired with small population sizes to more widely sample the solution space and avoid the tendency towards population homogeneity (Spears and De Jong, 1990; Eshelman, 1990). Mutations promote a greater traversal of solution space by a random walk. Without mutation, only b-vectors that were present in the initial population would exist in the simulations. Mutation was defined as exchanging a single b-vector for another in the full set. with a probability of 0.01. Mutation rate should be kept low to allow for the crossover operator to determine better q-space sampling schemes. Evaluation, selection, and breeding are repeated until stopping criteria are satisfied. The maximum number of generations was set to 250 to allow for sufficient evolutions while preventing overfitting of the signals. Additional stopping criteria were a stall time of 75 generations as well as convergence tolerances of 106 for non-Gaussianity (NG) and 0.01 for return-to-origin probability (RTOP). In this simulation, convergence was defined as the difference of the lowest individual MSE and the average MSE of the population being below the tolerance. These values correspond to approximately 0.1% error for the two metrics.
2. Methods 2.1. Genetic algorithms A genetic algorithm (GA) is a metaheuristic global optimization solver that utilizes evolutionary biology principles to reach a desired solution (Holland, 1973). GA are population-based solvers that make few assumptions about the problem and thus are considered “weak methods” (Whitley, 1994). Beneficial features are selected and recombined, approaching an optimal tuning over iterations. Even in noisy environments, GA are capable of approaching a global minimum, unlike most other optimizers (Fitzpatrick and Grefenstette, 1988; Beyer, 2000; Rakshit et al., 2017). Although a global optimum is not guaranteed, heuristic solvers can find an adequate solution in a reasonable amount of time. Consider the case of selecting an optimal subset of 95 diffusion-encoding gradient vectors from a collection of 489, which was used in the original € paper by Ozarslan et al. (2013). The binomial coefficient is on the order of 10103 possible combinations. Further, acceptable sampling schemes are sparse in solution space. Therefore, exhaustive and random-walk searches are not feasible. An overview of the GA is presented in Fig. 1, and a glossary of GA terminology is provided in Appendix B. For further details, the readers 238
D.V. Olson et al.
NeuroImage 199 (2019) 237–244
calculated from the MAP coefficients based on the respective q-space sampling schemes for the full set and each subsampled set. The GARTOP and GANG sampling schemes were taken directly from the output of the optimization without modification. For comparison, the reduced sampling scheme proposed by Avram et al. (2016) was also included in the analyses (referred to as Avram scheme). The diffusion-encoding directions for this set were selected through a Monte Carlo farthest neighbor search to approximate a uniformly distributed subset on each shell, but b-vectors were not optimized across shells to maximize angular incoherence on the unit sphere. A subsequent MAP fitting was applied to these simulated signals to calculate the zero-displacement and non-Gaussianity metrics. MSE were calculated for each subsampled set's metrics relative to the corresponding metrics from the full set.
2.2. Acquisition and processing of MAP MRI data for optimization In vivo data were acquired from a healthy volunteer (male, 26 years). The study was approved by the Institutional Review Board, and written informed consent was obtained from the participant. Data were collected on a General Electric (GE) SIGNA Premier 3 T MRI scanner with a 48channel receive head coil. Diffusion weighted images (DWI) were acquired using a single-shot SE-EPI sequence. The imaging volume was 20 slices with an acquisition matrix of 100 100. Voxel dimensions were 2.2 2.2 2 mm. The echo time (TE) was 86.8 ms, and the repetition time (TR) was 4200 ms. Bandwidth was 1953 Hz/pixel. Parallel imaging (ASSET) in-plane acceleration factor of 2 was used to reduce susceptibility artifacts, but no multiband acceleration was used. Acquisition time was 35 min. The data set consisted of 10 interspersed, non-diffusionweighted reference images and 489 DWI. The DWI were sampled in six shells with b ¼ {1000, 2000, 3000, 4000, 5000, 6000} s/mm2 along 19, 32, 56, 87, 125, and 170 directions, respectively. In each shell, diffusionencoding directions were uniformly distributed over the unit sphere (Jones et al., 1999). However, the angular incoherence of the b-vectors was not maximized across the shells. The maximum diffusion gradient amplitude was 69 mT/m with pulse width δ of 22.8 ms and a separation Δ of 41.3 ms. This data set will be called the full data set throughout the rest of this paper. Postprocessing of the data set included EPI phase correction on the scanner. Offline processing included motion and eddy current corrections using the EDDY function of the FSL software (Jenkinson et al., 2012; Andersson and Sotiropoulos, 2016). A brain mask was calculated with BET (Smith, 2002). Then, MAP fitting was performed using the DIPY Toolbox (Garyfallidis et al., 2014). The MAP was fit with radial order of 6. No denoising was applied to the data prior to MAP fitting, but a positivity constraint and Laplacian regularization with a weighting of 0.05 were applied during fitting (Fick et al., 2015). Those settings were recommended for robust fitting by the author of the DIPY MAPMRI toolbox. An initial FA map was also generated from the data set to identify white matter voxels. This FA map was thresholded at FA>0.3 to generate a white matter mask.
2.5. Validation of genetic algorithm results using emulation by down sampling To quantify the accuracy of the MAP metrics for each subsampling scheme in vivo, a similar analysis was performed on the acquired DWI. To emulate separate acquisitions, the full data set was downsampled to incorporate each of the subsampling schemes. The subsampled schemes were not acquired separately to minimize the variability between the sets introduced by subject motion and other physiological noise sources. From the white matter (FA > 0.3), voxels that were not included in the GA optimization were used to test the accuracy of the reduced q-space sampling schemes relative to the full set. MAP model fitting and MSE evaluations were performed in the same manner as with the numerical simulation. The distribution of RTOP and NG metrics are presented in Fig. 2. The distribution closely matched the values used in the GA optimization. 2.6. Second in vivo experiment for independent validation To test the generalizability of the proposed q-space sampling schemes, the same acquisitions were also performed on another healthy volunteer (female, 26 years), whose data was not used in optimization. Written informed consent was obtained from the participant before the study. Data were collected using identical scanning parameters as the first acquisition. The full MAP MRI data set was acquired as before, and the same postprocessing and MAP fitting pipelines were applied, except no FA masking was used.
2.3. Genetic algorithm optimization The genetic algorithm was developed in Python 2.7.14 to interface with the MAP functions of the DIPY Toolbox. Two optimizations were performed separately for return-to-origin probability and nonGaussianity. The resulting q-space sampling schemes are denoted GARTOP and GANG. Diffusion signals were measured from 400 voxels in the white matter. Poot et al. (2010) determined that white matter could be adequately represented with 400 voxels. These voxels were randomly selected from within the FA mask. In determining population size, there is a trade-off between sufficient population diversity to sample solution space and computational resources for more generations. In this optimization, the MAP model fitting is computationally intensive. Diversity has been shown to become independent of the individual size with population sizes of about 200 (Diaz-Gomez and Hougen, 2007). There were 200 individuals in each population for the GA, which consisted of 95 b-vectors from the full set. This corresponds to approximately a 10-min acquisition. The GA evaluation function was the MSE relative to the full set's metrics of either RTOP or NG. Selection was fitness-proportionate with Stochastic Uniform Sampling, as well as sigma scaling to maintain selection pressure. Uniform crossover was used to increase mixing. Elite survival rate, crossover rate, and mutation rate were 0.02, 0.98, and 0.01, respectively.
3. Results Signal-to-Noise Ratios (SNR) were estimated on diffusion data that was used in the optimization. Since the signals in DWI are dependent on the diffusion-encoding direction, SNR can vary substantially depending on the relative direction of the diffusion-encoding and tissue structure. In order to determine the worst-case SNR, the signal should be chosen in a region whose fiber direction corresponds with the diffusion-encoding direction
2.4. Numerical simulation with a digital brain MAP MRI phantom
Fig. 2. Distribution of RTOP1/3 and NG values for simulations. Training set (blue) consisted of 400 random voxels from the white matter (FA > 0.3) and was used in the genetic algorithm optimization. Test set (red) consisted of the remaining voxels from the white matter. The test set was used in the quantitative validation of q-space sampling schemes.
The reduced sampling schemes were first tested using a digital brain phantom. To generate the digital phantom, an initial MAP fitting was performed on the full set. Then, diffusion-weighted signals were 239
D.V. Olson et al.
NeuroImage 199 (2019) 237–244
(Descoteaux et al., 2011; Jones et al., 2012). This will show the most attenuated signal for a given b-value and more clearly indicate image quality. To calculate SNR values, regions of interest were manually drawn in the center of the splenium of the corpus callosum (which corresponded with the Gx direction) and the background. SNR levels were 111, 38.6, 18.0, 12.3, 8.6, 8.2, and 8.0 for respective b-values of 0–6000 s/mm2. DWI at b ¼ 6000 s/mm2 retained sufficient signal such that tissue anatomy was still clearly delineated. Interspersed non-diffusion-weighted images were acquired every 50 diffusion-weighted images. Motion artifacts were not observed in the DWI. The average temporal signal-to-noise ratio estimated from the non-diffusion-weighted images was 22.2 (Fig. 3). 3.1. Genetic algorithm and the optimized q-space sampling schemes The GA schemes acquired fewer directions on the middle shells than the innermost and outmost shells (Fig. 4). The final result shows that the number of samples from the b ¼ 1000 s/mm2 shell increased for both GARTOP and GANG schemes compared to the Avram scheme. This finding suggests that acquiring more samples on shells with lower b-values may be needed for robust fitting. On the other hand, collecting more directions on higher b-value shells has been generally recommended for more accurate MAP fitting. Still, no evidence of this requirement has been demonstrated to date. Our GA-based optimization clearly shows that more samples at high b-values are essential. In the original paper € that introduced the MAP MRI technique, Ozarslan et al. (2013) did ex vivo scans and went up to b ¼ 9800 s/mm2. In vivo human scans for MAP MRI typically goes up to b ¼ 6000 s/mm2 for practical reasons. Therefore, we also limited the outer shell to b ¼ 6000 s/mm2, which might be adequate if the radial order of the MAP MRI series expansion is kept at 6. Another consideration is the angular q-space coverage. Each shell of the full set contained isotropically distributed b-vectors determined by the electrostatic charge model, but diffusion gradient directions were not optimized across shells. The Avram scheme was determined by a Monte Carlo farthest neighbor search that found a subset of b-vectors within each shell with the highest angular difference. The distributions of bvectors on the unit sphere for each of the optimized sampling schemes tested here were approximately isotropic.
Fig. 4. Distribution of b-values for sampling schemes. Avram scheme (blue) increases the number of b-vectors per shell with b-value. The samplings optimized for RTOP1/3 (green) and NG (yellow) acquire more directions on the b ¼ 1000 s/mm2 shell, fewer directions in the middle shells, and a high number on the outer b ¼ 6000 s/mm2 shell.
were computed (Fig. 5). The ground truth for the MSE calculations were the metrics calculated from the respective full q-space sampling set. GA optimization produced schemes with lower errors for all metrics, compared to the heuristic subsampling scheme. GARTOP outperformed the other schemes for the zero-displacement metrics, as well as slightly outperforming GANG in the non-Gaussianity metrics. Both GA optimization schemes performed exceptionally well, particularly for nonGaussianity measures. 3.3. In vivo experiments for metric validation Validation was performed using diffusion data from the two subjects. For the first subject, only white matter voxels that were not used in optimization were used in the analyses to reduce bias. For independent verification of optimized sampling schemes, all white matter voxels from the second subject were used. Since those voxels were not used for optimization, a bias would not be expected. Both GA schemes reduced the MSE for all MAP metrics compared to the Avram scheme in white matter (Fig. 6). With in vivo data, GARTOP resulted in the lowest MSE values for the zero-displacement metrics, while GANG had the best performance for every non-Gaussianity metric. This white matter analysis
3.2. Numerical simulation using a digital MAP MRI brain phantom For each of the reduced sampling schemes, MSE of the MAP metrics
Fig. 5. Mean-squared error (MSE) of MAP metrics in white matter from numerical simulation. Avram q-space sampling (blue), genetic algorithm optimized for RTOP1/3 (green), genetic algorithm optimized for NG (yellow). MSE are relative to metrics calculated from the full q-space sampling set. Diffusion signals were calculated from MAP model coefficients from voxels in the white matter (FA > 0.3). GA schemes reduce MSE compared with the Avram scheme. The sampling optimized for RTOP provided the best overall results. Units of zero-displacement metrics are mm1.
Fig. 3. Temporal SNR of non-diffusion-weighted (b ¼ 0 s/mm2) volumes. Temporal SNR was calculated as the mean signal across the volumes normalized by the standard deviation. 240
D.V. Olson et al.
NeuroImage 199 (2019) 237–244
closely matched the predictions of the optimization. White matter MSE trends were similar between the two subjects, which implies generalizability of the optimized schemes. There were pronounced improvements in the non-Gaussianity metrics for both of the optimized sampling schemes. While the relative differences between the optimized sampling schemes and the Avram scheme in zerodisplacement probabilities (RTOP1/3, RTAP1/2 and RTPP) were not as high as the first set, there were still notable decreases in MSE for all three metrics. This analysis was repeated for the whole brain, since MAP values in structures such as deep brain gray matter could also be informative for various diseases. In the whole brain MSE calculations, voxels with RTOP1/3 values below 20 were excluded. These voxels likely contained CSF and were predominantly from the ventricles. The improvement in the non-Gaussianity metrics of the GA optimized schemes was less pronounced in whole brain MSE, yet it is still notable (Fig. 7). The trends of GARTOP performing best in the zero-displacement metrics and GANG in the non-Gaussianity measures are evident in whole brain MSE as well. RTAP1/2 errors were approximately halved in the whole brain analysis relative to the white matter analysis. This is likely an averaging effect from including gray matter, which has lower RTAP1/2 values than white matter. As a reminder, the optimization was run with an ensemble of voxels selected from the white matter. If a study requires more emphasis on, for instance, deep gray matter nuclei, the optimization could be run with additional voxels included from those regions. Qualitatively, each of the reduced schemes produced sufficiently accurate RTOP1/3 fits (Fig. 8, top). Slight improvement can be seen in the white matter with the GA sampling schemes. Similar to results found in Avram et al. (2016), fitting differences were observed around the edges of the ventricles, possibly due to partial volume effects. The low intensity areas observed in all RTOP1/3 maps near the gray matter could be CSF. The effects were most pronounced in the map from the full set. These regions were excluded from the MSE analysis with an RTOP1/3 > 20 threshold. As evidenced in the NG difference maps (Fig. 8, bottom), the GANG scheme most closely replicated the NG of the full set. GARTOP was relatively accurate in the white matter but slightly overestimated the gray matter NG. The Avram scheme performed similar to GARTOP except in the corpus callosum. Overall, both GARTOP and GANG performed well in highly organized fibers such as the splenium and genu producing values much closer to the full set.
Fig. 7. Mean-squared error of MAP metrics in whole brain. Avram sampling (blue), genetic algorithm optimized for RTOP1/3 (green), genetic algorithm optimized for NG (yellow). Each pair of bars represent results from the two subjects. For each pair, the left bar is from the first subject and the right bar is from the second subject. MSE are relative to metrics calculated from the full qspace sampling set. The sampling optimized for NG provided the best overall results. However, the GARTOP scheme did perform better for the zerodisplacement probabilities. Units of zero-displacement metrics are mm1.
The observed NG was slightly lower in the splenium and genu relative to other white matter tracts. In those fiber tracts, we expected to see high NG values as have usually been observed in DKI. Our initial interpretation was that Gibbs ringing might be affecting the NG estimates. Upon further inspection, however, we observed that the lower NG values
Fig. 6. Mean-squared error (MSE) of MAP metrics in white matter. Avram qspace sampling (blue), genetic algorithm optimized for RTOP1/3 (green), genetic algorithm optimized for NG (yellow). Each pair of bars represent results from the two subjects. For each pair, the left bar is from the first subject and the right bar is from the second subject. MSE are relative to metrics calculated from the full q-space sampling set. Diffusion signals were from voxels in the white matter (FA > 0.3). GA schemes reduce MSE compared with the Avram scheme. The sampling optimized for NG provided the best overall results. Units of zerodisplacement metrics are mm1.
Fig. 8. RTOP1/3 and NG maps from in vivo validation. RTOP1/3 metric maps (first row) and absolute difference images (second row) for each sampling scheme. Each reduced set performed similarly. NG metric maps (third row) and absolute difference images (bottom row) for each sampling scheme. GANG most closely matched the full set, especially in the major white matter tracts. The Avram scheme's largest deviation was in the corpus callosum. 241
D.V. Olson et al.
NeuroImage 199 (2019) 237–244
results. However, the same improvements were seen in the numerical simulation and second data set, which increases confidence that such artifacts did not play a significant role in these findings. While the Monte Carlo farthest-neighbor search to replicate the Avram scheme did fairly well at subsampling uniformly within each shell, it is possible that optimizing the angular incoherence of the b-vectors across shells would have improved fitting. However, in that case, the scheme would no longer be a subset of the full set. For the GA schemes, q-space interpolation could be used to synthesize diffusion signals beyond those acquired in the full set. Interpolation of qspace has been used in previous works including Hybrid Diffusion Imaging (HYDI) (Wu and Alexander, 2007) and a propagator estimation method by Ye et al. (2012). A future direction could be applying the MAP model to interpolated q-space signals to determine an optimal q-space sampling scheme with higher radial sampling density than the six shells acquired in this work. Diffusion signals should be synthesized in advance for a set of b-vectors. Still, more degrees of freedom rapidly expand the solution space, so the GA parameters would need to be adapted. Furthermore, the accuracy of such q-space interpolation needs to be validated before attempting such an approach. Computation time almost entirely depended on the MAP fitting. This was a limiting factor on the population size as well as the number of white matter voxels included. The MAP fitting must be performed for each individual on each voxel, which could be parallelized in future implementations. Beyond the fitting, minimal calculations are required for the algorithm. Fitness values and genetic operators are calculated once per generation, and computation time has negligible dependence on population size. Since this type of optimization is run only once for a particular study, computational time is not a major concern. This study focused on optimization for healthy volunteers. In patient populations where the diffusion properties are dramatically different than healthy tissue, the optimization can be repeated with patient data.
closely followed the whole splenium into the occipital lobe, where Gibbs ringing would not be expected. A similar trend was seen in the genu as well. There is further evidence that this trend actually reflects the underlying anatomy. For instance, the two lateral branches of the genu are known to have lower dispersion than in the midline (Ronen et al., 2014; Budde and Annese, 2013; Mollink et al., 2017). The NG images reflect this anatomical feature well in the full set, GARTOP, and GANG sets, while they are generally obscured in the Avram scheme. Therefore, there is enough evidence that the NG might be accurately reflecting the actual anatomical features with the optimized MAP acquisition schemes. In order to minimize the possibility that the measurements in the splenium and genu might be affected by Gibbs ringing, we reprocessed the data from the first subject using the DIFFPREP tool of TORTOISE, which has an option to remove Gibbs ringing from the DWI data (Pierpaoli et al., 2010). The measurements in the splenium and genu in the NG and RTAP1/2 images did not change after Gibbs ringing removal, supporting the view that the values in those regions reflect the anatomical features, not image artifacts. 4. Discussion The genetic algorithm implemented in this work determined two qspace sampling schemes, which could promote more widespread use of MAP MRI. Further reductions in the number of b-vectors is possible through this method, although accuracy and robustness would need to be carefully considered. It is important to note that this GA optimization merely tries to use fewer q-space samples and replicate the MAP metrics calculated from the full set. It does not find an optimal q-space sampling scheme for MAP MRI, in general. The methods and results presented here is a proof of concept to demonstrate that an optimal q-space sampling can be determined using a systematic approach. The goal was to minimize sampling redundancies without significant loss of accuracy for MAP MRI metric calculations. While GA provided such optimal acquisition schemes, it is possible that other systematic search methods could provide equally good or better results. The goal here was not to find the best performance among various optimization schemes but to show that it is possible to perform a systematic search and find a reduced set of q-space samples that delivered MAP MRI metrics that matched those of the full acquisition. It is possible that some of the voxels used in optimization or testing were affected by various sources of artifacts, such as cardiac pulsation, physiological motion, and eddy currents. Therefore, a digital brain phantom data was generated by calculating the diffusion signals using the MAP coefficients. This way, the effects of acquisition artifacts were mitigated by making the voxel data consistent for each sampling scheme. This simulation confirmed that the GA optimized subsampling schemes produced metrics that closely matched those of the full data set. Additionally, an independent data set was acquired from a different subject in order to test the robustness of the optimized sampling schemes. Data from a different subject with different anatomy and different time-varying artifacts would show that the optimized schemes could be generalized. The optimized schemes outperformed the Avram scheme for both in vivo data sets, indicating that the optimized schemes could be generalized and were not significantly affected by such sources of artifacts. A limitation of this work is that the reduced q-space samples were subsampled from the full set for analyses. This approach was chosen to avoid possible sources of error such as misregistration between scans, motion artifacts, and other effects of post-processing. If each subsampling scheme were acquired separately, the results would be skewed unfairly if the subject moved slightly more in one of the acquisitions compared to the others. It is still possible that one of the subsampling schemes coincidentally selected more of the volumes with artifacts and biased the
5. Conclusion MAP MRI is a novel technique that generates quantitative maps which could be sensitive and specific to various pathological changes in tissues. The metric maps are derived from the diffusion propagator that is estimated from a broader sampling of q-space. As a result, more subtle changes in tissues could be detected when compared to the ellipsoid assumption of DTI or second order effects modeled by kurtosis. It also does not make assumptions about the tissues under investigation. Acquisition time is a major impediment to more widespread adoption of MAP MRI. With reduced q-space sampling, data can be acquired in a clinically feasible time. However, the accuracy of metrics could be compromised if an optimal subsampling scheme was not devised. Among various optimization approaches, genetic algorithms are effective tools for q-space sampling optimization because they efficiently search a vast solution space by mimicking evolutionary biology principles. Here, we demonstrated that the GA approach can determine which diffusion-encodings could be eliminated with minimal influence on the final MAP MRI images to reproduce the metrics of the full set. While we demonstrated this optimization approach for MAP MRI, the methods described here can be adapted to other diffusion MRI techniques that require a large number of diffusion-encoding vectors. Acknowledgements The authors would like to thank Dr. Matthew Budde for helpful suggestions in the interpretation of results and preparation of the manuscript. This work is supported by funding provided by the Daniel M. Soref Charitable Trust.
242
D.V. Olson et al.
NeuroImage 199 (2019) 237–244
Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi.org/10.1016/j.neuroimage.2019.05.078. Appendix B. Glossary of genetic algorithm terminology Term
Definition
Breeding Crossover Crossover Rate Elite Survival Rate Elitism Evaluation Fitness Generation Individual Mutation Mutation Rate Offspring Population Selection
Application of crossover and mutation to produce new offspring Recombination of b-vectors from two individuals to create an offspring Proportion of the subsequent population created through crossover Proportion of the subsequent population created through elitism Copying the best individuals into the next generation Determination of the fitness of each individual Measure of an individual's performance Iteration of the genetic algorithm q-space sampling scheme Random replacement of a b-vector Probability that a given b-vector will be replaced Individual in the population of the subsequent generation Group of individuals in the generation Process of determining the probabilities that each individual will breed
Forrest, S., Mitchell, M., 1993. What makes a problem hard for a genetic algorithm? Some anomalous results and their explanation. Mach. Learn. 13 (2–3), 285–319. https:// doi.org/10.1007/BF00993046. Garyfallidis, E., Brett, M., Amirbekian, B., Rokem, A., van der Walt, S., Descoteaux, M., Nimmo-Smith, I., 2014. DIPY contributors, DIPY, a library for the analysis of diffusion MRI data. Front. Neuroinf. 8 (8) https://doi.org/10.3389/ fninf.2014.00008. Hasan, K.M., Narayana, P.A., 2003. Computation of the fractional anisotropy and mean diffusivity maps without tensor decoding and diagonalization: theoretical analysis and validation. Magn. Reson. Med. 50 (3), 589–598. https://doi.org/10.1002/ mrm.10552. Holland, J.H., 1973. Genetic algorithms and the optimal allocation of trials. SIAM J. Comput. 2 (2), 88–105. https://doi.org/10.1137/0202009. Hutchinson, E.B., Avram, A.V., Irfanoglu, M.O., Koay, C.G., Barnett, A.S., Komlosh € ME, M.E., Ozarslan, E., Schwerin, S.C., Juliano, S.L., Pierpaoli, C., 2017. Analysis of the effects of noise, DWI sampling, and value of assumed parameters in diffusion MRI models. Magn. Reson. Med. 78 (5), 1767–1780. https://doi.org/10.1002/ mrm.26575. Jenkinson, M., Beckmann, C.F., Behrens, T.E., Woolrich, M.W., Smith, S.M., 2012. FSL NeuroImage 62 (2), 782–790. https://doi.org/10.1016/j.neuroimage.2011.09.015. Jones, D.K., Horsfield, M.A., Simmons, A., 1999. Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging. Magn. Reson. Med. 42 (3), 515–525, 10.1002/(SICI)1522-2594(199909)42:3<515::AIDMRM14>3.0.CO;2-Q. Jones, D.K., Kn€ osche, T.R., Turner, R., 2012. White matter integrity, fiber count, and other fallacies: the do's and don’ts of diffusion MRI. Neuroimage 73, 239–254. https://doi.org/10.1016/j.neuroimage.2012.06.081. Koay, C.G., 2011. A simple scheme for generating nearly uniform distribution of antipodally symmetric points on the unit sphere. J. Comput. Sci. 2 (4), 377–381. https://doi.org/10.1016/j.jocs.2011.06.007. € Koay, G.C., Ozarslan, E., Johnson, K.M., Meyerand, M.E., 2012. Sparse and optimal acquisition design for diffusion MRI and beyond. Med. Phys. 39 (5), 2499–2511. https://doi.org/10.1118/1.3700166. Mitchell, M., 1998. An Introduction to Genetic Algorithms. MIT Press, Cambridge. Mollink, J., Kleinnijenhuis, M., van Cappellen van Walsum, A.-M., Sotiropoulos, S.N., Cottaar, M., Mirfin, C., Heinrich, M.P., Jenkinson, M., Pallebage-Gamarallage, M., Ansorge, O., Jbabdi, S., Miller, K.L., 2017. Evaluating fibre orientation dispersion in white matter: comparison of diffusion MRI, histology and polarized light imaging. Neuroimage 157, 561–574. https://doi.org/10.1016/j.neuroimage.2017.06.001. € _ Ozarslan, E., Koay, C.G., Shepherd, T.M., Komlosh, M.E., Irfano glu, M.O., Pierpaoli, C., Basser, P.J., 2013. Mean apparent propagator (MAP) MRI: a novel diffusion imaging method for mapping tissue microstructure. Neuroimage 78 (1), 16–32. https:// doi.org/10.1016/j.neuroimage.2013.04.016. Pierpaoli, C., Walker, L., Irfanoglu, M.O., Barnett, A.S., Basser, P.J., Chang, L.-C., Koay, C., Pajevic, S., Rohde, G., Sarlls, J.E., Wu, M., 2010. TORTOISE: an integrated software package for processing of diffusion MRI data. Proc. Intl. Soc. Mag. Reson. Med. 1597. Poot, D.H.J., den Dekker, A.J., Achten, E., Verhoye, M., Sijbers, J., 2010. Optimal experimental design for diffusion kurtosis imaging. IEEE Trans. Med. Imaging 29 (3), 819–829. https://doi.org/10.1109/TMI.2009.2037915. Rakshit, P., Konar, A., Das, S., 2017. Noisy evolutionary optimization algorithms – a comprehensive survey. Swarm Evol. Comput. 33, 18–45. https://doi.org/10.1016/ j.swevo.2016.09.002. Ronen, I., Budde, M., Ercan, E., Annese, J., Techawiboonwong, A., Webb, A., 2014. Microstructural organization of axons in the human corpus callosum quantified by diffusion-weighted magnetic resonance spectroscopy of N-acetylaspartate and postmortem histology. Brain Struct. Funct. 219 (5), 1773–1785. https://doi.org/ 10.1007/s00429-013-0600-0. Smith, S.M., 2002. Fast robust automated brain extraction. Hum. Brain Mapp. 17 (3), 143–155. https://doi.org/10.1002/hbm.10062.
References Andersson, J.L.R., Sotiropoulos, S.N., 2016. An integrated approach to correction for offresonance effects and subject movement in diffusion MR imaging. Neuroimage 125, 1063–1078. https://doi.org/10.1016/j.neuroimage.2015.10.019. Andr e, E.D., Grinberg, F., Farrher, E., Maximov, I.I., Shah, N.J., Meyer, C., Jaspar, M., Muto, V., Phillips, C., Balteau, E., 2014. Influence of noise correction on intra- and inter-subject variability of quantitative metrics in diffusion kurtosis imaging. PLoS One 9 (4), e94531. https://doi.org/10.1371/journal.pone.0094531. Assaf, Y., Basser, P.J., 2005. Composite hindered and restricted model of diffusion (CHARMED) MR imaging of the human brain. Neuroimage 27 (1), 48–58. https:// doi.org/10.1016/j.neuroimage.2005.03.042. € Avram, A.V., Sarlls, J.E., Barnett, A.S., Ozarslan, E., Thomas, C., Irfanoglu, M.O., Hutchinson, E., Pierpaoli, C., Basser, P.J., 2016. Clinical feasibility of using mean apparent propagator (MAP) MRI to characterize brain tissue microstructure. Neuroimage 127, 422–434. https://doi.org/10.1016/j.neuroimage.2015.11.027. Back, T., 1994. Selective pressure in evolutionary algorithms: a characterization of selection mechanisms. Proc. First IEEE Conf. Evol. Comput. 57–62. https://doi.org/ 10.1109/ICEC.1994.350042. Baker, J.E., 1987. Reducing bias and inefficiency in the selection algorithm. Proc. Second Int. Conf. Genet. Algorithm. Genet. Algorithm. Appl. 14–21. Baluja, S., Caruana, R., 1995. Removing the genetics from the standard genetic algorithm. Proc. Mach. Learn. 38–46. https://doi.org/10.1016/B978-1-55860-377-6.50014-1. Batchelor, P.G., Atkinson, D., Hill, D.L.G., Calamante, F., Connelly, A., 2003. Anisotropic noise propagation in diffusion tensor MRI sampling schemes. Magn. Reson. Med. 49 (6), 1143–1151. https://doi.org/10.1002/mrm.10491. Beyer, H.-G., 2000. Evolutionary algorithms in noisy environments: theoretical issues and guidance for practice. Comput. Methods Appl. Mech. Eng. 186 (2–4), 239–267. https://doi.org/10.1016/S0045-7825(99)00386-2. Budde, M.D., Annese, J., 2013. Quantification of anisotropy and fiber orientation in human brain histological sections. Front. Integr. Neurosci. 7 (3) https://doi.org/ 10.3389/fnint.2013.00003. Caruyer, E., Deriche, R., 2012. A computational framework for experimental design in diffusion MRI, CDMRI - MICCAI. In: Workshop Comput. Diffus. MRI. Nice, France, hal-00747700. Caruyer, E., Lenglet, C., Sapiro, G., Deriche, R., 2013. Design of multishell sampling schemes with uniform coverage in diffusion MRI. Magn. Reson. Med. 69 (6), 1534–1540. https://doi.org/10.1002/mrm.24736. De Jong, K.A., 1975. Analysis of the Behavior of a Class of Genetic Adaptive Systems. University of Michigan. Descoteaux, M., Deriche, R., Le Bihan, D., Mangin, J.F., Poupon, C., 2011. Multiple q-shell diffusion propagator imaging. Med. Image Anal. 15 (4), 603–621. https://doi.org/ 10.1016/j.media.2010.07.001. Diaz-Gomez, P.A., Hougen, D.F., 2007. Initial population for genetic algorithms: a metric approach. Proc. 2007 Int. Conf. Genet. Evol. Method. 43–49. Dubois, J., Poupon, C., Lethimonnier, F., Le Bihan, D., 2006. Optimized diffusion gradient orientation schemes for corrupted clinical DTI data sets. Magn. Reson. Mater. Phy. 19 (3), 134–143. https://doi.org/10.1007/s10334-006-0036-0. Eshelman, L.J., 1990. The CHC Adaptive Search Algorithm: how to have safe search when engaging in nontraditional genetic recombination. Found. Genet. Algorithms 1, 265–283. https://doi.org/10.1016/B978-0-08-050684-5.50020-3. Fick, R.H.J., Wassermann, D., Sanguinetti, G., Deriche, R., 2015. Laplacian-regularized MAP-MRI: improving axonal caliber estimation. Proc. IEEE Int. Symp. Biomed. Imag. https://doi.org/10.1109/ISBI.2015.7164084. Fitzpatrick, J.M., Grefenstette, J.J., 1988. Genetic algorithms in noisy environments. Mach. Learn. 3 (2–3), 101–120. https://doi.org/10.1007/BF00113893.
243
D.V. Olson et al.
NeuroImage 199 (2019) 237–244 Whitley, D., 1994. A genetic algorithm tutorial. Stat. Comput. 4 (2), 65–85. https:// doi.org/10.1007/BF00175354. Wu, Y.-C., Alexander, A.L., 2007. Hybrid diffusion imaging. Neuroimage 36 (3), 617–629. https://doi.org/10.1016/j.neuroimage.2007.02.050. Ye, W., Portnoy, S., Entezari, A., Blackband, S.J., Vemuri, B.C., 2012. An efficient interlaced multi-shell sampling scheme for reconstruction of diffusion propagators. IEEE Trans. Med. Imaging 31 (5), 1043–1050. https://doi.org/10.1109/ TMI.2012.2184551.
Spears, W.M., De Jong, K.A., 1990. An analysis of multi-point crossover. Found. Genet. Algorithms 1, 301–315. https://doi.org/10.1016/B978-0-08-050684-5.50022-7. Sprenger, T., Sperl, J.I., Fernandez, B., Golkov, V., Eidner, I., S€amann, P.G., Czisch, M., Tan, E.T., Hardy, C.J., Marinelli, L., 2016. Bias and precision analysis of diffusional kurtosis imaging for different acquisition schemes. Magn. Reson. Med. 76 (6), 1684–1696. https://doi.org/10.1002/mrm.26008. Syswerda, G., 1989. Uniform crossover in genetic algorithms. Proc. Int. Conf. Genet. Alg. 2–9.
244
Contents lists available at ScienceDirect
NeuroImage journal homepage: www.elsevier.com/locate/neuroimage
Optimization of q-space sampling for mean apparent propagator MRI metrics using a genetic algorithm Daniel V. Olson a, e, *, Volkan E. Arpinar b, d, L. Tugan Muftuler c, d a
Department of Biophysics, Medical College of Wisconsin, Milwaukee, WI, USA Department of Radiology, Medical College of Wisconsin, Milwaukee, WI, USA Department of Neurosurgery, Medical College of Wisconsin, Milwaukee, WI, USA d Center for Imaging Research, Medical College of Wisconsin, Milwaukee, WI, USA e Magnetic Resonance Imaging, GE Healthcare, Waukesha, WI, USA b c
A R T I C L E I N F O
A B S T R A C T
Keywords: Mean apparent propagator (MAP) Optimization q-space sampling Genetic algorithm
Mean Apparent Propagator (MAP) MRI is a recently introduced technique to estimate the diffusion probability density function (PDF) robustly. Using the estimated PDF, MAP MRI then calculates zero-displacement and nonGaussianity metrics, which might better characterize tissue microstructure compared to diffusion tensor imaging or diffusion kurtosis imaging. However, intensive q-space sampling required for MAP MRI limits its widespread adoption. A reduced q-space sampling scheme that maintains the accuracy of the derived metrics would make it more practical. A heuristic approach for acquiring MAP MRI with fewer q-space samples has been introduced earlier with scan duration of less than 10 minutes. However, the sampling scheme was not optimized systematically to preserve the accuracy of the model metrics. In this work, a genetic algorithm is implemented to determine optimal q-space subsampling schemes for MAP MRI that will keep total scan time under 10 min. Results show that the metrics derived from the optimized schemes more closely match those computed from the full set, especially in dense fiber tracts such as the corpus callosum.
1. Introduction Diffusion MRI (dMRI) is a powerful technique capable of probing tissue microstructure in vivo at a resolution that is much higher than the imaging resolution. Diffusion MRI achieves this by indirectly estimating the net displacement of water molecules in the biological environment. Effective diffusion modeling necessitates a series of diffusion-weighted images (DWI). However, acquiring a DWI set with sufficient diffusionencoding vectors (b-vectors) for advanced modeling is time-intensive. This makes advanced dMRI techniques impractical when working with certain patient populations and unfeasible for clinical use. The distribution of acquired b-vectors affects robustness to sample orientation and noise. At the same time, model fitting can be overdetermined such that a subset of DWI might yield similar results as the full set (Dubois et al., 2006). An ideal q-space sampling scheme balances robustness with minimizing redundancy in the measurements. For the most part, q-space sampling schemes have been designed to approximate uniform angular coverage on the unit sphere. Jones et al. formulated the optimization as a numerical minimization of electrostatic
charge distributions on a sphere (Jones et al., 1999). For 30 directions, this electrostatic solution approximates the icosahedral group which is rotationally invariant (Batchelor et al., 2003; Hasan and Narayana, 2003). An analytical model for distributing antipodally symmetric points on a sphere has been developed, which closely matches the electrostatic model (Koay, 2011). Extending the optimization to multiple shells is required for advanced diffusion models, as single-shell acquisitions are innately limited in the radial direction in q-space. Standard multi-shell acquisitions have uniformly distributed sample points on each shell. Recent works have also sought to maximize the angular incoherence between q-space samples taken across shells (Koay et al., 2012; Caruyer and Deriche, 2012; Caruyer et al., 2013). Generally, more samples are acquired on the high b-value shells (Wu and Alexander, 2007; Assaf and Basser, 2005; € Ozarslan et al., 2013). This results in more uniform sampling but also denser sampling where the diffusion signals are highly attenuated, which has raised concerns that the rectified noise floor will influence modeling (Andre et al., 2014; Hutchinson et al., 2017). Optimizing q-space sampling for a specific diffusion model can
* Corresponding author. Department of Biophysics, Medical College of Wisconsin, Milwaukee, WI, USA. E-mail address: [email protected] (D.V. Olson). https://doi.org/10.1016/j.neuroimage.2019.05.078 Received 4 March 2019; Received in revised form 28 April 2019; Accepted 29 May 2019 Available online 1 June 2019 1053-8119/© 2019 Elsevier Inc. All rights reserved.
D.V. Olson et al.
NeuroImage 199 (2019) 237–244
produce more accurate parameter estimates than general optimization. An example is the optimization for Diffusion Kurtosis Imaging (DKI) performed by Poot et al. (2010), which used simulated annealing as its search method. The optimal q-space sampling did not conform to a multi-shell distribution. Instead, the b-values had a varied radial distribution with a substantial portion of the gradients at the maximum allowed b-value of 2800 s/mm2. The b-vectors were approximately isotropically distributed on the unit sphere. A separate multi-shell scheme from that work performed best in an independent study (Sprenger et al., 2016). This suggests that a departure from multi-shell acquisitions to q-space sampling schemes of higher radial resolution could better characterize the diffusion probability density function. A recently proposed diffusion model called Mean Apparent Propagator (MAP) MRI is a mathematical framework connecting the diffusion€ weighted imaging q-space with the molecular displacements (Ozarslan et al., 2013). The basis functions are eigenfunctions of the quantum-mechanical analog of the simple harmonic oscillator Hamiltonian. The MAP MRI model yields novel metrics to characterize the diffusion properties. Derivation of the metrics from the harmonic co€ efficients is demonstrated in Ozarslan et al. (2013). While MAP MRI may offer better characterization of diffusion in complex tissue structures compared to DTI or DKI, an intensive acquisition impedes its widespread adoption, both clinically and in research. € The original MAP acquisition scheme proposed in Ozarslan et al. (2013) consisted of seven b-values of b ¼ {200, 800, 1800, 3200, 5000, 7200, 9800} s/mm2 along 5, 14, 32, 56, 87, 125, and 170 directions, respectively. An excised marmoset brain was scanned, so there were no concerns about long acquisition time or motion. A human adaptation proposed by Avram et al. (2016) reduced the acquisition time to 72 min. The sampling scheme consisted of six b-values b ¼ {1000, 2000, 3000, 4000, 5000, 6000} s/mm2 along 23, 43, 83, 131, 148, and 256 directions, respectively. A heuristically-derived q-space subsampling scheme was also proposed in that paper for clinical feasibility, which further reduced the number of directions per shell to 4, 7, 11, 17, 23, and 31 directions for an acquisition of about 10 min. In this work, we describe an optimization framework implemented for MAP MRI. A genetic algorithm is employed, which determined clinically feasible subsets of q-space sampling schemes that minimize the mean-squared error (MSE) relative to the MAP metrics produced by the full q-space. The optimized q-space sampling schemes were validated using both a digital phantom and also in vivo data. In the next section, genetic algorithms are reviewed first to provide a context for our implementation.
Fig. 1. Overview of genetic algorithm.
should refer to the citations provided. The GA consisted of four parts: initialization, fitness evaluation, selection, and breeding. Initialization is the construction of the first population. For the initial population, individuals were assigned random subsamples of b-vectors from the full set without prohibition of duplicates. Future generations are determined through iterations of evaluation, selection, and breeding. Each individual was then evaluated with MAP fitting, and the MSE was calculated relative to the true value of the metric from the full data set, which was similar to € the 489 directions scheme used in Ozarslan et al. (2013). The selection method was Stochastic Universal Sampling (SUS), a fitness-proportionate method that minimizes deviation of the selections from their reproduction expectation values (Baker, 1987; Mitchell, 1998). Fitness values were adaptively scaled with sigma scaling (Forrest and Mitchell, 1993). In populations with high fitness variance, sigma scaling limits the impact of individuals with outlying fitness values, as would occur in early generations. Conversely, low fitness variance amplifies the differences between individuals in later generations to maintain a relatively constant selection pressure over the duration of the algorithm (Back, 1994). A final selection concept is elitism, where a small number of the best individuals are guaranteed survival (De Jong, 1975). GA with elitism have been shown to outperform those without it in almost every trial (Baluja and Caruana, 1995). Breeding is comprised of the genetic operators of crossover and mutation. Ninety-eight percent of the offspring were generated through crossover. In uniform crossover, each b-vector is inherited independently from any other for maximal disruption without loss of b-vectors (Syswerda, 1989). A disruptive crossover operator is desired with small population sizes to more widely sample the solution space and avoid the tendency towards population homogeneity (Spears and De Jong, 1990; Eshelman, 1990). Mutations promote a greater traversal of solution space by a random walk. Without mutation, only b-vectors that were present in the initial population would exist in the simulations. Mutation was defined as exchanging a single b-vector for another in the full set. with a probability of 0.01. Mutation rate should be kept low to allow for the crossover operator to determine better q-space sampling schemes. Evaluation, selection, and breeding are repeated until stopping criteria are satisfied. The maximum number of generations was set to 250 to allow for sufficient evolutions while preventing overfitting of the signals. Additional stopping criteria were a stall time of 75 generations as well as convergence tolerances of 106 for non-Gaussianity (NG) and 0.01 for return-to-origin probability (RTOP). In this simulation, convergence was defined as the difference of the lowest individual MSE and the average MSE of the population being below the tolerance. These values correspond to approximately 0.1% error for the two metrics.
2. Methods 2.1. Genetic algorithms A genetic algorithm (GA) is a metaheuristic global optimization solver that utilizes evolutionary biology principles to reach a desired solution (Holland, 1973). GA are population-based solvers that make few assumptions about the problem and thus are considered “weak methods” (Whitley, 1994). Beneficial features are selected and recombined, approaching an optimal tuning over iterations. Even in noisy environments, GA are capable of approaching a global minimum, unlike most other optimizers (Fitzpatrick and Grefenstette, 1988; Beyer, 2000; Rakshit et al., 2017). Although a global optimum is not guaranteed, heuristic solvers can find an adequate solution in a reasonable amount of time. Consider the case of selecting an optimal subset of 95 diffusion-encoding gradient vectors from a collection of 489, which was used in the original € paper by Ozarslan et al. (2013). The binomial coefficient is on the order of 10103 possible combinations. Further, acceptable sampling schemes are sparse in solution space. Therefore, exhaustive and random-walk searches are not feasible. An overview of the GA is presented in Fig. 1, and a glossary of GA terminology is provided in Appendix B. For further details, the readers 238
D.V. Olson et al.
NeuroImage 199 (2019) 237–244
calculated from the MAP coefficients based on the respective q-space sampling schemes for the full set and each subsampled set. The GARTOP and GANG sampling schemes were taken directly from the output of the optimization without modification. For comparison, the reduced sampling scheme proposed by Avram et al. (2016) was also included in the analyses (referred to as Avram scheme). The diffusion-encoding directions for this set were selected through a Monte Carlo farthest neighbor search to approximate a uniformly distributed subset on each shell, but b-vectors were not optimized across shells to maximize angular incoherence on the unit sphere. A subsequent MAP fitting was applied to these simulated signals to calculate the zero-displacement and non-Gaussianity metrics. MSE were calculated for each subsampled set's metrics relative to the corresponding metrics from the full set.
2.2. Acquisition and processing of MAP MRI data for optimization In vivo data were acquired from a healthy volunteer (male, 26 years). The study was approved by the Institutional Review Board, and written informed consent was obtained from the participant. Data were collected on a General Electric (GE) SIGNA Premier 3 T MRI scanner with a 48channel receive head coil. Diffusion weighted images (DWI) were acquired using a single-shot SE-EPI sequence. The imaging volume was 20 slices with an acquisition matrix of 100 100. Voxel dimensions were 2.2 2.2 2 mm. The echo time (TE) was 86.8 ms, and the repetition time (TR) was 4200 ms. Bandwidth was 1953 Hz/pixel. Parallel imaging (ASSET) in-plane acceleration factor of 2 was used to reduce susceptibility artifacts, but no multiband acceleration was used. Acquisition time was 35 min. The data set consisted of 10 interspersed, non-diffusionweighted reference images and 489 DWI. The DWI were sampled in six shells with b ¼ {1000, 2000, 3000, 4000, 5000, 6000} s/mm2 along 19, 32, 56, 87, 125, and 170 directions, respectively. In each shell, diffusionencoding directions were uniformly distributed over the unit sphere (Jones et al., 1999). However, the angular incoherence of the b-vectors was not maximized across the shells. The maximum diffusion gradient amplitude was 69 mT/m with pulse width δ of 22.8 ms and a separation Δ of 41.3 ms. This data set will be called the full data set throughout the rest of this paper. Postprocessing of the data set included EPI phase correction on the scanner. Offline processing included motion and eddy current corrections using the EDDY function of the FSL software (Jenkinson et al., 2012; Andersson and Sotiropoulos, 2016). A brain mask was calculated with BET (Smith, 2002). Then, MAP fitting was performed using the DIPY Toolbox (Garyfallidis et al., 2014). The MAP was fit with radial order of 6. No denoising was applied to the data prior to MAP fitting, but a positivity constraint and Laplacian regularization with a weighting of 0.05 were applied during fitting (Fick et al., 2015). Those settings were recommended for robust fitting by the author of the DIPY MAPMRI toolbox. An initial FA map was also generated from the data set to identify white matter voxels. This FA map was thresholded at FA>0.3 to generate a white matter mask.
2.5. Validation of genetic algorithm results using emulation by down sampling To quantify the accuracy of the MAP metrics for each subsampling scheme in vivo, a similar analysis was performed on the acquired DWI. To emulate separate acquisitions, the full data set was downsampled to incorporate each of the subsampling schemes. The subsampled schemes were not acquired separately to minimize the variability between the sets introduced by subject motion and other physiological noise sources. From the white matter (FA > 0.3), voxels that were not included in the GA optimization were used to test the accuracy of the reduced q-space sampling schemes relative to the full set. MAP model fitting and MSE evaluations were performed in the same manner as with the numerical simulation. The distribution of RTOP and NG metrics are presented in Fig. 2. The distribution closely matched the values used in the GA optimization. 2.6. Second in vivo experiment for independent validation To test the generalizability of the proposed q-space sampling schemes, the same acquisitions were also performed on another healthy volunteer (female, 26 years), whose data was not used in optimization. Written informed consent was obtained from the participant before the study. Data were collected using identical scanning parameters as the first acquisition. The full MAP MRI data set was acquired as before, and the same postprocessing and MAP fitting pipelines were applied, except no FA masking was used.
2.3. Genetic algorithm optimization The genetic algorithm was developed in Python 2.7.14 to interface with the MAP functions of the DIPY Toolbox. Two optimizations were performed separately for return-to-origin probability and nonGaussianity. The resulting q-space sampling schemes are denoted GARTOP and GANG. Diffusion signals were measured from 400 voxels in the white matter. Poot et al. (2010) determined that white matter could be adequately represented with 400 voxels. These voxels were randomly selected from within the FA mask. In determining population size, there is a trade-off between sufficient population diversity to sample solution space and computational resources for more generations. In this optimization, the MAP model fitting is computationally intensive. Diversity has been shown to become independent of the individual size with population sizes of about 200 (Diaz-Gomez and Hougen, 2007). There were 200 individuals in each population for the GA, which consisted of 95 b-vectors from the full set. This corresponds to approximately a 10-min acquisition. The GA evaluation function was the MSE relative to the full set's metrics of either RTOP or NG. Selection was fitness-proportionate with Stochastic Uniform Sampling, as well as sigma scaling to maintain selection pressure. Uniform crossover was used to increase mixing. Elite survival rate, crossover rate, and mutation rate were 0.02, 0.98, and 0.01, respectively.
3. Results Signal-to-Noise Ratios (SNR) were estimated on diffusion data that was used in the optimization. Since the signals in DWI are dependent on the diffusion-encoding direction, SNR can vary substantially depending on the relative direction of the diffusion-encoding and tissue structure. In order to determine the worst-case SNR, the signal should be chosen in a region whose fiber direction corresponds with the diffusion-encoding direction
2.4. Numerical simulation with a digital brain MAP MRI phantom
Fig. 2. Distribution of RTOP1/3 and NG values for simulations. Training set (blue) consisted of 400 random voxels from the white matter (FA > 0.3) and was used in the genetic algorithm optimization. Test set (red) consisted of the remaining voxels from the white matter. The test set was used in the quantitative validation of q-space sampling schemes.
The reduced sampling schemes were first tested using a digital brain phantom. To generate the digital phantom, an initial MAP fitting was performed on the full set. Then, diffusion-weighted signals were 239
D.V. Olson et al.
NeuroImage 199 (2019) 237–244
(Descoteaux et al., 2011; Jones et al., 2012). This will show the most attenuated signal for a given b-value and more clearly indicate image quality. To calculate SNR values, regions of interest were manually drawn in the center of the splenium of the corpus callosum (which corresponded with the Gx direction) and the background. SNR levels were 111, 38.6, 18.0, 12.3, 8.6, 8.2, and 8.0 for respective b-values of 0–6000 s/mm2. DWI at b ¼ 6000 s/mm2 retained sufficient signal such that tissue anatomy was still clearly delineated. Interspersed non-diffusion-weighted images were acquired every 50 diffusion-weighted images. Motion artifacts were not observed in the DWI. The average temporal signal-to-noise ratio estimated from the non-diffusion-weighted images was 22.2 (Fig. 3). 3.1. Genetic algorithm and the optimized q-space sampling schemes The GA schemes acquired fewer directions on the middle shells than the innermost and outmost shells (Fig. 4). The final result shows that the number of samples from the b ¼ 1000 s/mm2 shell increased for both GARTOP and GANG schemes compared to the Avram scheme. This finding suggests that acquiring more samples on shells with lower b-values may be needed for robust fitting. On the other hand, collecting more directions on higher b-value shells has been generally recommended for more accurate MAP fitting. Still, no evidence of this requirement has been demonstrated to date. Our GA-based optimization clearly shows that more samples at high b-values are essential. In the original paper € that introduced the MAP MRI technique, Ozarslan et al. (2013) did ex vivo scans and went up to b ¼ 9800 s/mm2. In vivo human scans for MAP MRI typically goes up to b ¼ 6000 s/mm2 for practical reasons. Therefore, we also limited the outer shell to b ¼ 6000 s/mm2, which might be adequate if the radial order of the MAP MRI series expansion is kept at 6. Another consideration is the angular q-space coverage. Each shell of the full set contained isotropically distributed b-vectors determined by the electrostatic charge model, but diffusion gradient directions were not optimized across shells. The Avram scheme was determined by a Monte Carlo farthest neighbor search that found a subset of b-vectors within each shell with the highest angular difference. The distributions of bvectors on the unit sphere for each of the optimized sampling schemes tested here were approximately isotropic.
Fig. 4. Distribution of b-values for sampling schemes. Avram scheme (blue) increases the number of b-vectors per shell with b-value. The samplings optimized for RTOP1/3 (green) and NG (yellow) acquire more directions on the b ¼ 1000 s/mm2 shell, fewer directions in the middle shells, and a high number on the outer b ¼ 6000 s/mm2 shell.
were computed (Fig. 5). The ground truth for the MSE calculations were the metrics calculated from the respective full q-space sampling set. GA optimization produced schemes with lower errors for all metrics, compared to the heuristic subsampling scheme. GARTOP outperformed the other schemes for the zero-displacement metrics, as well as slightly outperforming GANG in the non-Gaussianity metrics. Both GA optimization schemes performed exceptionally well, particularly for nonGaussianity measures. 3.3. In vivo experiments for metric validation Validation was performed using diffusion data from the two subjects. For the first subject, only white matter voxels that were not used in optimization were used in the analyses to reduce bias. For independent verification of optimized sampling schemes, all white matter voxels from the second subject were used. Since those voxels were not used for optimization, a bias would not be expected. Both GA schemes reduced the MSE for all MAP metrics compared to the Avram scheme in white matter (Fig. 6). With in vivo data, GARTOP resulted in the lowest MSE values for the zero-displacement metrics, while GANG had the best performance for every non-Gaussianity metric. This white matter analysis
3.2. Numerical simulation using a digital MAP MRI brain phantom For each of the reduced sampling schemes, MSE of the MAP metrics
Fig. 5. Mean-squared error (MSE) of MAP metrics in white matter from numerical simulation. Avram q-space sampling (blue), genetic algorithm optimized for RTOP1/3 (green), genetic algorithm optimized for NG (yellow). MSE are relative to metrics calculated from the full q-space sampling set. Diffusion signals were calculated from MAP model coefficients from voxels in the white matter (FA > 0.3). GA schemes reduce MSE compared with the Avram scheme. The sampling optimized for RTOP provided the best overall results. Units of zero-displacement metrics are mm1.
Fig. 3. Temporal SNR of non-diffusion-weighted (b ¼ 0 s/mm2) volumes. Temporal SNR was calculated as the mean signal across the volumes normalized by the standard deviation. 240
D.V. Olson et al.
NeuroImage 199 (2019) 237–244
closely matched the predictions of the optimization. White matter MSE trends were similar between the two subjects, which implies generalizability of the optimized schemes. There were pronounced improvements in the non-Gaussianity metrics for both of the optimized sampling schemes. While the relative differences between the optimized sampling schemes and the Avram scheme in zerodisplacement probabilities (RTOP1/3, RTAP1/2 and RTPP) were not as high as the first set, there were still notable decreases in MSE for all three metrics. This analysis was repeated for the whole brain, since MAP values in structures such as deep brain gray matter could also be informative for various diseases. In the whole brain MSE calculations, voxels with RTOP1/3 values below 20 were excluded. These voxels likely contained CSF and were predominantly from the ventricles. The improvement in the non-Gaussianity metrics of the GA optimized schemes was less pronounced in whole brain MSE, yet it is still notable (Fig. 7). The trends of GARTOP performing best in the zero-displacement metrics and GANG in the non-Gaussianity measures are evident in whole brain MSE as well. RTAP1/2 errors were approximately halved in the whole brain analysis relative to the white matter analysis. This is likely an averaging effect from including gray matter, which has lower RTAP1/2 values than white matter. As a reminder, the optimization was run with an ensemble of voxels selected from the white matter. If a study requires more emphasis on, for instance, deep gray matter nuclei, the optimization could be run with additional voxels included from those regions. Qualitatively, each of the reduced schemes produced sufficiently accurate RTOP1/3 fits (Fig. 8, top). Slight improvement can be seen in the white matter with the GA sampling schemes. Similar to results found in Avram et al. (2016), fitting differences were observed around the edges of the ventricles, possibly due to partial volume effects. The low intensity areas observed in all RTOP1/3 maps near the gray matter could be CSF. The effects were most pronounced in the map from the full set. These regions were excluded from the MSE analysis with an RTOP1/3 > 20 threshold. As evidenced in the NG difference maps (Fig. 8, bottom), the GANG scheme most closely replicated the NG of the full set. GARTOP was relatively accurate in the white matter but slightly overestimated the gray matter NG. The Avram scheme performed similar to GARTOP except in the corpus callosum. Overall, both GARTOP and GANG performed well in highly organized fibers such as the splenium and genu producing values much closer to the full set.
Fig. 7. Mean-squared error of MAP metrics in whole brain. Avram sampling (blue), genetic algorithm optimized for RTOP1/3 (green), genetic algorithm optimized for NG (yellow). Each pair of bars represent results from the two subjects. For each pair, the left bar is from the first subject and the right bar is from the second subject. MSE are relative to metrics calculated from the full qspace sampling set. The sampling optimized for NG provided the best overall results. However, the GARTOP scheme did perform better for the zerodisplacement probabilities. Units of zero-displacement metrics are mm1.
The observed NG was slightly lower in the splenium and genu relative to other white matter tracts. In those fiber tracts, we expected to see high NG values as have usually been observed in DKI. Our initial interpretation was that Gibbs ringing might be affecting the NG estimates. Upon further inspection, however, we observed that the lower NG values
Fig. 6. Mean-squared error (MSE) of MAP metrics in white matter. Avram qspace sampling (blue), genetic algorithm optimized for RTOP1/3 (green), genetic algorithm optimized for NG (yellow). Each pair of bars represent results from the two subjects. For each pair, the left bar is from the first subject and the right bar is from the second subject. MSE are relative to metrics calculated from the full q-space sampling set. Diffusion signals were from voxels in the white matter (FA > 0.3). GA schemes reduce MSE compared with the Avram scheme. The sampling optimized for NG provided the best overall results. Units of zerodisplacement metrics are mm1.
Fig. 8. RTOP1/3 and NG maps from in vivo validation. RTOP1/3 metric maps (first row) and absolute difference images (second row) for each sampling scheme. Each reduced set performed similarly. NG metric maps (third row) and absolute difference images (bottom row) for each sampling scheme. GANG most closely matched the full set, especially in the major white matter tracts. The Avram scheme's largest deviation was in the corpus callosum. 241
D.V. Olson et al.
NeuroImage 199 (2019) 237–244
results. However, the same improvements were seen in the numerical simulation and second data set, which increases confidence that such artifacts did not play a significant role in these findings. While the Monte Carlo farthest-neighbor search to replicate the Avram scheme did fairly well at subsampling uniformly within each shell, it is possible that optimizing the angular incoherence of the b-vectors across shells would have improved fitting. However, in that case, the scheme would no longer be a subset of the full set. For the GA schemes, q-space interpolation could be used to synthesize diffusion signals beyond those acquired in the full set. Interpolation of qspace has been used in previous works including Hybrid Diffusion Imaging (HYDI) (Wu and Alexander, 2007) and a propagator estimation method by Ye et al. (2012). A future direction could be applying the MAP model to interpolated q-space signals to determine an optimal q-space sampling scheme with higher radial sampling density than the six shells acquired in this work. Diffusion signals should be synthesized in advance for a set of b-vectors. Still, more degrees of freedom rapidly expand the solution space, so the GA parameters would need to be adapted. Furthermore, the accuracy of such q-space interpolation needs to be validated before attempting such an approach. Computation time almost entirely depended on the MAP fitting. This was a limiting factor on the population size as well as the number of white matter voxels included. The MAP fitting must be performed for each individual on each voxel, which could be parallelized in future implementations. Beyond the fitting, minimal calculations are required for the algorithm. Fitness values and genetic operators are calculated once per generation, and computation time has negligible dependence on population size. Since this type of optimization is run only once for a particular study, computational time is not a major concern. This study focused on optimization for healthy volunteers. In patient populations where the diffusion properties are dramatically different than healthy tissue, the optimization can be repeated with patient data.
closely followed the whole splenium into the occipital lobe, where Gibbs ringing would not be expected. A similar trend was seen in the genu as well. There is further evidence that this trend actually reflects the underlying anatomy. For instance, the two lateral branches of the genu are known to have lower dispersion than in the midline (Ronen et al., 2014; Budde and Annese, 2013; Mollink et al., 2017). The NG images reflect this anatomical feature well in the full set, GARTOP, and GANG sets, while they are generally obscured in the Avram scheme. Therefore, there is enough evidence that the NG might be accurately reflecting the actual anatomical features with the optimized MAP acquisition schemes. In order to minimize the possibility that the measurements in the splenium and genu might be affected by Gibbs ringing, we reprocessed the data from the first subject using the DIFFPREP tool of TORTOISE, which has an option to remove Gibbs ringing from the DWI data (Pierpaoli et al., 2010). The measurements in the splenium and genu in the NG and RTAP1/2 images did not change after Gibbs ringing removal, supporting the view that the values in those regions reflect the anatomical features, not image artifacts. 4. Discussion The genetic algorithm implemented in this work determined two qspace sampling schemes, which could promote more widespread use of MAP MRI. Further reductions in the number of b-vectors is possible through this method, although accuracy and robustness would need to be carefully considered. It is important to note that this GA optimization merely tries to use fewer q-space samples and replicate the MAP metrics calculated from the full set. It does not find an optimal q-space sampling scheme for MAP MRI, in general. The methods and results presented here is a proof of concept to demonstrate that an optimal q-space sampling can be determined using a systematic approach. The goal was to minimize sampling redundancies without significant loss of accuracy for MAP MRI metric calculations. While GA provided such optimal acquisition schemes, it is possible that other systematic search methods could provide equally good or better results. The goal here was not to find the best performance among various optimization schemes but to show that it is possible to perform a systematic search and find a reduced set of q-space samples that delivered MAP MRI metrics that matched those of the full acquisition. It is possible that some of the voxels used in optimization or testing were affected by various sources of artifacts, such as cardiac pulsation, physiological motion, and eddy currents. Therefore, a digital brain phantom data was generated by calculating the diffusion signals using the MAP coefficients. This way, the effects of acquisition artifacts were mitigated by making the voxel data consistent for each sampling scheme. This simulation confirmed that the GA optimized subsampling schemes produced metrics that closely matched those of the full data set. Additionally, an independent data set was acquired from a different subject in order to test the robustness of the optimized sampling schemes. Data from a different subject with different anatomy and different time-varying artifacts would show that the optimized schemes could be generalized. The optimized schemes outperformed the Avram scheme for both in vivo data sets, indicating that the optimized schemes could be generalized and were not significantly affected by such sources of artifacts. A limitation of this work is that the reduced q-space samples were subsampled from the full set for analyses. This approach was chosen to avoid possible sources of error such as misregistration between scans, motion artifacts, and other effects of post-processing. If each subsampling scheme were acquired separately, the results would be skewed unfairly if the subject moved slightly more in one of the acquisitions compared to the others. It is still possible that one of the subsampling schemes coincidentally selected more of the volumes with artifacts and biased the
5. Conclusion MAP MRI is a novel technique that generates quantitative maps which could be sensitive and specific to various pathological changes in tissues. The metric maps are derived from the diffusion propagator that is estimated from a broader sampling of q-space. As a result, more subtle changes in tissues could be detected when compared to the ellipsoid assumption of DTI or second order effects modeled by kurtosis. It also does not make assumptions about the tissues under investigation. Acquisition time is a major impediment to more widespread adoption of MAP MRI. With reduced q-space sampling, data can be acquired in a clinically feasible time. However, the accuracy of metrics could be compromised if an optimal subsampling scheme was not devised. Among various optimization approaches, genetic algorithms are effective tools for q-space sampling optimization because they efficiently search a vast solution space by mimicking evolutionary biology principles. Here, we demonstrated that the GA approach can determine which diffusion-encodings could be eliminated with minimal influence on the final MAP MRI images to reproduce the metrics of the full set. While we demonstrated this optimization approach for MAP MRI, the methods described here can be adapted to other diffusion MRI techniques that require a large number of diffusion-encoding vectors. Acknowledgements The authors would like to thank Dr. Matthew Budde for helpful suggestions in the interpretation of results and preparation of the manuscript. This work is supported by funding provided by the Daniel M. Soref Charitable Trust.
242
D.V. Olson et al.
NeuroImage 199 (2019) 237–244
Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi.org/10.1016/j.neuroimage.2019.05.078. Appendix B. Glossary of genetic algorithm terminology Term
Definition
Breeding Crossover Crossover Rate Elite Survival Rate Elitism Evaluation Fitness Generation Individual Mutation Mutation Rate Offspring Population Selection
Application of crossover and mutation to produce new offspring Recombination of b-vectors from two individuals to create an offspring Proportion of the subsequent population created through crossover Proportion of the subsequent population created through elitism Copying the best individuals into the next generation Determination of the fitness of each individual Measure of an individual's performance Iteration of the genetic algorithm q-space sampling scheme Random replacement of a b-vector Probability that a given b-vector will be replaced Individual in the population of the subsequent generation Group of individuals in the generation Process of determining the probabilities that each individual will breed
Forrest, S., Mitchell, M., 1993. What makes a problem hard for a genetic algorithm? Some anomalous results and their explanation. Mach. Learn. 13 (2–3), 285–319. https:// doi.org/10.1007/BF00993046. Garyfallidis, E., Brett, M., Amirbekian, B., Rokem, A., van der Walt, S., Descoteaux, M., Nimmo-Smith, I., 2014. DIPY contributors, DIPY, a library for the analysis of diffusion MRI data. Front. Neuroinf. 8 (8) https://doi.org/10.3389/ fninf.2014.00008. Hasan, K.M., Narayana, P.A., 2003. Computation of the fractional anisotropy and mean diffusivity maps without tensor decoding and diagonalization: theoretical analysis and validation. Magn. Reson. Med. 50 (3), 589–598. https://doi.org/10.1002/ mrm.10552. Holland, J.H., 1973. Genetic algorithms and the optimal allocation of trials. SIAM J. Comput. 2 (2), 88–105. https://doi.org/10.1137/0202009. Hutchinson, E.B., Avram, A.V., Irfanoglu, M.O., Koay, C.G., Barnett, A.S., Komlosh € ME, M.E., Ozarslan, E., Schwerin, S.C., Juliano, S.L., Pierpaoli, C., 2017. Analysis of the effects of noise, DWI sampling, and value of assumed parameters in diffusion MRI models. Magn. Reson. Med. 78 (5), 1767–1780. https://doi.org/10.1002/ mrm.26575. Jenkinson, M., Beckmann, C.F., Behrens, T.E., Woolrich, M.W., Smith, S.M., 2012. FSL NeuroImage 62 (2), 782–790. https://doi.org/10.1016/j.neuroimage.2011.09.015. Jones, D.K., Horsfield, M.A., Simmons, A., 1999. Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging. Magn. Reson. Med. 42 (3), 515–525, 10.1002/(SICI)1522-2594(199909)42:3<515::AIDMRM14>3.0.CO;2-Q. Jones, D.K., Kn€ osche, T.R., Turner, R., 2012. White matter integrity, fiber count, and other fallacies: the do's and don’ts of diffusion MRI. Neuroimage 73, 239–254. https://doi.org/10.1016/j.neuroimage.2012.06.081. Koay, C.G., 2011. A simple scheme for generating nearly uniform distribution of antipodally symmetric points on the unit sphere. J. Comput. Sci. 2 (4), 377–381. https://doi.org/10.1016/j.jocs.2011.06.007. € Koay, G.C., Ozarslan, E., Johnson, K.M., Meyerand, M.E., 2012. Sparse and optimal acquisition design for diffusion MRI and beyond. Med. Phys. 39 (5), 2499–2511. https://doi.org/10.1118/1.3700166. Mitchell, M., 1998. An Introduction to Genetic Algorithms. MIT Press, Cambridge. Mollink, J., Kleinnijenhuis, M., van Cappellen van Walsum, A.-M., Sotiropoulos, S.N., Cottaar, M., Mirfin, C., Heinrich, M.P., Jenkinson, M., Pallebage-Gamarallage, M., Ansorge, O., Jbabdi, S., Miller, K.L., 2017. Evaluating fibre orientation dispersion in white matter: comparison of diffusion MRI, histology and polarized light imaging. Neuroimage 157, 561–574. https://doi.org/10.1016/j.neuroimage.2017.06.001. € _ Ozarslan, E., Koay, C.G., Shepherd, T.M., Komlosh, M.E., Irfano glu, M.O., Pierpaoli, C., Basser, P.J., 2013. Mean apparent propagator (MAP) MRI: a novel diffusion imaging method for mapping tissue microstructure. Neuroimage 78 (1), 16–32. https:// doi.org/10.1016/j.neuroimage.2013.04.016. Pierpaoli, C., Walker, L., Irfanoglu, M.O., Barnett, A.S., Basser, P.J., Chang, L.-C., Koay, C., Pajevic, S., Rohde, G., Sarlls, J.E., Wu, M., 2010. TORTOISE: an integrated software package for processing of diffusion MRI data. Proc. Intl. Soc. Mag. Reson. Med. 1597. Poot, D.H.J., den Dekker, A.J., Achten, E., Verhoye, M., Sijbers, J., 2010. Optimal experimental design for diffusion kurtosis imaging. IEEE Trans. Med. Imaging 29 (3), 819–829. https://doi.org/10.1109/TMI.2009.2037915. Rakshit, P., Konar, A., Das, S., 2017. Noisy evolutionary optimization algorithms – a comprehensive survey. Swarm Evol. Comput. 33, 18–45. https://doi.org/10.1016/ j.swevo.2016.09.002. Ronen, I., Budde, M., Ercan, E., Annese, J., Techawiboonwong, A., Webb, A., 2014. Microstructural organization of axons in the human corpus callosum quantified by diffusion-weighted magnetic resonance spectroscopy of N-acetylaspartate and postmortem histology. Brain Struct. Funct. 219 (5), 1773–1785. https://doi.org/ 10.1007/s00429-013-0600-0. Smith, S.M., 2002. Fast robust automated brain extraction. Hum. Brain Mapp. 17 (3), 143–155. https://doi.org/10.1002/hbm.10062.
References Andersson, J.L.R., Sotiropoulos, S.N., 2016. An integrated approach to correction for offresonance effects and subject movement in diffusion MR imaging. Neuroimage 125, 1063–1078. https://doi.org/10.1016/j.neuroimage.2015.10.019. Andr e, E.D., Grinberg, F., Farrher, E., Maximov, I.I., Shah, N.J., Meyer, C., Jaspar, M., Muto, V., Phillips, C., Balteau, E., 2014. Influence of noise correction on intra- and inter-subject variability of quantitative metrics in diffusion kurtosis imaging. PLoS One 9 (4), e94531. https://doi.org/10.1371/journal.pone.0094531. Assaf, Y., Basser, P.J., 2005. Composite hindered and restricted model of diffusion (CHARMED) MR imaging of the human brain. Neuroimage 27 (1), 48–58. https:// doi.org/10.1016/j.neuroimage.2005.03.042. € Avram, A.V., Sarlls, J.E., Barnett, A.S., Ozarslan, E., Thomas, C., Irfanoglu, M.O., Hutchinson, E., Pierpaoli, C., Basser, P.J., 2016. Clinical feasibility of using mean apparent propagator (MAP) MRI to characterize brain tissue microstructure. Neuroimage 127, 422–434. https://doi.org/10.1016/j.neuroimage.2015.11.027. Back, T., 1994. Selective pressure in evolutionary algorithms: a characterization of selection mechanisms. Proc. First IEEE Conf. Evol. Comput. 57–62. https://doi.org/ 10.1109/ICEC.1994.350042. Baker, J.E., 1987. Reducing bias and inefficiency in the selection algorithm. Proc. Second Int. Conf. Genet. Algorithm. Genet. Algorithm. Appl. 14–21. Baluja, S., Caruana, R., 1995. Removing the genetics from the standard genetic algorithm. Proc. Mach. Learn. 38–46. https://doi.org/10.1016/B978-1-55860-377-6.50014-1. Batchelor, P.G., Atkinson, D., Hill, D.L.G., Calamante, F., Connelly, A., 2003. Anisotropic noise propagation in diffusion tensor MRI sampling schemes. Magn. Reson. Med. 49 (6), 1143–1151. https://doi.org/10.1002/mrm.10491. Beyer, H.-G., 2000. Evolutionary algorithms in noisy environments: theoretical issues and guidance for practice. Comput. Methods Appl. Mech. Eng. 186 (2–4), 239–267. https://doi.org/10.1016/S0045-7825(99)00386-2. Budde, M.D., Annese, J., 2013. Quantification of anisotropy and fiber orientation in human brain histological sections. Front. Integr. Neurosci. 7 (3) https://doi.org/ 10.3389/fnint.2013.00003. Caruyer, E., Deriche, R., 2012. A computational framework for experimental design in diffusion MRI, CDMRI - MICCAI. In: Workshop Comput. Diffus. MRI. Nice, France, hal-00747700. Caruyer, E., Lenglet, C., Sapiro, G., Deriche, R., 2013. Design of multishell sampling schemes with uniform coverage in diffusion MRI. Magn. Reson. Med. 69 (6), 1534–1540. https://doi.org/10.1002/mrm.24736. De Jong, K.A., 1975. Analysis of the Behavior of a Class of Genetic Adaptive Systems. University of Michigan. Descoteaux, M., Deriche, R., Le Bihan, D., Mangin, J.F., Poupon, C., 2011. Multiple q-shell diffusion propagator imaging. Med. Image Anal. 15 (4), 603–621. https://doi.org/ 10.1016/j.media.2010.07.001. Diaz-Gomez, P.A., Hougen, D.F., 2007. Initial population for genetic algorithms: a metric approach. Proc. 2007 Int. Conf. Genet. Evol. Method. 43–49. Dubois, J., Poupon, C., Lethimonnier, F., Le Bihan, D., 2006. Optimized diffusion gradient orientation schemes for corrupted clinical DTI data sets. Magn. Reson. Mater. Phy. 19 (3), 134–143. https://doi.org/10.1007/s10334-006-0036-0. Eshelman, L.J., 1990. The CHC Adaptive Search Algorithm: how to have safe search when engaging in nontraditional genetic recombination. Found. Genet. Algorithms 1, 265–283. https://doi.org/10.1016/B978-0-08-050684-5.50020-3. Fick, R.H.J., Wassermann, D., Sanguinetti, G., Deriche, R., 2015. Laplacian-regularized MAP-MRI: improving axonal caliber estimation. Proc. IEEE Int. Symp. Biomed. Imag. https://doi.org/10.1109/ISBI.2015.7164084. Fitzpatrick, J.M., Grefenstette, J.J., 1988. Genetic algorithms in noisy environments. Mach. Learn. 3 (2–3), 101–120. https://doi.org/10.1007/BF00113893.
243
D.V. Olson et al.
NeuroImage 199 (2019) 237–244 Whitley, D., 1994. A genetic algorithm tutorial. Stat. Comput. 4 (2), 65–85. https:// doi.org/10.1007/BF00175354. Wu, Y.-C., Alexander, A.L., 2007. Hybrid diffusion imaging. Neuroimage 36 (3), 617–629. https://doi.org/10.1016/j.neuroimage.2007.02.050. Ye, W., Portnoy, S., Entezari, A., Blackband, S.J., Vemuri, B.C., 2012. An efficient interlaced multi-shell sampling scheme for reconstruction of diffusion propagators. IEEE Trans. Med. Imaging 31 (5), 1043–1050. https://doi.org/10.1109/ TMI.2012.2184551.
Spears, W.M., De Jong, K.A., 1990. An analysis of multi-point crossover. Found. Genet. Algorithms 1, 301–315. https://doi.org/10.1016/B978-0-08-050684-5.50022-7. Sprenger, T., Sperl, J.I., Fernandez, B., Golkov, V., Eidner, I., S€amann, P.G., Czisch, M., Tan, E.T., Hardy, C.J., Marinelli, L., 2016. Bias and precision analysis of diffusional kurtosis imaging for different acquisition schemes. Magn. Reson. Med. 76 (6), 1684–1696. https://doi.org/10.1002/mrm.26008. Syswerda, G., 1989. Uniform crossover in genetic algorithms. Proc. Int. Conf. Genet. Alg. 2–9.
244