Reliable Estimation of the Intra-Voxel Incoherent Motion Parameters of Brain Diffusion Imaging using θ-Teaching-Learning-Based Optimization

Intra-voxel incoherent motion (IVIM) imaging can characterize diffusion and perfusion of tissues. Traditionally, the least-square method has been used to determine IVIM parameters consisting of pure diffusion coefficient (D), pseudo-diffusion coefficient (D) and the microvascular volume fraction (f). This paper proposes an accurate estimation method for IVIM parameters in human brain tissues using θ-teaching-learning-based-optimization (θ-TLBO). θ-TLBO as an evolutionary algorithm provides high quality solutions for parameter estimations in curve fitting problems. Evaluation of the proposed method was performed on simulated data with different levels of noise and experimental data. The estimated parameters were compared with the results of TLBO and three conventional algorithms: Segmented-Unconstrained (“SU”), Segmented-Constrained (“SC”) and “Full”. The results show that the proposed θ-TLBO has higher accuracy, precision and robustness than other methods in estimating parameters of simulated and experimental data in human brain images especially in low SNR images according to analysis of variance (ANOVA), coefficient of variation (CV), relative bias and relative root mean square errors.


Introduction
Diffusion-weighted magnetic resonance imaging (DW-MRI) is a method to measure microscopic motion associated with diffusion of water molecules and a fast diffusion component associated with microcirculation of blood in the capillary network (perfusion). Sensitivity to diffusion which is known as the 'b-value' will be changed by changing the strength and duration of the diffusion gradients in the acquisition process. When the b-value is low, micro-perfusion causes rapid signal decay. Le Bihan et al. [1] proposed the bi-exponential intra-voxel incoherent motion (IVIM) model (1) to separate these two movements which are called diffusion and perfusion. This model is characterized by three parameters: the pure diffusion coefficient (D), pseudo-diffusion coefficient (D * ) related to the incoherent microcirculation and the micro-vascular volume fraction (f).
where S(b) and S 0 are the signal intensities in DW-MRI with diffusion sensitivity of b and 0, respectively, measured in each individual voxel.
The IVIM parameters are calculated for each voxel by fitting the intensity decay of DW-MRI images in different b-values to the bi-exponential IVIM model which is considered as an optimization problem. Some researchers have shown that the fitting algorithm may have impact on the obtained results [2][3][4][5]. In recent years, several studies have been conducted to propose IVIM curve fitting methods. Least-square minimization methods such as Levenberg-Marquardt (LM) algorithm have been widely used for fitting IVIM signal data [4], [6][7][8][9][10]. Barbieri et al. [3] compared variability, precision and accuracy of six common fitting algorithms for IVIM parameters estimation in upper abdominal organs which included LM ("Full") algorithm, Trust-Region (TR) algorithm, Fixed-D algorithm, Segmented-Unconstrained ("SU") algorithm, Segmented-Constrained ("SC") algorithm and Bayesian-probability based algorithm. They reported the variability between the mentioned algorithms in all abdominal regions. Young Cho et al. [11] have compared two common fitting algorithms named "Full" and "SU" algorithms and showed that the second one has higher precision and accuracy in comparison with "Full" in breast cancer. A similar study was devoted by Suo et al. [6] to assess the performance of three curve-fitting methods based on the LM algorithm for breast cancer. They reported that IVIM parameters obtained using different methods are different depending on the calculation methods. They concluded that the "SU" algorithm was more valid compared with the other two algorithms. Recently, Jalnefjord et al. [4] compared four fitting methods for estimating D and f parameters in tumor and healthy liver tissues including segmented IVIM fitting, leastsquares fitting and Bayesian fitting using marginal posterior modes or posterior means. They concluded that the segmented approach gives better results by considering numerical complexity and computational time.
Although the LM algorithm has a good performance for nonlinear least-squares curve-fitting problems, it is sensitive to initialization [12]. In other words, the algorithm may converge to the global minimum only if the initial guess is close to the final solution; otherwise, the algorithm only is able to find the local minimum.
Recently, evolutionary algorithms (EA) inspired by the biological evolution have gained much focus in solving optimization problems. These algorithms often start with generating an initial population within a given search space. Then, the fitness function is computed for this generated population individually to determine the quality of solution. After that, mutation of the population takes place. The best individual has a higher chance for reproduction in order to optimize solutions. This process is iteratively repeated until a termination criterion is satisfied. Since these algorithms generate high quality solutions for difficult problems, they have been widely applied to solve scientific and engineering problems. The EAs have reasonable calculation time and easy implementation, robustness and require little information about the problem being solved [13].
The teaching-learning-based optimization (TLBO) algorithm proposed by Rao et al. [14] is an evolutionary algorithm that does not require any algorithm-specific parameter except the common control parameters. This algorithm is based on the teaching and learning mechanism in a class in which learning is achieved either by the teacher's teaching (known as teacher phase) or by interaction between learners (known as learner phase). Therefore, the procedure of the TLBO has two phases: i) the teacher phase and ii) the learner phase. The θ-TLBO algorithm is a modified version of the TLBO algorithm which is presented by Niknam et al. [15]. In this algorithm, the phase angles are allocated to the design variables contrary to use of design variables themselves in the TLBO optimization process. Since the feasible searching space is more compact than the TLBO algorithm, the searching process is more precise and the convergence time is shorter especially in problems with large searching space [15].
The purpose of this study is to propose a method on the basis of the θ-TLBO algorithm to estimate the IVIM parameters in human brain tissue. The proposed method is evaluated using simulated data with different noise levels and experimental DW-MRI data from selected region of interest (ROI) in the brain. The estimated parameters were compared with those obtained using "SU", "SC" and "Full" algorithms as well as applying TLBO algorithm.

Background
Levenberg-Marquardt algorithm: The Levenberg-Marquardt algorithm [16], [17] is a standard technique for solving nonlinear least-squares problems and has been widely used for parameters estimation in the context of IVIM analysis [3,5,6,18,19]. Curve-fitting with LM algorithm interpolates between the Gauss-Newton algorithm and the method of gradient descent. The algorithm iteratively reduces the sum of the squares of the errors between the bi-exponential IVIM function and the measured data points in each voxel. The minimization is started by providing an initial guess for the parameters. Since the LM algorithm can have many local minima, it is not guaranteed that algorithm converges to global minima in the parameter space if the initial guess is far from it [20].

Segmented-Unconstrained algorithm:
In this algorithm, by assuming small contribution of the pseudo-diffusion term in the high b-values (b > 200 s/mm 2 ), its effect was neglected and diffusion (D) was obtained by fitting the simplified mono-exponential equation: where S int is the value in that S 0 would satisfy if there is not any pseudo-diffusion phenomenon in the underlying tissue.
Then, the f parameter was estimated as follows: by owing D and f parameters, D * was fitted using (1) over all b-values applying the LM algorithm with least-squares minimization. Since the D * parameter is greater than D (at least 10 times), the initial guess for D * in the LM algorithm was considered ten times greater than estimated D [3], [6].

Segmented-Constrained algorithm:
As in the segmented-unconstrained algorithm, the diffusion D was estimated by mono-exponential (2). Then, D * and f were obtained using the LM algorithm. In this algorithm, the initial value of f parameter was determined according to (3) and the initial value of the D * parameter was taken like the "SU" algorithm.

Full algorithm:
In the "Full" algorithm, three IVIM parameters were estimated using the LM algorithm in each voxel simultaneously based on (1). The values obtained from (2) and (3) were considered as the initialization values of D and f parameters, respectively. The initial value for D * was considered ten times greater than D.

θ-Teaching-Learning-Based Optimization (θ-TLBO) algorithm:
The TLBO algorithm is originally proposed by Rao et al. [14] and modified to θ-TLBO algorithm by Niknam et al. [15]. It is a population-based iterative algorithm that can simulate the teaching-learning process in a given classroom. The algorithm searches for an optimized solution without using any function derivative. The most important property of this algorithm is that it does not have any adjusting tuning parameter and it only needs to control the usual parameters such as the population size and the number of iterations [21].
In both TLBO and θ-TLBO algorithms, after generating the population randomly as X = [x 1 ,x 2 ,…,x n ] within feasible range [x min , x max ], the fitness function f(x i ) is computed for each individual x i . The optimization process in the TLBO algorithm is done on the design variables, while in the θ-TLBO algorithm, phase angles are assigned to the real value of the design variables and then the optimization process is performed on these phase angles. In θ-TLBO algorithm, all the design variables are mapped to θ-space and conversely as follows: With the bijective mapping and one-to-one relation- ship between x and  in (4), they can easily map together.
Then, the fitness function f(x i ) is computed for each individual population x i and the phase angle corresponding with the best solution is considered as the teacher (θ teacher ) and other populations are assumed as learners. In order to tend the population towards the global optimum and break the members out of a local optimum, the update phase of TLBO-based algorithms is divided into two steps: 'Teacher Phase' and 'Learner Phase.' This means that the learners may improve their knowledge either by learning directly from the teacher (teacher phase) or by reinforcing learning among themselves through influencing each other (learner phase). In θ-TLBO algorithm these two phases are formulated as follows [15]: Teacher Phase: In this phase, each phase angle corresponding to the i th learners (θ i ) gain knowledge from the phase angle of the teacher (θ teacher ) and from the phase angle of the mean value of the class (θ mean ): where θ old,i , θ teacher and θ mean are mapping of individual x i , x teacher and x mean to phase angle, respectively. They are mapped to (-π/2, π/2). Also, T F and r are random factors which have been placed to make a diverse population. T F is a teaching factor which can be either 1 or 2, randomly and r is a random number in the interval range of 0 and 1.

Learner phase:
In this phase, another phase angle corresponding to the j th learner (θ j ) is selected randomly such that (j ≠ i) to exchange its knowledge with the phase angle corresponding to the i th learner (θ i ). If θ j has more knowledge than θ i , θ i is moved toward θ j ; otherwise, it is moved away from θ j according to (6): After applying each phase, the fitness function is computed for each new individual x new,i related to θ new,i . x old,i will be replaced only if x new,i has a better fitness function than the x old,i . These two phases are iteratively continued until the termination criterion is satisfied.

Proposed Method
In this paper, a curve fitting method based on the θ-TLBO algorithm for estimating IVIM parameters from each voxel is proposed. Since low signal to noise ratio (SNR) is observed in DW-MRI images, a 3 × 3 Gaussian filter was applied on DW-MRI images to smooth them as a preprocessing step. Then f, D * and D parameters are estimated using the θ-TLBO algorithm according to the general constrained minimization problem as written in (7).
where f(X) and F(X) are the fitness function and the objective function of the problem, respectively. The parameters l and m are the numbers of equality g(X) and inequality h(X) constraints of the problem, respectively. The optimal value of the fitness function f(X) subjected to constraints g a (X) and h b (X) is known as the best solution. In our optimization problem, the penalty coefficients λ 1 and λ 2 are considered zero because there is not any constraint in our problem. Hence: In order to solve the IVIM curve fitting optimization problem, the θ-TLBO fitness function is defined as (9):

Data Simulations
In order to assess the estimation errors of the proposed method, a simulated heterogeneous tumor map with 100 × 100 pixels consisting of three regions with known ground truth were generated based on (1), S 0 = 200 and b = 0,5,50,100,200,270,400,600,800 s/mm 2 . Three sets of following reference parameters were used for three regions in the ground truth maps: x 1 = (0.15,0.01,0.001) for inner part, x 2 = (0.25,0.02,0.002) for middle part and x 3 = (0.35,0.03,0.003) for border part of maps [9]. To achieve more realistic simulation data, the noise factor in DW-MRI was modeled. Since both real and imaginary parts of complex MRI images are degraded by white Gaussian noise with the same standard deviations (SD), the nature of the magnitude signal follows Rician distribution [22]:

MR Imaging
DW-MRI was performed on a 1.5-T whole-body MRI scanner (Ingenia CX, Philips Healthcare, Best, The Netherlands) which is installed in Tabesh Medical Imaging Center, Shiraz, Iran, using a 66 mT/m gradient system and a 15-channel head-spine coil. DW-MRI was acquired with SENSE approach by single-shot spin-echo echo-planar imaging sequence and seven b-values, b = 0,50,100,200,400,600,800 s/mm 2 . Axial DW-MRI was obtained with slice thickness of 5 mm, inter-slice gap of 1 mm, repetition time (TR) of 3.5 s, and 88 ms echo time (TE). The EPI factor was set to 75, field of view (FOV) was adjusted to 230 × 230 mm 2 , flip angle was 90, acquisition matrix was 152 × 103 and receiver bandwidth to 1610 Hz/pixel. Total imaging time was 93 s. Seventeen subjects (11 healthy and 6 tumor cases with the mean age of 38.1±10 years and 52.5±10.2 years, respectively) were participated in the present study and their oral consents were obtained. Voxels within the White Matter (WM) regions for healthy subjects and within the tumor tissue for tumor subjects on the axial sections of the images were sampled by a radiologist as the ROI. ROIs on tumor were drawn on most cellular parts of tumors including areas of necrosis and hemorrhage.

Statistical Analysis
A repeated measures analysis of variance (ANOVA) with a Greenhouse-Geisser correction was conducted on all in vivo data to determine whether there is an overall difference between estimated parameter values obtained by different algorithms. Post hoc tests using the Bonferroni correction was applied to pairwise comparison between algo-rithms. The accuracy of the fitting algorithms on simulated data was assessed by relative bias criteria. It is obtained by calculating the difference between the true IVIM parameters and the estimated IVIM parameters from different fitting algorithms, normalized to the true value. The precision of the parameter estimation on simulated data was assessed by relative root mean squared (RRMS) error for each parameter and algorithm separately. It is calculated by dividing the root mean squared error by the true value and represented as a percentage. The variability of the obtained parameters for each patient in considered ROI was computed by coefficient of variation (CV) for each fitting algorithm. The statistical significance of the difference in the variability of the parameter estimates obtained by different fitting algorithms was examined using a two-tailed paired student's t-test (p < 0.05). All algorithms have been implemented in MATLAB (R2013a); also, SPSS statistics (IBM, v.22) has been applied for statistical analysis.

Simulation Results
In order to assess the accuracy and precision of the estimated parameters from different methods, computer simulations were conducted. In this study, eight synthetic noisy images based on a reference noise free heterogeneous tumor sample were generated using (1) and (10). The IVIM parametric maps from simulated data were estimated using the proposed θ-TLBO-based method, TLBO-based and three other algorithms: "SU", "SC" and "Full". Figure 2 illustrates the ground truth maps and the estimated IVIM parametric maps using different algorithms in the simulated image with SNR = 20. The differences between the true IVIM parameters (ground truth) and the estimated values using relative bias and RRMS were computed for each parameter and each considered fitting algorithm (Fig. 3). Figures 3a-c show the accuracy of the estimators by relative bias values for the simulations in different SNRs for D, D * and f parameters, respectively. The results show that the estimated D and f values obtained by the proposed θ-TLBO method were more accurate than other algorithms in different noise levels, especially in low SNRs. All fitting algorithms underestimated the D parameter and overestimated the D * and f parameters. Although in lower SNRs the accuracy of all algorithms was decreased, the proposed θ-TLBO method was more robust to noise than other algorithms in all parameters estimation. Regarding accuracy, the D * parameter obtained by the "Full" algorithm has more significant bias than other algorithms. Figures 3d-f show the RRMS error for estimated parameters from simulated images as a function of SNR. As shown, both θ-TLBO and TLBO-based methods proved to be robust against noise. Although in high SNRs "SU" and "SC" algorithms have better performance than θ-TLBO, however, in lower SNRs the TLBO-based algorithms were more precise than other ones. Quantitatively, the proposed θ-TLBO method reduces the overall RRMS . Rows two to six show parametric maps obtained using the proposed θ-TLBO method, TLBO, "SU", "SC" and "Full" algorithms, respectively. Visually, θ-TLBO method estimates more similar parametric maps to the ground truth maps than other algorithms. All images are presented in the same gray-value range.
of the D parameter about 15%, of the D * parameter about 40% and of the f parameter about 20% as compared to "SU" and "SC" algorithms in SNR = 12.5, respectively.
Also, according to both criteria outlined in Fig. 3, the D parameter has the least variability and the most reliability than D * and f parameters.

In-vivo DW-MRI Experiment Results
In order to demonstrate the accuracy of the proposed θ-TLBO method for IVIM parameter estimation from each voxel of in vivo DW-MRI data, we assessed the estimated parameters from voxels within the WM of healthy subjects and within the brain tumor (Tab. 1). As shown, the estimated D * parameter using all the algorithms has higher SD in comparison with other estimated parameters and D parameter has the least variation. However, the proposed θ-TLBO method provides lower variance in comparison with other algorithms in both healthy and patient subjects.
Repeated measures ANOVA with a Greenhouse-Geisser correction (p < 0.05) was conducted on all experimental data (Tab. 2). Significant differences are in bold. The test concluded that there was a significant main effect of algorithm selection on the estimation of D value (F (1.198,19.171)=31.072, p < .001), D * value (F(2.350, 37.595)=17.697, p < .001) and f value (F(1.091,17.450)= 28.788, p < .001). Post hoc tests using the Bonferroni correction was applied to pairwise comparison between the θ-TLBO method and four other algorithms. It revealed that all IVIM parameters values measured by the θ-TLBO method significantly differed from three conventional algorithms "SU", "SC" and "Full", except between the D * value obtained from "SC" algorithm (p = 0.164) and the f value obtained from "Full" algorithm (p = 0.179). Figure 4 illustrates the variability of each IVIM parameters as a bar plot of the CVs over all subjects for each fitting algorithm. Error bars represent 95% confidence interval for mean. The CV value related to each patient over the considered ROI was calculated separately and then averaged over all healthy and tumor groups.
Statistical significance of the difference between the estimated parameters across all fitting algorithms was calculated by one-way ANOVA test (Post hoc comparisons using the Tukey (honestly significant difference) HSD test, p < .05). The lowest CV was achieved for all parameters by applying the proposed θ-TLBO method and the "Full" algorithm have higher CV compared with other algorithms. The CVs for estimated parameters in the healthy group using the proposed θ-TLBO method were 7%, 46% and 16%, respectively. The CVs for estimated D, D * and f parameters for tumor group were 9%, 35% and 18%, respectively. The CV related to D parameter obtained from the proposed θ-TLBO method was significantly different from four other algorithms in both groups (all p < 0.001).
Also, the CVs of D * parameter obtained from the proposed θ-TLBO method in both groups were significantly lower compared with four other algorithms (p < 0.001) as well as the CVs of f parameter (p < 0.001). No significant differences were found among the CVs of all three parameters obtained with the "SC" and the "SU" algorithms (p = 1.00 for D, p = 0.99 for D * and p = 0.993 for f). Also, the CVs related to all three parameters obtained from TLBO algorithm significantly differed from three algorithms "SU", "SC" and "θ-TLBO" in both groups, while there is no significant differences between the CV related to D and D * values obtained from "Full" and "TLBO" algorithms.

Discussion
In this study, we proposed a θ-TLBO-based method to estimate IVIM parameters. The performance of the proposed method was compared with the TLBO algorithm and three widely used fitting algorithms "SU", "SC" and "Full" in healthy and brain tumor subjects. Furthermore, extensive evaluations based on the simulated DW-MRI data were performed to assess the accuracy and precision of the proposed method. According to the results obtained based on the simulated data, the proposed θ-TLBO-based method and TLBO-based algorithm provide more accurate estimation than three conventional algorithms in estimating D and f parameters regarding relative bias criteria. Although in high SNRs, "SU" and "SC" algorithms provide less RRMS error for simulated images, however, the proposed θ-TLBO method was more robust to noise than other algorithms in all estimated parameters. Also, our experiments showed that parameter estimation using θ-TLBO and TLBO-based algorithms reduce RRMS, especially in low SNRs as compared to the conventional fitting algorithms. According to the estimated parameters based on experimental data, considerable variability was seen between f and D * parameters obtained by two segmented algorithms and TLBO-based algorithms. This may be due to the fact that in the segmented algorithms, the parameter estimation process is performed in several steps and separately, while in TLBObased algorithms, the parameter estimation is performed simultaneously for all three parameters. Moreover, f and D * parameters have a higher standard deviation than D parameter. This could be because of the non-linear nature of the IVIM model that makes some IVIM parameters more reliable than others [5]. Repeated measures ANOVA with a Greenhouse-Geisser correction on experimental data showed a significant difference between estimated IVIM parameters using different algorithms (p < 0.05). The evidence of applying post hoc tests supports the hypothesis that the average of the estimated D values obtained by the proposed θ-TLBO or by TLBO methods are significantly different in comparison to three other algorithms (p < 0.001).
However, no difference was seen among the estimated D values obtained by using the proposed θ-TLBO and TLBO methods. On the other hand, there was no difference between the average of D * value that was estimated by the proposed θ-TLBO and TLBO methods as well as the proposed θ-TLBO method and "SC" algorithms. However, "SU" and "Full" algorithms showed a significant difference with the proposed θ-TLBO method for the average D * value. No significant difference between the estimated f value obtained by the proposed θ-TLBO method and the "Full" algorithm was seen, while the f value obtained by the proposed θ-TLBO method and other fitting algorithms differed significantly. The variability of IVIM parameters obtained from the proposed θ-TLBO method was compared with other algorithms in vivo data (healthy and tumor subjects) by CV criteria as the quantitative precision measure. For all three IVIM parameters, the proposed θ-TLBO method demonstrates higher precision in comparison with four other algorithms. Also, our results support findings of previous studies [6], [9] which showed the D * parameter has the least precision, while estimated parameter D has the most precision in all algorithms.
On the other hand, using θ-TLBO instead of the TLBO algorithm in the IVIM problem has some advantages. First, since each of the IVIM parameters has different limitation range; feasible searching space will be large. By using θ-TLBO and mapping individuals from design variables space to the θ-space, phase angles of the individuals are calculated and the ranges of all variables will be the same and map to (-π/2, π/2). Second, since the searching space in the θ-space is much more compact than TLBO, so the local optimum is very close to the global one and the chance of entrapping into local optima is much less than the TLBO. Therefore the θ-TLBO algorithm is more accurate and faster than the TLBO one. Third, searching procedure in the θ-TLBO algorithm is more accelerated than TLBO one because in θ-space by sweeping one solution (angle), a sector is searched while in TLBO point-topoint searching process is done which needs more time and movements to search the feasible space [15].
Although extensive evaluations were performed to assess the accuracy of the proposed method based on the synthetic data, evaluation of experimental DW-MRI was restricted by the number of patients. Since the focus of this paper is developing a new method based on EA for reliable estimation of diffusion parameters, we showed improvements in estimated parameters in comparison with estimated parameters using conventional algorithms based on synthetic data. We are continuously working on augmenting the number of subjects in our dataset.

Conclusion
In this paper, an evolutionary-based method for estimation of the IVIM parameters in the human brain was proposed. Synthetic data with different SNR levels and in vivo DW-MRI data were used to assess the performance of the proposed methods. Based on the simulated and in vivo data, it was shown that the diffusion and perfusion IVIM parameters obtained by using the proposed θ-TLBO method were more accurate and precise than those obtained by using the conventional methods such as "SU", "SC" and "Full" algorithms.