Our study group consisted of 13 patients whose ages ranged from 18 to 57 years (mean 40.2 years). Our sample size was 18 mass lesions that consisted of nine malignant lesions and nine benign entities, which were randomly selected among patients with BIRAD 3 and BIRAD 4 classification by the expert radiologist. The inclusion criteria consisted of there being at least one mass (with BIRAD 3 or BIRAD 4 classification) in the subject’s breast and the availability of patient’s MR images with standard imaging protocol. On the other hand, pregnant patients and allergic patients to gadolinium contrast agent were excluded from our study. We did not have any missing value in our database. Informed consent was obtained from all patients or their relatives. All patients underwent a surgical biopsy. Their basic information is presented in
Table 1. Lesion Characteristics
Lesion Number Gender Pathology BIRAD Classification Age Type 1 Female Benign 4b 33 Fibroadenoma 2 Female Benign 4 35 Fibroadenoma 3 Female Benign 4b 18 Intraductal Papilloma 4 Female Benign 4 42 Intraductal Papilloma 5 Female Benign 4b 37 Invasive Ductal Carcinoma 6 Female Benign 4 42 Intraductal Papilloma 7 Female Benign 4a 57 Fibroadenoma 8 Female Benign 4a 40 Intraductal Papilloma 9 Female Benign 4 36 Fibroadenoma 10 Female Malignant 4b 37 Invasive Ductal Carcinoma 11 Female Malignant 4b 37 Invasive Ductal Carcinoma 12 Female Malignant 4b 43 Invasive Ductal Carcinoma 13 Female Malignant 4b 43 Invasive Ductal Carcinoma 14 Female Malignant 4a 53 Invasive Ductal Carcinoma 15 Female Malignant 3 53 Invasive Ductal Carcinoma 16 Female Malignant 4b 43 Invasive Ductal Carcinoma 17 Female Malignant 4c 36 Invasive Ductal Carcinoma 18 Female Malignant 4b 38 Invasive Ductal Carcinoma 3.2. Data Acquisition
Magnetic resonance imaging was performed at the Pardisnoor imaging center during years 2015 and 2016, using a calibrated 1.5 Tesla ACHIEVA Philips MR scanner in the prone position with a specific breast coil. The dynamic study was performed after injection of 0.1 (mmol/kg) of gadopentetate dimeglumine. The DCE portion of the MRI was based on a standard clinical protocol that consisted of a series of T1 weighted, fat-suppressed, three-dimensional gradient echo acquisitions in the axial plane with the following parameters: repetition time, 4.9 - 5.4 ms; echo time, 2 - 2.5 ms; 10° flip angle, field of view, 240 × 240 (mm
2), voxel size, 1 × 1 × 1 (mm 3), matrix, 320 × 320 and temporal resolution of 80 - 90 seconds ( 14). 3.3. Time Series Extraction
The extracted time series from tumor margin, which can be considered as a nonlinear system contains valuable diagnostic information. To capture this information, conventional methods, such as Fourier transform, were not successful due to the nonlinearity of system dynamics. Chaos analysis is a suitable alternative to describe nonlinear and complex systems. In this paper, we considered breast tumors as a nonlinear complex system, which can be described by a strange attractor in phase space (PS). The time evolution of these observables in the state-space establishes a trajectory (
8). We didn’t have enough information about the dynamics of the system. We had only one TS measurement ( 8). In such case, it is not possible to find the exact PS of the system. Therefore, a pseudo-PS may still be constructed. This pseudo-PS is called the reconstructed phase space (RPS) ( 8). In this research, Time Delay Embedding (TDE) method was used for PS reconstruction.
In order to extract the time series from tumor boundary, we selected the best part of MR images for each tumor. The center of gravity of the tumor was calculated in each MR image and tumor margin was traced utilizing Fuzzy C-Mean (FCM) segmentation method. The FCM is a fast and easy to use segmentation method, which can be suitable to segment ill-defined breast tumors from the background (
15, 16). Then, the distance of each pixel on the tumor margin from the center of gravity was measured and used for extracting the time series for all breast lesions. Using this method, we converted the tumor shape into a time series. Figure 1 shows a typical MR image after pre-processing and filtering by the median filter. Figure 2 shows the segmented tumor and its depicted contour using Fuzzy C-Mean segmentation algorithm and “bwtraceboundary” in Matlab.
Figure 1. Filtered Image Using Median Filter
Figure 2. Segmented Tumor Using Fuzzy C-Mean Method
Figure 3 shows time series extracted from the tumor contour of Figure 2. We extracted chaos features from the time series by calculating the LLE of the extracted time series. Figures 4 and 5 demonstrate the extracted benign and malignant mass lesions, respectively, and their depicted contour using the FCM method. Also, the centers of mass of all tumors were shown in these images.
Figure 3. Time Series Extracted From Tumor Contour
Figure 4. Extracted Benign Mass Lesions and Their Depicted Contour Using Fuzzy C-Mean Segmentation Method
Figure 5. Extracted Malignant Mass Lesions and Their Depicted Contour Using Fuzzy C-Mean Segmentation Method
3.4. Embedding Dimension and Time Delay
We utilized time delay embedding to reconstruct the phase space. A chaotic time series can be embedded into an RPS with an embedding dimension m and a time delay J. In the TDE method, two parameters are necessary to be selected and optimized, the embedding dimension, m, and the time delay J. These parameters should be optimal enough due to their impact on calculating the LEs (
8). The inherent system dimension was considered as d. According to Takens’ theorem, if the system dimension is d, we can establish an RPS that is equivalent to the original PS by embedding with a dimension m (greater than 2d + 1) ( 17).
In order to be able to identify an attractor, the time delay should be selected optimally. Many methods have been proposed for selecting an optimal time delay but it is still not possible to apply a unique method for all types of data (
9). We used the autocorrelation indicator for time delay estimation. 3.5. Calculation of Largest Lyapunov Exponent
To calculate the LLE, we considered a discrete system with a 1-D map x
k + 1 = f (xk), which evolves when it is started at two initial conditions of x 0 and (x 0 + ε 0) ( 10, 11). The parameter ε 0 requires a small value to show the two initial states are very close to each other. The Lyapunov exponent is defined when the two trajectories are diverged by a distance ε n after n iterations, as follows ( Equation 1).
Where λ is the Lyapunov exponent (
A practical numerical technique for calculating the LLE is the method developed by Rosenstein et al. (
9). Due to its robustness to time delay and embedding dimension changes, this well-known method was applied in this research to calculate the LLE of the MRI-based time series ( 18). We changed embedding dimension (m) in our program from 2 to 5 and the best value was m = 2. The autocorrelation method was used for time delay estimation. The implementation of this method is easy because it utilizes a simple measure of exponential divergence ( 9, 18). 3.6. Fractal Analysis
The concept of fractals was proposed by Mandelbrot to describe objects with irregular structures (
13). For quantifying the complexity and self-similarity of the structure of an object, a measure known as the fractal dimension (FD) can be utilized ( 12). A self-similar structure was considered, which consists of a number of self-similar pieces at the reduction factor 1/s ( 12). FD can be defined as follows ( Equation 2):
Then, we have (
The most commonly used method for estimating FD is the box-counting method. We used this method for calculation of fractal dimension and other fractal parameters.
3.7. Feature Extraction
We extracted a couple of chaos and fractal features from breast tumor MR images. In addition to LLE, we calculated the mean value of each time series as a measure of tumor size. Also the variance of time series was calculated. The LLE was calculated as a measure of chaos in the extracted time series. We calculated the variance of Lyapunov exponents (VarLE) to obtain another TS feature and quantify the spread of Lyapunov exponents. Also, we utilized fractal analysis to extract significant fractal features besides the time series descriptors.
For classification of benign and spiculated masses, we used roughness in the boundary of masses as one of the main fractal features. Spiculated masses have rough variation in boundaries whereas the benign masses are round and with smooth variation (
13). Therefore, the variation of FDs in different scales was utilized to extract important information for classification. N2 as shown in Equation 4 describes the amount of changes in FDs in different scales with respect to the maximum number of FDs, which is measured on the smallest scale. This feature, especially in low scales, has the information of spiculation with high resolution ( 13).
Where NB indicates FDs of the boundary of the mass.
Another feature, used in the classification of different masses, was circularity of the mass. We proposed two circularity measures in this research, Circ and NonCirc. The former shows the circularity of tumor and the latter indicates the non-circularity of the mass. These features can be defined as
Equations 5 and 6:
Where CL indicates the length of mass contour and TS is the extracted time series of mass contour. The mean value of TS indicates the average radius of the tumor and the difference between TS maximum and minimum can quantify the degree of deviation from the circular shape.
It should be noted that all the implementations of feature extraction process were performed on a 2.53-GHz laptop computer, with 4 GB random access memory, in the MATLAB 188.8.131.525 software environment (MathWorks, Inc).
3.8. Statistical Analysis
In this cross-sectional study we selected 18 samples using a systematic sampling method. We used Pierson’s correlation analysis to present the significance of each parameter. Due to the small size of our database, it cannot have a normal distribution. We performed bootstrap method based on 5,000 bootstrap samples to solve the problem of small sample size. The obtained results are reported in
Table 2. Bootstrap Analysis Based on 5000 Bootstrap Samples
Pathology LLE N2 Mean Var varLE circ1 noncirc Pathology (Pierson’s correlation) 1 0.745 -0.336 -0.307 -0.352 0.419 0.381 -0.494 P Value 0 0.172 0.215 0.152 0.083 0.119 0.037 Bootstrap Bias 0 0.035 0.002 0.009 -0.034 0.068 -0.002 0.009 Std. Error 0 0.089 0.211 0.203 0.133 0.138 0.201 0.184 95%Confidence Interval Lower 1 0.614 -0.695 -0.639 -0.647 0.249 -0.061 -0.822 Upper 1 0.943 0.118 0.160 -0.105 0.796 0.718 -0.094 3.9. Statistical Software
Statistical analysis was performed using the SPSS software for windows, version 19.0.0 (SPSS Inc., Chicago, IL).
This study was considered by the ethical committee of our university, and as it used available information of patients for regular diagnostic procedures, it was approved by this committee. We did not change any standard diagnostic procedure or imaging protocols.