1 Introduction
The Intravoxel Incoherent Motion (IVIM) model uses data collected with the pulsed gradient spin echo sequence, often used to measure diffusionweighted MRI contrasts, to derive information about blood microcirculation. In IVIM, information about perfusion of blood is estimated as a pseudodiffusion process from a PGSE experiment, in which several low bvalues 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 perfusiondriven 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 bicompartmental 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 nonorthogonal, which makes it hard to project it along the realaxis by taking an integral and (2) The inherent resolution limit in modeling the exponential decay of multiexponential models (Istratov and Vyvenko, 1999). The present work tackles these challenges by simultaneously fitting the model parameters without an empirical fixed threshold on the bvalues. 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 tradeoff 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 diffusionweighted component and a perfusionweighted component, and their exponents represent the apparent diffusion coefficient (ADC) and a pseudodiffusion 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 multiexponential 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 PadeLaplace 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 multicompartment 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 bvalues 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 stateofthe 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 Methods
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 nonstochastic in nature in comparison with evolutional algorithms. The calculated homology groups perform optimization in a derivativefree 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:
(1) 
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:
(2) 
Thus, from equations (1) and (2) we can write our objective function to fit the IVIM model as follows:
(3) 
Now, as mentioned (Golub and Pereyra, 2003), we will project out the volume fractions with the help of a projector function. This is a MoorePenrose inverse of the nonlinear components of our model defined in equation (2). We can write the projection function as:
(4) 
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).
(5) 
This objective function is then used in the differential global optimization stage (via SH or DE) described in further detail in the sections (2.3 and 2.4 respectively).
2.3 Simplicial Homology Global Optimization
The Simplicial Homology (SH) optimizer is a topographic optimizer that takes a derivativefree 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 quasiequilibrium 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 nonlinear 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 BroydenFletcher–Goldfarb–Shanno (LBFGSB) 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 nonlinear 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 linearleast squares. The constraint imposed on the volume fractions is and is used in the objective function of the convex optimizer as follows:
(6) 
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 interdependency 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 speedup 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 goodnessoffit per voxel using a 2fold 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 heldout 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 multistage NLLS method with and without fixing the parameter. It has also been shown in the literature that for the liver and pancreas (GurneyChampion 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 pseudodiffusion/ 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 goodnessoffit at lower bvalues, 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.
4 Discussion
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 kfold 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 biexponetial models in Diffusion MRI such as Free Water Diffusion Tensor Imaging (FWDTI). 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).5 Conclusion
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.
6 Acknowledgments
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.
References
 Bapat (1989) Bapat, R. B. (1989). A constructive proof of a permutationbased generalization of sperners lemma. Mathematical Programming, 44(13):113–120.
 Bell (1965) Bell, E. L. (1965). Fitting multicomponent 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 LavalJeantet, 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 hypothesisgenerating 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 largescale boundconstrained minimization problems. SIAM Journal on Scientific Computing, 21(1):1–23.
 Byrd et al. (1996) Byrd, R., Peihuang, L., and Nocedal, J. (1996). A limitedmemory algorithm for boundconstrained 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 twodimensional subspaces. Mathematical Programming, 4040(13):247–263.
 Diamond and Boyd (2016) Diamond, S. and Boyd, S. (2016). CVXPY: A Pythonembedded 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 NimmoSmith, 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.
 GurneyChampion et al. (2018) GurneyChampion, 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 intravoxel incoherent motion model of diffusionweighted magnetic resonance imaging data of pancreatic cancer patients. Plos One, 13(4).

Halvorson (1992)
Halvorson, H. R. (1992).
[3] 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 metaheuristics 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 diffusionweighted 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 diffusionklms 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). Derivativefree 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 diffusionweighted mri. Magnetic Resonance in Medicine, 78(6):2373–2387.
Comments
There are no comments yet.