Fitting intravoxel incoherent motion model to diffusion MR signals of the human breast tissue using particle swarm optimization

Article history: Received: 25 July 2018 Accepted: 25 January 2019 Available Online: 15 February 2019 Intravoxel incoherent motion (IVIM) modeling offers the parameters f, D and D from diffusion MR signals as biomarkers for different lesion types and cancer stages. Challenges in fitting the model to the signals using the available optimization algorithms motivate new studies for improved parameter estimations. In this study, one thousand value sets of f, D, D for human breast tissue are assembled and used to generate five thousand diffusion MR signals considering noise-free and noisy situations exhibiting signal-to-noise ratios (SNR) of 20, 40, 80 and 160. The estimates of f, D, D are obtained using Levenberg-Marquardt (LM), trust-region (TR) and particle swarm (PS) algorithms. On average, the algorithms provide the highest fitting performance for the noise-free signals (Radj=1.000) and great fitting performances for the noisy signals with SNR>20 (Radj>0.988). TR algorithm performs slightly better for SNR=20 (Radj=0.947). TR and PS algorithms achieve the highest parameter estimation performance for all the parameters, while LM algorithm reveals the highest performance for f and D only on the noise-free signals (r=1.00). For the noisy signals, performances increase with an increase in SNR. All algorithms accomplish poor performances for D (r=0.01-0.20) while TR and PS algorithms perform same for f (r=0.48-0.97) and D (r=0.85-0.99) but remarkably better than LM algorithm for f (r=0.08-0.97) and D (r=0.53-0.99). Overall, TR and PS algorithms demonstrate better but indistinguishable performances. Without requiring any user-given initial value, PS algorithm may facilitate improved estimation of IVIM parameters of the human breast tissue. Further studies are needed to determine its benefit in clinical practice.


Introduction
Diffusion-weighted MR imaging traces random displacements of water within living tissues using the exponential decay of the diffusion signal amplitude with respect to the degree of field gradient encoding exposed to the tissues during imaging [1].This imaging technique uses no ionizing radiation and requires no contrast agent administration to the patient for quantitative tissue characterization and therefore has become a very popular medical imaging technique in diagnosis and treatment of cancer [2,3].To assess the volume fraction of incoherently flowing blood in the tissue, the diffusion coefficient of water in the tissue and the sum of the pseudo-diffusion coefficient associated to the motion effect and the diffusion coefficient of water in blood, the diffusion MR signal of the tissue needs be processed using intravoxel incoherent motion (IVIM) model [4].This process conventionally involves nonlinear least squares fitting of a bi-exponential decay function to the diffusion MR signal [5].Among the optimization strategies proposed for the fittings, the Levenberg-Marquardt and the trustregion optimization algorithms stand forward due to their easy implementation.The Levenberg-Marquardt algorithm [6] uses a search direction that is a cross between the Gauss-Newton direction and the steepest descent direction that reveals increased robustness.However, the algorithm does not consider any boundary constraints and therefore an estimate may be out of physiologically acceptable range.The trust-region algorithm [7] inherited from the Levenberg-Marquardt algorithm employs a search space restricted to a subset of the domain of a cost function.The algorithm offers estimates within acceptable ranges by incorporating boundary constraints that can be determined easily from the recent research on breast diffusion MR imaging.Both algorithms suffer from the same concern: the initial value of any model parameter to be estimated should be pre-selected correctly to minimize the possible disruptive effect of local minima during fitting.The selection process is quite complicated and requires decent experience on excessive trials with different values.There is a need for an optimization algorithm that requires no user given initial values and that utilizes boundary constraints.Such an algorithm holds a priceless potential to improve the use of IVIM modeling in distinguishing different lesion types and cancer stages of the human breast tissue.The particle swarm optimization algorithm [8] performs search behavior of a swarm of particles hovering through a multidimensional search space.The algorithm eliminates the need for user given initial values while making use of boundary constraints.It has been suggested to optimize several tasks in medicine, such as medical image segmentation, image enhancement and image registration [9][10][11].However, to the best of our knowledge, its use in fitting diffusion MR signals to the IVIM model has not yet been explored.In the current study, a particle swarm optimization algorithm is proposed and compared with the Levenberg-Marquardt and the trust-region optimization algorithms in least squares fitting of the IVIM model to breast diffusion MR signals.

Synthetic diffusion MR signal generation using intravoxel incoherent motion model
The intravoxel incoherent motion (IVIM) model expresses the attenuation in the diffusion MR signal strength acquired for a specific diffusion weighting determined by a b-value, s(b), with respect to the signal strength captured without any diffusion weighting, s(0).The model is formulated using a bi-exponential decay function giving a normalized diffusion signal [12]: Here f, D and D * are the three free model parameters.

Nonlinear least squares fitting
The diffusion MR signals obtained are fitted to the IVIM model using the bi-exponential decay function in a nonlinear least squares fashion.The process involves minimizing the difference between the signal fed to fitting and the signal predicted during fitting: Here, β = [f, D, D * ]. i s and ˆ() i s  denote the signal strength given and the signal strength predicted using a set of values for β for the i-th b-value respectively.n is the total number of b-values considered for the diffusion signal.In the current study, the nonlinear least squares fittings are performed by using three different optimization algorithms, namely Levenberg-Marquardt, trust-region and particle swarm.

Levenberg-Marquardt algorithm
The Levenberg-Marquardt (LM) algorithm [6] uses a search direction that is a cross between the Gauss-Newton direction and the steepest descent direction and therefore offers increased robustness: Here Jk and λk are the Jacobian matrix of derivatives of the residuals with respect to β and the Marquardt parameter both for the k-th iteration.dk is a direction of descent satisfying βk+1= βk+dk.λ is updated from iteration to iteration; increasing the value has the effect of changing both the direction and the length of the shift vector.In the current implementation, λ is initially set to 0.01.After an iteration, if the difference takes a lower value, λ is decreased by a factor of 10; otherwise, it is increased by a factor of 10 for the next iteration.The iteration is stopped when the gradient of the difference reaches its minimum of 10 -12 or the change in model parameters for finite difference gradients reaches its minimum of 10 -9 .No bound constraints are considered during iterations; however, the optimization is started with 0.10, 1.38µm 2 /ms and 110µm 2 /ms as the initial values for f, D and D * respectively.

Trust-region algorithm
The trust-region (TR) algorithm [7] is inherited from LM algorithm and potentially offers better and faster solutions especially when the fitting process is far from the correct solution.Let ∆k be the two-norm of the solution to Eq. 3 updated by each iteration according to standard rules, then dk may be the solution as follows: Fitting intravoxel incoherent motion model to diffusion MR signals of the human breast tissue In the current implementation, the TR algorithm considers the same stopping criteria and the same initial values set for the LM algorithm.Besides, the algorithm also makes use of lower and upper limits for the parameters set as 0≤ f≤ 0.25, 0.25≤ D≤ 3.50 (×µm 2 /ms) and 25≤ D * ≤ 250 (×µm 2 /ms).

Particle swarm algorithm
The particle swarm (PS) algorithm [8] performs search behavior of a swarm of particles hovering through a multidimensional search space to find a solution.
During this search, a particle iteratively adjusts its velocity (v) and position (x) according to Here, w is an inertia weight that balances the global search and local search while c1 and c2 are two scaling factors and r1 and r2 are the randomly generated numbers uniformly distributed between 0 and 1. x b and x ba denote respectively the best previous position of the particle and the best position among all the particles that gives the minimum difference.μ presents the flying time for the particle.In the current implementation, the number of particles in the swarm is set to 100 while μ= 1, w= 1.1 and c1=c2= 1.49.The lower and the upper limits for the parameters are set to the values considered for the TR algorithm.The iteration is stopped when the gradient of the difference reaches its minimum of 10 -12 .

Performance assessment
Success of each optimization algorithm studied is assessed with respect to the model fitting and the parameter estimation performances.The model fitting performance is measured by goodness-of-fit given by adjusted R-squared [15]: Here, n and m represent respectively the number of diffusion signal strength measurements and the number of parameters within the model fitted (for the current study, n= 10 and m= 3).s denotes the average of the diffusion signal strength for different b-values.R 2 adj ranges from 0 to 1 with a higher value indicating a better fit (i.e.better similarity between the diffusion signal predicted by the estimated parameter values and the diffusion signal given by the ground truth parameter values for a specific noise level).To assess the parameter estimation performance of the optimization algorithms, correlations between the estimates by the algorithms and the true simulated values are measured by Person's correlation coefficient (r).The coefficient value ranges between -1 and +1 and a higher absolute value designates a better estimate.The fitting algorithms are numerically implemented and performance metrics are computed using the existing libraries of MATLAB (v8.2;Natick, MA) on a standard PC (Intel i5-4460 3.20GHz processor, 6GB memory and 64-bit OS).Though, the TR algorithm offers the best performance among all of the algorithms (R 2 adj=0.951-0.999).However, for the signals having SNR≥ 40, very good fitting performances are attained by all the algorithms (R 2 adj≥ 0.954).The LM, TR and PS algorithms require respectively 35ms, 39ms and 353ms on average for fitting the IVIM model to a noise-free signal.For a noisy signal, the average computation time for the LM, TR and PS algorithms increases respectively from 40ms to 71ms, 42ms to 58ms and 356ms to 360ms with the decrease in SNR from 160 to 20.The PS algorithm always requires a considerably longer computation time when compared to the LM and the TR algorithms.

A
Estimates of f, D and D* from all diffusion MR signals obtained by the LM, TR and PS algorithms are given in Table 3.For f and D, the same statistical measures (i.e.Mean±SD, the minimum and the maximum) are provided by all the algorithms.However, for D*, the LM algorithm often provides remarkably different statistical measures when compared to the TR and the PS algorithms that reveal indistinguishable measures.In addition, noise has a certain impact on the estimates by the algorithms; when the level of noise increases (i.e.SNR decreases), higher values for f, lower values for D and remarkably higher values for D* are estimated.Table 4 shows the results of parameter estimation performances of the optimization algorithms (i.e. the correlations between the estimated and the true simulated values for f, D and D*).For the noise-free signals, the strongest positive correlation exists for f and D (r=1.00) regardless of the optimization algorithm.However, for D*, very strong positive correlations are recognizable between the true simulated values and their estimates by the TR and the PS algorithms only (r=1.00).Besides, a very weak positive correlation is obtained for the LM algorithm (r=0.08).For the noisy signals, varying positive correlations are observable for f, D and D*.When compared to the LM algorithm, the TR algorithm always offers better correlations.Besides, the TR and the PS algorithms demonstrate almost the same correlations.These algorithms permit very strong In a general sense, the TR and the PS algorithms provide better parameter estimates than the LM algorithm.However, the PS algorithm shows almost the same estimation performance as the TR algorithm.
In contrast to the TR algorithm, the PS algorithm does not require any user-given initial value, therefore it offers a priceless tool in nonlinear least squares fitting of IVIM model to the diffusion MR signals of the human breast tissue especially for low SNRs.Scatter plots of the estimates by the PS algorithm against the true simulated values are presented in Figure 1.

Conclusion
The intravoxel incoherent motion (IVIM) modeling for the diffusion MR signals of the human breast tissue enables assessment of volume fraction of the incoherently flowing blood in the tissue (f), the diffusion coefficient of water in the tissue (D) and the sum of the pseudo-diffusion coefficient associated to the motion effect and the diffusion coefficient of water in blood (D * ) as the potential biomarkers for different types of breast lesions and different stages of breast cancer.In the current study, the model fitting and the model parameter estimation performances of three optimization algorithms, namely Levenberg-Marquart (LM), trust-region (TR) and particle swarm (PS), in nonlinear least squares based fitting the IVIM model to the diffusion MR signals are investigated.Our results from five thousand breast diffusion MR signals generated synthetically show that all the optimization algorithms achieve very good model fitting performances for both noise-free and noisy diffusion signals.However, computation time spent by all the algorithms for fitting the IVIM model increase with the increase in noise level.When compared to the LM and the TR algorithms, the PS optimization algorithm always requires a considerably longer computation time.On the other hand, our results also demonstrate that the optimization algorithms play a very important role in getting reliable and precise estimates of the three model parameters f, D and D * .On the noise-free signals, the LM, TR and PS algorithms all exhibit the highest performance in estimating f and D. However, for D * , most of the estimates by the LM algorithm are out of the physiologically acceptable range and therefore the algorithm shows a remarkably lower performance then the TR and PS algorithms both of which reveal the highest performance possible.In presence of noise, the TR algorithm performs better than the LM algorithm; the PS algorithm possesses the same estimation performance as the TR algorithm.The TR and the PS algorithms perform remarkably better in estimating f and D, especially in presence of high noise levels.However, their performances for f are not as good as for D. In addition to these, all the three algorithms suffer from very poor performances in estimating D * regardless of the noise level.These findings verify the wide use of the estimates of D in recent research on quantitative diffusion weighted imaging, arising the question that f and D * may also offer valuable information under low noise conditions.Overall, the TR algorithm performs better than the LM algorithm while the PS algorithm shows almost the same performance as the TR algorithm.Revealing a very good performance without requiring any usergiven initial value, the PS algorithm can be preferably used in nonlinear least squares fitting of the IVIM model to the diffusion MR signals of the human breast tissue, especially in presence of high levels of noise.The current study has some limitations.First of all, the results are from synthetic MR signals generated for the IVIM model parameters with certain mean and standard deviation values that satisfy ranges reported from a single-center study performed with a 3.0T MR scanner.Ranges recognized by different scanners at different centers may differ and therefore dissimilar MR signals can be of concern.Second, noisy MR signals are obtained for SNRs of 20, 40, 80 and 160 using Gaussian noise.3T MR scanners with optimized imaging protocols easily provide these high SNRs.However, lower SNRs may be experienced in some cases and for signals with SNR<5, it might be questionable to use of Gaussian noise instead of Rician noise.Another issue is that the mean values of the IVIM model parameters generated are later used as the usergiven initial values for the LM and the TR algorithms.For most of the fittings by these algorithms, the initial values might be very close to the solution of the model and this might minimize the disruptive effect of local minima during fitting while leading to convergence to a precise solution.Consequently, performances of the LM and the TR algorithms might be over assessed.In addition, current implementation of the PS algorithm consumes considerably large amount of time in fitting the IVIM model.The use of faster computers or programming techniques dedicated to PS optimizations might lead to much shorter computation time.On the other hand, adjusted R-squared measure gives almost the same model fitting performances for all the algorithms.The use of an alternative measure such as the Akaike Information Criterion or Bayesian Information Criterion may give remarkably different fitting performances.There are some issues awaiting further exploration and improvement.Noise has undesirable impact on most of the IVIM model parameter estimates limiting their use in the detection of different types of breast lesions and stages of breast cancer.To cope with this, an additional parameter estimation may be performed from a "reference" healthy tissue to get a normalized parameter that may be less corrupted by noise.Another solution may be the use of alternative fitting techniques less sensitive to noise then the least squares fitting.For instance, artificial neural networks can be developed to approximate the IVIM model parameters by approximating the bi-exponential decay function.On the other hand, improved methods based on least squares fitting such as segmented least squares fitting [16] or weighted least squares fitting may also be used to minimize the effect of noise [17].In conclusion, by facilitating the estimation of IVIM model parameters from diffusion MR signals of the human breast, the particle swarm optimization algorithm holds a priceless potential to improve quantitative characterization of human breast tissue.Further studies with real clinical data are needed to determine its benefit in distinguishing types of breast tissue and lesions in the diagnosis of breast cancer.

Table 1 .
total of one thousand different value sets for f, D and D * are synthetically generated to mimic the diffusion characteristics of the human breast tissue regarding to the IVIM model.Statistical summary for each parameter is as shown in Table 1.The value sets are used to simulate diffusion MR signals considering noise-free and noisy imaging conditions with four different levels of SNR (i.e. 20, 40, 80 and 160) resulting in the generation of five thousand synthetic breast diffusion MR signals in total.Statistics for the IVIM model parameters f, D and D * simulated.

Table 2 .
Model fitting performances of the optimization algorithms assessed by R 2 adj

Table 3 .
Estimates of f, D and D * by the optimization algorithms.Units of D and D* are [µm 2 /ms] and f is dimensionless.Results are the same for all the algorithms (*) and as the ones estimated by the TR algorithm (**).

Table 4 .
Parameter estimation performances the optimization algorithms assessed by r Results are the same for all the algorithms (*) and as the ones estimated by the TR algorithm (**).intravoxel incoherent motion model to diffusion MR signals of the human breast tissue 109 correlations for D from the noisy signals with SNR≥ 40 (r≥ 0.95) meanwhile for f from the noisy signals with SNR≥ 80 (r≥ 0.90).However, for D*, they demonstrate very weak correlations (r= 0.02-0.20)as the LM algorithm does (r=0.01-0.14)regardless of the SNR.