Fitting IVIM with Variable Projection and Simplicial Optimization
Fitting multi-exponential models to Diffusion MRI (dMRI) data has always been challenging due to various underlying complexities. In this work, we introduce a novel and robust fitting framework for the standard two-compartment IVIM microstructural model. This framework provides a significant improvement over the existing methods and helps estimate the associated diffusion and perfusion parameters of IVIM in an automatic manner. As a part of this work we provide capabilities to switch between more advanced global optimization methods such as simplicial homology (SH) and differential evolution (DE). Our experiments show that the results obtained from this simultaneous fitting procedure disentangle the model parameters in a reduced subspace. The proposed framework extends the seminal work originated in the MIX framework, with improved procedures for multi-stage fitting. This framework has been made available as an open-source Python implementation and disseminated to the community through the DIPY project.
Variable projection simplicial homology differential evolution model fitting IVIM microstructure perfusion diffusion projection functional optimization
The Intravoxel Incoherent Motion (IVIM) model uses data collected with the pulsed gradient spin echo sequence, often used to measure diffusion-weighted MRI contrasts, to derive information about blood micro-circulation. In IVIM, information about perfusion of blood is estimated as a pseudo-diffusion process from a PGSE experiment, in which several low b-values are added to sensitize the measurement both to diffusion, as well as to perfusion (Bihan, 2019). The IVIM model has found important applications in a wide variety of perfusion-driven imaging methods since its introduction (Bihan et al., 1988). For example, IVIM is applicable to the physiology of brain tumors (Federau et al., 2013) along with monitoring the response of treatment to these tumors (Kim et al., 2013). Over the past twenty years, IVIM has shown promising results in imaging of breast (Iima and Kataoka, 2018), nasopharyngeal (Lai et al., 2017), pancreatic (Kang et al., 2013), and liver (Li et al., 2017) tumors. In addition, the IVIM model can be used in virtual MR elastography (Bihan et al., 2017) and MR angiography (Bihan, 2019).
The IVIM signal model is expressed as a mixture of two exponential components, representing the standard bi-compartmental isotropic signal decay model for perfusion and diffusion. This makes it difficult to fit, for the following reasons: (1) The exponential components in the signal model are non-orthogonal, which makes it hard to project it along the real-axis by taking an integral and (2) The inherent resolution limit in modeling the exponential decay of multi-exponential models (Istratov and Vyvenko, 1999). The present work tackles these challenges by simultaneously fitting the model parameters without an empirical fixed threshold on the b-values. In previous work, a variety of different approaches were taken to fitting the IVIM model (Jalnefjord et al., 2018). However, there has been no consensus about the best method to use (Bihan, 2019). Furthermore, the existing literature on fitting methods underscores the trade-off between precision and bias in model estimation.
In the signal equation of this model, two exponentials are entangled: their linear coefficients represent volume fractions of a diffusion-weighted component and a perfusion-weighted component, and their exponents represent the apparent diffusion coefficient (ADC) and a pseudo-diffusion coefficient of these components (Bihan et al., 1988). Fitting this model can be solved as a simple least squares problem, by minimizing the norm of the residuals between the observed and the predicted values. Unfortunately, because of the challenges mentioned above, fitting the linear coefficients and the exponents simultaneously is a hard nonlinear problem. Several methods have been proposed to analyze and fit multi-exponential functions, such as graphical analysis (peeling method) (Bell, 1965) (Ferrante et al., 1987), Nonlinear Least Squares (NLS) (Byrd et al., 1988), algebraic techniques (Prony’s method) (Osborne and Smyth, 1995) and Pade-Laplace algorithm (Halvorson, 1992).
Here, we propose to use an approach previously used in the MIX framework (Farooq et al., 2016) by improving the global optimization strategy. The MIX framework showed that the Variable Projection approach was beneficial to fit multi-compartment microstructure models for crossings. For the IVIM model at hand, it is particularly useful to separate the treatment of linear parameters (i.e. the perfusion and diffusion fractions) from that of the nonlinear parameters (i.e. the perfusion and diffusion coefficients) via the Variable Projection method (Golub and Pereyra, 2003). The main advantage of using Variable Projection to segregate the linear parameters from the nonlinear ones lies in simultaneous fitting without the need for a fixed threshold to split b-values for perfusion and diffusion components.
While fitting the volume fractions is a simple convex optimization problem, we still need a good estimate to fit the nonlinear exponent parameters. For this purpose, we make use of state-of-the art global optimization methods to get good initial estimates for the entangled exponential functions. We provide tools to efficiently switch between different types of global optimization: simplicial homology (SH) (Endres et al., 2018) (Paulavičius and Žilinskas, 2014) and differential evolution (DE). Our results show that the SH and DE methods for global optimization give similar results, but SH reduces the number of function evaluations required and is therefore faster as compared to DE. The results obtained from either of the global optimizers and convex optimizer give a good initial estimate for the model parameters in their separated subspaces. We then make use of these initial estimates to fit the full IVIM model using the Trust Region Reflective optimizer via nonlinear least squares (Golub and Pereyra, 2003).
2.1 Overview of the fitting framework
The proposed framework using variable projection (Golub and Pereyra, 2003; Farooq et al., 2016) aims to disentangle the linear parameters (here, the associated perfusion and diffusion fractions) by projecting and solving the problem in a reduced subspace of the nonlinear parameters. This is done by taking an epigraphical projection of the measured signal onto the range of the which is assumed to have a closed form solution. The next step is to optimize for which is done via either the Simplicial Homology (SH) or Differential Evolution (DE) algorithm (explained in detail in Section 2.3 and Section 2.4 respectively).
Simplicial homology based topological optimization works by constructing homology groups of the hypersurface of the objective function. The SH optimizer is particularly useful as it is non-stochastic in nature in comparison with evolutional algorithms. The calculated homology groups perform optimization in a derivative-free setting (Rios and Sahinidis, 2012) so that heuristics to switch between the local and global search space are not needed. SH has not only shown to have comparative comparisons with DE and Basinhopping (Endres et al., 2018), but also gives the user a unique capability to track the optimization process via the constructed homology groups. In our hands, SH gives the same results as the DE optimizer with an improved speed of searching the hypersurface of the objective function.
Evolutionary computation is known to be particularly useful for fitting exponential functions. In the proposed framework, DE was chosen as the stochastic optimizer as it is known to outperform other metaheuristic optimization methods, such as Particle Swarm Optimization, Genetic Algorithms (GA), Harmony Search (HS) and Seeker Optimization Algorithms (SOA) (Ismail and Halim, 2017). The particular advantage of using DE as an optimizer is that it gives the same result even if the model parameters are varied which is useful to speed up the search process and to avoid parameter tuning under different settings. Due to low values in the standard deviation of the fit, the results obtained from this approach improve the predictability in fitting. However, our experiments depict that using SH is a better strategy than DE for this type of a problem. Comparisons of SH against DE has been done in greater detail in global optimization literature
(Endres et al., 2018).
After the global optimization step, convex optimization is used to get better initial estimates of the projected volume fractions. We know that the volume fractions are combined as a convex combination (i.e. ). Thus a simple search for these linear parameters in the reduced subspace improves the quality of the fit and helps obtain refined estimates for the respective volume fractions for the last stage. The nonlinear least squares solver is used as the final step for finding the parameter estimates, but now with improved initial estimates of the space. We show that the parameters obtained from this estimation significantly improve the quality of fitting.
2.2 Variable Projection (VarPro) for IVIM model
The volume fraction parameters of IVIM model appear linearly, i.e. as linear combinations of variables that are nonlinearly parameterized (exponential components). Therefore, for every choice of the nonlinear parameters, we can solve for the linear coefficients as a convex optimization problem.
Problem formulation: The IVIM model equation can be written as follows:
Here denotes the perfusion fraction and the fraction associated with diffusion component. From the model setup (Bihan et al., 1988), we know that . and are the nonlinear parameters that are entangled with the volume fractions. denotes the diffusion coefficent of water in the blood. Note that the parameter together are just treated as during fitting as the fraction of is much smaller as compared to the tissue water content.
To take the epigraphical projection of the measured signal, let us define the signal model as: and (isotropic models for the exponential decay of each of the individual components). Now we can define this set of nonlinear functions as our atom functions of a finite dictionary (Mitra and Bhatia, 2014) that we will combine as a linear sum of the volume fractions (). We can set up our dictionary with the atom functions (here, signal models for perfusion and diffusion) as:
Now, as mentioned (Golub and Pereyra, 2003), we will project out the volume fractions with the help of a projector function. This is a Moore-Penrose inverse of the nonlinear components of our model defined in equation (2). We can write the projection function as:
We can make use of the measured signal and the projector function to now project out the linear components onto the range of as: . Using the above projected component for the volume fraction, we can reformulate the objective function of the IVIM model as in equation (3) where we only minimize over in a reduced subspace using the projector function described in equation (4).
2.3 Simplicial Homology Global Optimization
The Simplicial Homology (SH) optimizer is a topographic optimizer that takes a derivative-free approach to optimization by constructing a simplicial homology over the surface of the objective function mentioned in equation (3). The SH attempts to find the quasi-equilibrium solutions to extract all the local minima from the sampled search space. While the algorithm guarantees convergence, it also provides the capability of tracking the progress of optimization and providing non-linear constraints to the objective function (Endres et al., 2018).
The SH algorithm (depicted briefly in Fig. 2) can be divded into four basic phases: In the first phase, we generate the sample points in the constrained search space. These are the vertices of the simplicial complex to be constructed. In the second phase, the directed graph for the simplicial homology is constructed by triangulating these points. In the third phase, minimizer pools are generated from the mesh of simplicial complexes generated in the previous stage. This makes use of Sperner’s lemma (Bapat, 1989) in a repeated manner to generate the pools for local minimization on top of the points sampled in the first phase from the objective function. In the last stage, local minimization is performed in each of the pools to thus reach a stationary global minimum. To perform the local minimization in each pool, we used the Sequential Least SQuares Programming (SLSQP) (Jones et al., 2001) (Golub and Saunders, 1969) optimization method.
One key aspect for setting up the SH optimizer is the sampling scheme to construct the homology groups from the objective function’s hypersurface. For this purpose, we make use of the Sobol sequential sampling method (Sobol, 1967) within the SH global optimizer. We notice that the SH optimizer also gives a significant speedup as compared to the DE optimizer.
2.4 Differential Evolution Optimization
In the MIX approach to optimization, elitism based Genetic Algorithm (GA) was used as the optimizer for the nonlinear parameters. While GA are a very good approach to fit exponential functions, the stochastic nature of parameter tuning gives rise to a certain level of randomness in the fitting, which in turn affects reproducibility. We extend the previous work on MIX (Farooq et al., 2016) with an improvement in fitting using simplicial homology (SH) and differential evolution (DE) algorithms (Price et al., 2014).
The IVIM model consists of exponential functions where parameters are hard to find (Istratov and Vyvenko, 1999). To mitigate this difficulty, we make use of the latin hypercube initialization for the DE step. This improves sampling over the space, compared to random sampling. We make use of the ’best1bin’ strategy within the DE optimizer, as defined in (Price et al., 2014). At each step of the DE, we make use of the Limited memory Broyden-Fletcher-–Goldfarb-–Shanno (L-BFGS-B) optimization (Byrd et al., 1996) to improve the minimization process of the algorithm. The differential evolution this gives us an initial estimation of the model parameters: and .
2.5 Convex Optimization
Based on these initial optimization steps, we have good initial estimates of the non-linear parameters . Using these initial estimates obtained from the DE, finding the is a convex problem and can be solved via a convex optimizer (Boyd and Vandenberghe, 2011) (Diamond and Boyd, 2016) via constrained linear-least squares. The constraint imposed on the volume fractions is and is used in the objective function of the convex optimizer as follows:
2.6 Nonlinear Least Squares Fitting
As a final step, we make use of the initial estimates of the volume fractions and the nonlinear components as initialization for the full IVIM model. This ensures that the inter-dependency among the linear and nonlinear parameters is taken into account, with the final fit resulting in a stationary point. We make use of the Trust Region Reflective Algorithm (Branch et al., 1999) to perform this type of optimization, due to its robustness and capability to handle sparsity in the data. It makes use of the reflective transformation and large scale approximation to find the optima of the search space (Jones et al., 2001).
3 Experiments and Results
The effort of this work has been to provide a good simultaneous fitting method via newer and more advanced optimization approaches. For the first time, we provide the Simplicial Homology (SH) optimization approach to improve the speed of global optimization process (2X speed-up shown in Fig. 6) initially proposed in MIX (Farooq et al., 2016). To demonstrate that the proposed framework improves the quality of estimates amd to ensure the stability and robustness of fitting, we conducted experiments with both simulated (with and without noise) and real data (Peterson, 2016) We simulated the signal from the IVIM model in DIPY (Garyfallidis et al., 2014) and fitted it with the proposed framework. Fig. 3 depicts different parameters used for simulating the data and their respective estimations. We can clearly see that the model gives almost exact estimates for the signal for each case. Additionally, signals have been simulated and fitted with both Rician and Gaussian noise at varying SNR level (Fig. 3), demonstrating the high quality of estimation and its robustness to noise.
Furthermore, we evaluate the goodness-of-fit per voxel using a 2-fold cross validation procedure (Rokem et al., 2015) and quantified via a (coefficient of determination) statistic. To do so, we fit the IVIM model to a part of the data (learning set) and then use the model to predict a held-out set (test set). Prediction of the held out data is done and recorded iteratively over the whole dataset. At the end of these iterations, a prediction of all of the data is compared directly to all of the data in that voxel. The score of the proposed method is always higher compared against other fitting methods.
Additionally, in Fig. 6 we show that the results obtained from the VarPro estimation have a lower MSE in predicting the as compared to the pervasive multi-stage NLLS method with and without fixing the parameter. It has also been shown in the literature that for the liver and pancreas (Gurney-Champion et al., 2018), setting the perfusion to improves the quality of fitting. We compare our results against the same and show that MSE obtained from the model fitting for the predicted S0 is lower, with the advantage of not having to set any empirical constraints.
As demonstrated in Fig. 5, when we fix the pseudo-diffusion/ perfusion parameter, the maps may end up with an artifact, as the estimation cannot find values for regions that go out of bounds (see comparison in bottom row). Furthermore, the figure also demonstrates that our method improves the estimation of the perfusion and diffusion fractions with more informative contrast maps. This is primarily due to the correct estimation of the associated fractions of perfusion and diffusion. VarPro fitting finds a good segregation between the parameter subspaces. In Fig. 4, we also show the goodness-of-fit at lower b-values, where the IVIM perfusion effect is more apparent.
In conclusion, we show that the proposed fitting method for parameter estimation is a promising method for IVIM. From Fig. 5, we see that the perfusion fraction is more interpretable, as compared to that obtained with other methods. We also see that the perfusion coefficient () has an agreement with the perfusion fraction obtained from fixing the to a certain value. While we can see that the MSE of the predicted is lower than that of the other methods (see Fig. 6), the diffusion coefficient () is also properly estimated without artifacts from fixing other parameters. Therefore, we quantitatively and qualitatively show that the estimation is more precise than other methods.
We provide new and improved strategies for estimating the IVIM parameters. For the first time, we make use of the Simplicial Homology Global Optimizer (SH) (Endres et al., 2018) to speed up the fitting process. We provide means of model selection via k-fold cross validation and show that the results obtained from this fitting are better via the statistic. Our results show that the estimation maps obtained from the SH method of global optimization give the same results as DE and we therefore recommend using SH to gain speed performance. Our results show that the quality of estimation is better than other contemporary methods, and also overcomes the problem of bias that occurs in Bayesian (While, 2017) approaches through simultaneous fitting. While the proposed method has been applied and tested for IVIM, it can be easily extended to other bi-exponetial models in Diffusion MRI such as Free Water Diffusion Tensor Imaging (FW-DTI). Future work will extend this framework in DIPY to provide alternatives to fitting microstructural models via more advanced optimization and machine learning methods (Fadnavis et al., 2019).
This paper introduces a novel method for IVIM model fitting via multistage optimization for automatic parameter estimation. The key features of this method include a combination of simplicial homology based topological global optimization and variable projection. Results show that the quality of estimation is improved with an increased accuracy in the fitting.
NIH CRCNS grant R01EB027585 supported the work of Shreyas Fadnavis, Eleftherios Garyfallidis and Ariel Rokem. We would also like to thank the DIPY community (https://dipy.org) for distributing the method and providing meaningful feedback.
- Bapat (1989) Bapat, R. B. (1989). A constructive proof of a permutation-based generalization of sperners lemma. Mathematical Programming, 44(1-3):113–120.
- Bell (1965) Bell, E. L. (1965). Fitting multi-component exponential decay curves by digital computer.
- Bihan (2019) Bihan, D. L. (2019). What can we see with ivim mri? NeuroImage, 187:56–67.
- Bihan et al. (1988) Bihan, D. L., Breton, E., Lallemand, D., Aubin, M. L., Vignaud, J., and Laval-Jeantet, M. (1988). Separation of diffusion and perfusion in intravoxel incoherent motion mr imaging. Radiology, 168(2):497–505.
- Bihan et al. (2017) Bihan, D. L., Ichikawa, S., and Motosugi, U. (2017). Diffusion and intravoxel incoherent motion mr imaging–based virtual elastography: A hypothesis-generating study in the liver. Radiology, 285(2):609–619.
- Boyd and Vandenberghe (2011) Boyd, S. P. and Vandenberghe, L. (2011). Convex optimization. Cambridge Univ. Pr.
- Branch et al. (1999) Branch, M. A., Coleman, T. F., and Li, Y. (1999). A subspace, interior, and conjugate gradient method for large-scale bound-constrained minimization problems. SIAM Journal on Scientific Computing, 21(1):1–23.
- Byrd et al. (1996) Byrd, R., Peihuang, L., and Nocedal, J. (1996). A limited-memory algorithm for bound-constrained optimization.
- Byrd et al. (1988) Byrd, R. H., Schnabel, R. B., and Shultz, G. A. (1988). Approximate solution of the trust region problem by minimization over two-dimensional subspaces. Mathematical Programming, 40-40(1-3):247–263.
- Diamond and Boyd (2016) Diamond, S. and Boyd, S. (2016). CVXPY: A Python-embedded modeling language for convex optimization. Journal of Machine Learning Research, 17(83):1–5.
- Endres et al. (2018) Endres, S. C., Sandrock, C., and Focke, W. W. (2018). A simplicial homology algorithm for lipschitz optimisation. Journal of Global Optimization, 72(2):181–217.
- Fadnavis et al. (2019) Fadnavis, S., Reisert, M., Farooq, H., Afzali, M., Hu, C., Amirbekian, B., and Garyfallidis, E. (2019). Microlearn: Framework for machine learning, reconstruction, optimization and microstructure modeling.
- Farooq et al. (2016) Farooq, H., Xu, J., Nam, J. W., Keefe, D. F., Yacoub, E., Georgiou, T., and Lenglet, C. (2016). Microstructure imaging of crossing (mix) white matter fibers from diffusion mri. Scientific Reports, 6(1).
- Federau et al. (2013) Federau, C., Meuli, R., Obrien, K., Maeder, P., and Hagmann, P. (2013). Perfusion measurement in brain gliomas with intravoxel incoherent motion mri. American Journal of Neuroradiology, 35(2):256–262.
- Ferrante et al. (1987) Ferrante, J., Ottenstein, K. J., and Warren, J. D. (1987). The program dependence graph and its use in optimization. ACM Transactions on Programming Languages and Systems, 9(3):319–349.
- Garyfallidis et al. (2014) Garyfallidis, E., Brett, M., Amirbekian, B., Rokem, A., Walt, S. V. D., Descoteaux, M., and Nimmo-Smith, I. a. (2014). Dipy, a library for the analysis of diffusion mri data. Frontiers in Neuroinformatics, 8.
- Golub and Pereyra (2003) Golub, G. and Pereyra, V. (2003). Separable nonlinear least squares: the variable projection method and its applications. Inverse Problems, 19(2).
- Golub and Saunders (1969) Golub, G. H. and Saunders, M. A. (1969). Linear least squares and quadratic programming.
- Gurney-Champion et al. (2018) Gurney-Champion, O. J., Klaassen, R., Froeling, M., Barbieri, S., Stoker, J., Engelbrecht, M. R. W., Wilmink, J. W., Besselink, M. G., Bel, A., Laarhoven, H. W. M. V., and et al. (2018). Comparison of six fit algorithms for the intra-voxel incoherent motion model of diffusion-weighted magnetic resonance imaging data of pancreatic cancer patients. Plos One, 13(4).
- Halvorson (1992) Halvorson, H. R. (1992).  padé-laplace algorithm for sums of exponentials: Selecting appropriate exponential model and initial estimates for exponential fitting. Methods in Enzymology Numerical Computer Methods, page 54–67.
- Iima and Kataoka (2018) Iima, M. and Kataoka, M. (2018). Ivim mri of the breast. Intravoxel Incoherent Motion (IVIM) MRI, page 173–194.
- Ismail and Halim (2017) Ismail, I. and Halim, A. H. (2017). Comparative study of meta-heuristics optimization algorithm using benchmark function. International Journal of Electrical and Computer Engineering (IJECE), 7(3):1643.
- Istratov and Vyvenko (1999) Istratov, A. A. and Vyvenko, O. F. (1999). Exponential analysis in physical phenomena. Review of Scientific Instruments, 70(2):1233–1257.
- Jalnefjord et al. (2018) Jalnefjord, O., Andersson, M., Montelius, M., Starck, G., Elf, A.-K., Johanson, V., Svensson, J., and Ljungberg, M. (2018). Comparison of methods for estimation of the intravoxel incoherent motion (ivim) diffusion coefficient (d) and perfusion fraction (f). Magnetic Resonance Materials in Physics, Biology and Medicine, 31(6):715–723.
- Jones et al. (2001) Jones, E., Oliphant, T., Peterson, P., et al. (2001). SciPy: Open source scientific tools for Python.
- Kang et al. (2013) Kang, K. M., Lee, J. M., Yoon, J. H., Kiefer, B., Han, J. K., and Choi, B. I. (2013). Intravoxel incoherent motion diffusion-weighted mr imaging for characterization of focal pancreatic lesions. Radiology, page 122712.
- Kim et al. (2013) Kim, H., Suh, C., Kim, N., Choi, C.-G., and Kim, S. (2013). Histogram analysis of intravoxel incoherent motion for differentiating recurrent tumor from treatment effect in patients with glioblastoma: Initial clinical experience. American Journal of Neuroradiology, 35(3):490–497.
- Lai et al. (2017) Lai, V., Lee, V. H. F., Lam, K. O., Huang, B., Chan, Q., and Khong, P. L. (2017). Intravoxel incoherent motion mr imaging in nasopharyngeal carcinoma: comparison and correlation with dynamic contrast enhanced mr imaging. Oncotarget, 8(40).
- Li et al. (2017) Li, Y. T., Cercueil, J.-P., Yuan, J., Chen, W., Loffroy, R., and Wáng, Y. X. J. (2017). Liver intravoxel incoherent motion (ivim) magnetic resonance imaging: a comprehensive review of published data on normal values and applications for fibrosis and tumor evaluation. Quantitative Imaging in Medicine and Surgery, 7(1):59–78.
- Mitra and Bhatia (2014) Mitra, R. and Bhatia, V. (2014). The diffusion-klms algorithm. 2014 International Conference on Information Technology.
- Osborne and Smyth (1995) Osborne, M. R. and Smyth, G. K. (1995). A modified prony algorithm for exponential function fitting. SIAM Journal on Scientific Computing, 16(1):119–138.
- Paulavičius and Žilinskas (2014) Paulavičius, R. and Žilinskas, J. (2014). Simplicial global optimization. SpringerBriefs in Optimization.
- Peterson (2016) Peterson, E. (2016). Ivim dataset.
- Price et al. (2014) Price, K., Storn, R. M., and Lampinen, J. A. (2014). Differential evolution a practical approach to global optimization.
- Rios and Sahinidis (2012) Rios, L. M. and Sahinidis, N. V. (2012). Derivative-free optimization: a review of algorithms and comparison of software implementations. Journal of Global Optimization, 56(3):1247–1293.
- Rokem et al. (2015) Rokem, A., Yeatman, J. D., Pestilli, F., Kay, K. N., Mezer, A., van der Walt, S., and Wandell, B. A. (2015). Evaluating the accuracy of diffusion MRI models in white matter. PLOS ONE, 10(4):e0123272.
- Sobol (1967) Sobol, I. (1967). On the distribution of points in a cube and the approximate evaluation of integrals. USSR Computational Mathematics and Mathematical Physics, 7(4):86–112.
- While (2017) While, P. T. (2017). A comparative simulation study of bayesian fitting approaches to intravoxel incoherent motion modeling in diffusion-weighted mri. Magnetic Resonance in Medicine, 78(6):2373–2387.