1. Introduction

Breast cancer, which is common among women, is an international concern. There are about 200 000 cases a year, and it is estimated that there are more than 1 million women with undetected breast cancer. Among various techniques to detect breast cancer, infrared (IR) imaging has been widely used since the late 1950’s. In general, the temperature of the human body on the surface of the skin depends on the metabolic activity, the blood flow, and the temperature of the surroundings. Any abnormality in the tissue, such as the presence of a tumor, alters the normal temperature on the skin surface due to increased metabolic activity of the tumor. Therefore, abnormal skin temperature profiles are an indication of diseases such as mastitis, benign tumors, fibrocystic breast disease, and cancer. Lawson [1] was the first to propose the use of thermographic detection of breast cancer, when he observed that the local temperatures of the skin over a tumor were significantly higher (about 2-3 degrees) than normal skin temperatures. Lawson and Chughtai [2] established that the regional temperature difference over an embedded tumor was due to convection effects associated with increased blood perfusion, and increased metabolism around the tumor.

Thermography, also known as thermal or infrared imaging, is a procedure to determine if an abnormality is present in the breast tissue surface temperature distribution. This abnormality in temperature distribution might indicate the presence of an embedded tumor. Amalu [3] and Keyserlingk et al. [4] have reviewed the role of infrared thermography with respect to the role of IR imaging in the detection of breast cancer and risk indicator. They [3, 4] concluded that the IR imaging has a significant role to play in breast cancer screening due to its ability to image metabolic aspects of the breast and along with sensitivity, specificity, and prognostics. Further aspects of IR imaging techniques of breast cancer and detection methods from infrared images are described in detail and compiled in a recent book by Diakides and Bronzino [5].

In order to complement the IR imaging techniques of breast tumor detection, models have been proposed to estimate the temperature distributions over breast with and without tumors. For example, Sudarshan et al. [6] studied the temperature distributions in breast cancer by considering a two-dimensional model first proposed by Romrel and Bland [7]. In this model, Pennes bioheat equation was used to describe the flow of thermal energy. The breast model consisted of a subcutaneous fat layer, followed by a gland layer and a deep muscle layer adjacent to the thoracic wall. Using this model, the effects of the tumor size, depth, and the blood perfusion rate on the surface temperature magnitude and profile were determined. Evolutionary algorithms based on neural networks and genetic algorithms have been proposed to detect abnormalities in several applications (such as ECG, tomography, ultrasound); see, for example, Diakides and Bronzino [5]. Recently, Mital and Scott [8] developed a method to determine parameters of an embedded heat source representing a tumor using infrared imaging. It appears that based on literature there are no methods available to estimate the tumor parameters based on the temperature distribution obtained from IR imaging techniques.

The objective of this study is the development of a methodology to estimate the depth, size, and heat generation rate of an embedded tumor in the breast, based on temperature profile on the surface that may be obtained by infrared thermography. This study differs from other works related to breast cancer in that inverse methods are used in the estimation methodology. The methodology is demonstrated with a simple model in order to illustrate the details of the procedures involved. Several cases of various embedded tumor conditions were presented to illustrate the generality of the method.

2. Estimation Methodology and Procedure

The estimation methodology involves evolutionary algorithms using artificial neural network (ANN) and genetic algorithm (GA) as shown in Figure 1. The ANN is used to map the relationship of tumor depth, tumor size, and the heat generation to the temperature profile over the idealized breast model. The relationship obtained from ANN is compared and validated to that obtained by standard numerical model-based finite element analysis software (ANSYS). GA is used to estimate the tumor parameters (depth, size, and heat generation) by minimizing a fitness function involving the temperature profiles obtained from simulated or clinical data to those obtained by ANN. The fitness function is given as a sum of the squares’ function that relates the data from the observed thermographic temperature profile to the temperature profile for a given set of estimated parameters. It is defined as

Figure 1: Overall estimation methodology for prediction of breast tumor parameters.

where Y and T are vectors containing the observed and estimated temperatures, respectively, and superscript T denotes the transpose of the vector. The observed temperature vector Y contains surface temperatures that would be obtained experimentally from thermal imaging. The estimated temperatures contained in vector T are obtained by the outputs of the trained neural network model using estimates for the unknown parameters (from the population pool of the GA). The vector containing the estimates for unknown parameters is denoted by . In this study, the unknown parameter vector contains the depth, diameter, and the heat generation rate of the source. Once the objective function was formulated, the estimation scheme using GA was used to determine the unknown parameters.

3. Artificial Neural Network Model to Map Breast Temperature Profiles

A multilayer feed-forward ANN with backpropagation learning algorithm is developed to represent the surface temperature profile of the breast as a function of depth, diameter, and heat generation rate of the tumor. The network is developed with the MATLAB Neural Network Toolbox using the “newff” function. The learning rate was set reasonably low at 0.05 to ensure convergence of the algorithm. There are many transfer functions available in MATLAB software. After some experimentation, the “logsig” transfer function was chosen for the hidden layers and the output neurons, due to nature of the desired outputs. A schematic of the network is shown in Figure 2. The network has 3 input nodes for the three parameters, and 31 output nodes that define one half of the symmetrical temperature profile on the skin surface. There are 2 hidden layers in the network, the first with 21 nodes, and the second with 19 nodes. This configuration was reached after a few iterations of a single hidden layer network proved to give less than adequate results. Two layer deep-nested loops were used in an exhaustive search that varied the number of neurons in the first and second hidden layers, and found the appropriate combination (21-19) that best trained the network. Therefore, a 3-21-19-31 neural network architecture was developed and trained and tested to validate the model. More details about the neural networks and their concepts can be found in Russell and Norvig [9].

Figure 2: Schematic of the neural network.

In this study, a simplified 2D model of the cross-section of the breast similar to that used by Sudershan et al. [6] was considered, as shown in Figure 3. The Pennes bioheat equation has to be solved to estimate the temperature distribution over the breast model. Due to the complexity of the breast geometry, there are no known analytical solutions to the bioheat equation. Therefore, solutions are usually obtained using computational techniques such as the finite element method. A finite element model created using a commercially available finite element package, ANSYS, solves the Pennes bioheat equation to find the surface temperatures, and this information is used in training/validating ANN model. In this study, the breast is assumed to be of 0.072 m radius hemispherical shape, with a subcutaneous fat layer of 0.005 m uniform thickness, followed by gland, muscle, and the thoracic wall. The tumor is situated symmetrically in the gland layer and assumed to be circular in shape for simplicity. Note that the entire breast is modeled even though the 2D geometry is symmetrical about the Y axis. This is done for future studies where an offcenter tumor location may be considered. The heat transfer in each layer is governed by (1). The values of thermal conductivity, metabolic heat generation rate, and the perfusion source term for each layer, adopted from Werner and Buse [10], are presented in Table 1. The thoracic wall is assumed to be maintained at a constant core body temperature of 37^°C. Based on estimate by Osman and Afify [11], the surface of the breast is subjected to a combined heat transfer coefficient (convection, radiation, and evaporation) of 13.5 -K, a reasonable value in the absence of forced convection. The clinical environment under which the thermographic measurements are made is assumed to be at a constant temperature of 21^°C.

Table 1: Values of thermal conductivity, metabolic heat generation, and perfusion term.

Figure 3: A 2D model of the breast based on Sudershan et al. [6].

Using the finite element analysis of the breast model, a total of 17 different cases were analyzed and used for training the ANN. The inputs parameters are scaled between 0 and 1 using upper and lower bounds for each parameter. The diameter of the tumor is allowed to vary between 0.005 and 0.015 m so that it is in mammary gland region. The bounds for depth are 0.014 and 0.041 m (based on the assumption that the tumor will most likely be present in the gland). The lower and upper bounds on the heat generation rate of the tumor are 20 times (14,000 W/m³) and 100 times (70,000 W/m³) that of normal tissue. The output surface temperatures of the network are scaled values between 0 and 1, using a lower limit of 26^°C and upper limit of 33^°C. The network was trained for 400 epochs or training cycles using the data generated from finite element analysis results. Three different cases were used for validating the trained network.

4. Genetic Algorithms in Breast Tumor Parameters Estimation

In the implementation of the GA, the potential solution (chromosome) contains the depth, diameter, and the tumor heat generation rate as the 3 genes and 10 bits are used to represent the value of each gene in binary format. Therefore, each chromosome has 30 bits. The population size has been chosen to be 20. The initial population pool was generated using a random number generator.

The parameters contained in each chromosome are changed from their binary representation to floating point representation using (for the ith parameter) where P_i is the decimal value of ith parameter, a and b are the lower and upper limits of the search interval, respectively, nb is the number of bits used to represent the parameter (e.g., 10 in this case), and m is the decimal value of the parameter in binary form. The lower limit a and upper limit b for each of the three parameters were mentioned earlier.

The fitness of each chromosome is evaluated using the sum of the squares objective function given by (2). The fitness f of a chromosome is related to the sum of the squares function using where was arbitrarily chosen to be 10 so it always stays greater than the maximum value of sum of squares possible, and S is the sum of the squares generated using a parameter set (chromosome) .

The fitness is determined for each parameter set (chromosome) and then ranked. Chromosomes of higher fitness are more likely to be selected to reproduce and contribute their genetic material to the next generation. A probability of mutation (random changing of a chromosome) of 0.10 and a probability of crossover (formation of children) equal to 0.90 are used to determine how members of the population will reproduce to bring forth the next generation. A flowchart of the solution method is presented in Figure 4. More details about GA procedures and applications can be found in Goldberg [12] and Davis [13].

Figure 4: Flowchart of the GA estimation procedure.

Several cases were simulated in order to validate the ability of the estimation procedure to accurately predict the depth, size, and heat generation rate of the embedded tumor. For each parameter set, two runs of the GA were conducted by introducing 5% and 10% noises in the simulation data to see the sensitivity of the prediction of breast tumor parameters.

5. Results and Discussion

A convergence plot of the mean squared error of the neural network output during training is shown in Figure 5. In general, the results obtained from ANN training were in good agreement with those obtained from finite element analysis model. After the ANN is trained, the network is validated for three different cases as shown in Figure 6. It can be seen from Figure 6 that, for each validation case, the output of the neural network model shows good agreement with the simulation results with values of 0.936, 0.982, and 0.988, respectively. To further validate the ANN results, the temperature profiles obtained for a specific case 3 (depth of 0.029 m, diameter of 0.008 m, and heat generation of 30000 W/m³) are compared with the finite element and ANN models as shown in Figure 7. Once again a good comparison is seen.

Figure 5: Convergence plot with mean squared error of the network output during training.

Figure 6: Comparison of output of trained ANN with finite element simulation.

Figure 7: Comparison of temperature profile of trained ANN with direct numerical simulation.

The parameter estimates for the depth, diameter, and heat generation rate of the tumor for the four cases, corresponding to 5% and 10% random noise in the temperature rise on the skin surface, are presented in Tables 2 and 3, respectively. MATLAB built in function “rand” was used for generating random noise in the temperature rise on the skin surface. For each case, the GA was run twice with different random initial populations yielding similar results. With 5% random noise, the results show good agreement between actual and predicted parameters. In this case, the absolute errors in depth and size were within 0.001 m, and absolute errors in heat generation rate were within 3000 W/m³. With 10% noise, the accuracy of the estimated parameters deteriorates for deep seated tumors with low heat generation rate. In this case, the absolute errors in depth were within 0.005 m, error in size within 0.002 m, and in heat generation rate within 4870 W/m³.

Table 2: Results from parameter estimation procedure to determine embedded tumor parameters using temperature data with 5% noise.

Table 3: Results from parameter estimation procedure to determine embedded tumor parameters using temperature data with 10% noise.

6. Conclusions

An estimation methodology-based evolutionary algorithm using neural networks and genetic algorithms was developed to estimate the breast tumor parameters based on surface temperature profile that may be obtained by infrared thermography. The methodology is demonstrated with a simple model in order to illustrate the details of the procedures involved. For the cases studied, results from simulations show that it is possible to determine the depth, diameter, and heat generation rate from the surface temperature data (with 5% random noise) with good accuracy for the simplified 2D model. With 10% noise, the accuracy of the estimated parameters deteriorates for deep-seated tumors with low heat generation rate. With the recent advances in the accuracy of IR images and noise reduction techniques, it may be possible to reduce noise to within acceptable levels (5%) for the methodology to be effective. In general, the breast shape is complex, and this should be considered for realistic studies. However, the estimation methodology is general and can be applied even to realistic breast geometry. This may increase the computation time but the methodology can be easily extended. In order to further develop the proposed methodology so that it can be used in a clinical scenario, several aspects such as complexities of 3D breast geometry and the effects of nonuniform cooling should be considered in future investigations.