Novel Meta-Heuristic Model for Discrimination between Iron Deficiency Anemia and B-Thalassemia with CBC Indices Based on Dynamic Harmony Search

In recent decades, attention has been directed at anemia classification for various medical purposes, such as thalassemia screening and predicting iron deficiency anemia (IDA). In this study, a new method has been successfully tested for discrimination between IDA and \b{eta}-thalassemia trait (\b{eta}-TT). The method is based on a Dynamic Harmony Search (DHS). Complete blood count (CBC), a fast and inexpensive laboratory test, is used as the input of the system. Other models, such as a genetic programming method called structured representation on genetic algorithm in non-linear function fitting (STROGANOFF), an artificial neural network (ANN), an adaptive neuro-fuzzy inference system (ANFIS), a support vector machine (SVM), k-nearest neighbor (KNN), and certain traditional methods, are compared with the proposed method.


I. INTRODUCTION
Iron deficiency anemia (IDA) and β-thalassemia trait (β-TT) are the most common types of anemia. They have similar effects on red blood cell indices such as red blood cell count (RBC), hemoglobin concentration (Hb), mean corpuscular volume (MCV), mean corpuscular hemoglobin (MCH), mean corpuscular hemoglobin concentration (MCHC), red blood cell distribution width (RDW), and hematocrit (HCT). Thalassemia disease is a genetic trait that directly impacts on reducing the life span of red blood cells through regulating the formation of hemoglobin (Hb) as a fundamental component of red blood cells [1]. In contrast, IDA is not inherited and can generally be treated by dietary changes and iron supplements. Complete blood count (CBC) is the primary clinical test for diagnosis of anemia. If CBC indices cannot discriminate between IDA and β-TT, complementary tests, such as hemoglobin electrophoresis or genetic tests, are needed. Complementary tests are time consuming, expensive and not available in all laboratories and hospitals. As a result, scientists have conducted many studies on the discrimination between IDA and β-TT anemia classification using the indices of the CBC test [2]. Historically, mathematical formulas have been used to solve the problem [3], [4], [13]- [16], [5]- [12]. Recently, there has been increasing interest in the use of artificially intelligent automated medical diagnostic systems [17]- [22]. There exists engineering applications for anemia classification including image analysis [23], statistical analysis [24], and clustering techniques [25], rule-based expert systems [26]- [28] , hybrid neural rule-based method [29]. One of the first efforts to diagnosis anemia using an ANN was performed by Brindorf et al. in 1996 in which a hybrid expert system of rule-based and ANN has been developed. They reported 96.5% precision in classifying microcytic anemia and the anemia of chronic disease. In 2002, Nikolaev and de-Menezes introduced an image recognition approach for classification of anemia [30]. Amendolia et al. in 2003 published a comparative study of K-nearest neighbor, support vector machine (SVM) and multi-layer perceptron (MLP) methods for thalassemia screening [31], [32]. The results revealed that MLP slightly improved the classification sensitivity (92%) in comparison with SVM (83%) while the classification accuracies were equal to 95%. Another comparison has been performed by Azarkhish et al. in 2012 [33]. They compare three different methods, an adaptive neuro-fuzzy inference system (ANFIS), an ANN and a logistic regression model in diagnosing iron deficiency anemia (IDA), and reported that the ANN is superior to the others. All of these methods have been successfully tested on various types of anemia. Moreover, to narrow the domain of anemia specifically on β-TT and IDA, an alternative diagnostic tool was compulsory. This study proposes such a tool. It compares the existing methods for discriminating β-TT and IDA, including traditional methods and the newer, AI-based diagnostic systems, with a novel meta-heuristic model. The main goal of this study is to reduce the diagnosis time and cost for β-TT and IDA subjects by increasing their discrimination precision through the analysis of CBC indices. This paper is organized as follows. Section 2 describes anemia and its effects on blood indices. Section 3 introduces the harmony search (HS) algorithm. In section 4, materials and methods are presented. Implementation details, experimental tests, and the results of the proposed method are discussed in section 5. Further tests and comparisons with other methods and models are presented in section 6. Finally, section 7 presents the study's conclusions.

II. EFFECTS OF ANEMIA ON CBC INDICES
There are several different, task-specific blood cell types, such as white blood cells (WBCs), red blood cells (RBCs), and platelets. A red blood cell transmits oxygen to all cells of the body and receives the CO2 excreted by them. These tasks are carried out by hemoglobin. Each red blood cell contains approximately 300 million molecules of hemoglobin. Hemoglobin consists of four polypeptide chains (globin) and heme. Hence, a change in the structure of globin affects the structure and functionality of the red blood cells. Each molecule of hemoglobin has four heme molecules. A heme molecule is a cofactor consisting of an Fe2+ (ferrous) ion in the center of a heterocyclic organic ring called porphyrin. The iron atom can bind an oxygen molecule through ion-induced dipole forces. In 1825, J.F. Engelhard discovered that the ratio of Fe to protein is identical in the hemoglobin of several species [34]. Iron deficiency anemia (IDA) leads to a reduction in the quantity and the quality of hemoglobin. There are four types of globin chains: alpha (α), beta (β), gamma (γ), and delta (δ). Based on the type of globin chains in its structure, the hemoglobin belongs to one of three major types, HbA (α2β2), HbA2 (α2δ2), or HbF (α2γ2). HbF, or fetal hemoglobin, is the main component of fetal red blood cells, which gradually diminishes in proportion to HbA after birth. The proportions of HbA, HbA2, and HbF in a mature, healthy person's blood are approximately 97%, 2-3%, and less than 1%, respectively. Based on the proportions of the hemoglobin types in a person, HBA (α2β2) is a major component of red blood cells. As a result, the most vital globin types are α-globin and βglobin. The first contains 141 amino acids regulated by genes on chromosome 16 and the second contains 161 amino acids regulated by genes on chromosome 11. Thalassemia results from defects on the genes that regulate the formation of globin chains. The type of thalassemia and the severity of the disease depend on the types and the number of defective genes, respectively. A person must have at least two defective genes of the same type to have the disease. If only one defective gene or two different defective genes are present, he has thalassemia trait. Persons with thalassemia trait do not have the disease but inherit genes that cause the disease. The type of defective genes determines whether the subject has α-thalassemia or β-thalassemia. Generally, βthalassemia is more common than α-thalassemia.

A. IRON DEFICIENCY ANEMIA (IDA)
As previously discussed, IDA leads to a reduction in heme. It is clear that a reduction in heme affects the concentration of hemoglobin in blood. Therefore, in a person with IDA, the blood parameters dependent on the hemoglobin concentration are affected. Some of these parameters, which are generally used as diagnosis indices, are: hemoglobin concentration (Hb), which denotes the quantity of hemoglobin in blood, red blood cell count (RBC), mean corpuscular volume (MCV), mean corpuscular hemoglobin (MCH), mean corpuscular hemoglobin concentration (MCHC), red blood cell distribution width (RDW), and hematocrit (HCT).

B. CHARACTERISTIC OF THALASSEMIA
In α-thalassemia or β-thalassemia, regardless of whether the subject has the disease or the trait, hemoglobin concentration becomes lower than normal. This reduction is a result of a decreased concentration of α-globin (α-thalassemia) or βglobin (β-thalassemia) chains. When the concentration of α or β-globin reduces, the quantity and quality of hemoglobin are naturally lower than normal. As a result, similar to an IDA subject, someone with thalassemia disease or trait has reduced Hb, MCV, MCH, MCHC, and HCT values. For thalassemia disease, diagnosis is not difficult because the indices are too low. But in thalassemia trait, the indices are only a little below the normal range, leading to difficulty discriminating between IDA and thalassemia trait. One of the most important differences between IDA and β-TT is the red blood cell count (RBC). Generally, to compensate for the lower Hb with β-TT, the RBC increases, while the RBC index is lower than normal for IDA. In summary, we can conclude that IDA and β-TT both lead to a reduced concentration of hemoglobin, but the reduction for β-TT is greater than for IDA. Despite these differing amounts of reduction, in some circumstances it is too difficult to discriminate between these types of anemia.

III. HARMONY SEARCH ALGORITHM (HS)
Mathematical based optimization algorithms exploit the benefit of substantial gradient data. In addition, there exist the wide range of programming approaches including Linear and non-linear programming to solve optimization problems [35]. The initial population should have been selected carefully for optimal convergence of the optimization algorithm. In contrast, Non-gradient meta-heuristic methods are other popular solutions which does not have inherent difficulties of gradient-based methods [35]- [43]. Generally, meta-heuristic algorithms comes from nature and are guided random search algorithm [43]. Harmony search (HS) is inspired from jazz musicians and has variety of applications in engineering introduced by Geem Z-W. in 2001 [35] and will be described in the remainder of this section.

A. AN OPTIMIZATION PROBLEM
Optimization problems either try to maximize the fitness or minimize the cost. An optimization problem which tries to minimize the cost function can be defined in Eq. (1) as follows: where ( ⃗⃗⃗ ) is the cost function, ⃗⃗⃗ is a vector describing the design parameters, and ( ) and ( ) are the lower and upper bound of the ' ℎ parameter, respectively.

B. EXPLANATION OF HARMONY SEARCH (HS) ALGORITHM
Each solution in HS algorithm is called a harmony and described as an n-dimensional vector. First, the initial population of HS should be randomly generated from the search space. Then, a new harmony is constructed by three major operators in HS namely, memory consideration, pitch adjustment, and random re-initialization. The algorithm proceeds by comparing the new harmony ( ⃗⃗⃗ ) with the existing worst harmony in the HM ( ⃗⃗⃗ ). ⃗⃗⃗ will be substituted by ⃗⃗⃗ provided that the new harmony has better cost than the worst harmony. The process of generating ⃗⃗⃗ and possible substitution with ⃗⃗⃗ continues until convergence occurs. HS has some parameters which should be preallocated at the beginning of optimization. These parameters which tune HS for optimized quality of convergence are harmony memory consideration rate (HMCR), harmony memory size (HMS), number of improvisations (NI), pitch adjustment rate (PAR), and bandwidth (BW).

C. INITIALIZING HARMONY MEMORY (HM)
Harmony memory (HM) is a matrix in which the rows indicate the harmonies. Suppose that ⃗⃗⃗ = ( (1), … , ( )) is the ' ℎ row (harmony) in HM. Then, each row of HM can be initialized by Eq. (2) as follows: Where is a random number in range of [0,1). Hence, HM is a matrix of HMS size with solution vectors in each row as shown in Eq. (3). (3)
Where is a random number in the range of [0, 1). 2 ∈ [0, 1) is a randomly generated number provided that ( ) was selected from HM. If 2 is less than pitch adjustment rate (PAR), an adjustment using Eq. (5) must be performed: Where is a random number in the range of [0,1) and ( ) is the corresponding bandwidth of the ' ℎ parameter in ⃗⃗⃗ .

E. UPDATING THE HARMONY MEMORY (HM)
The HM matrix should be updated upon creating a new harmony. HM rows are sorted ascendingly based on the cost. Thus, the last row of HM incorporates the worst harmony ( ⃗⃗⃗ ). The new recently created harmony, ⃗⃗⃗ , will be substituted by existing ⃗⃗⃗ if ⃗⃗⃗ has cheaper cost.in comparison with ⃗⃗⃗ . The sorting of HM matrix should be updated accordingly.

IV. MATERIALS AND METHODS
We have used a harmony search algorithm to minimize error in anemia classification. For this purpose, a new approximation method has been developed as the optimization tool for minimizing the approximator error, based on the harmony search algorithm (HS). The method borrows some principles from a successfully tested structure representation on genetic algorithms for nonlinear function fitting, or STROGANOFF, method [44]. The proposed method is explained in the following subsections.

A. THE PRINCIPLES OF STROGANOFF
In STROGANOFF, the input-output mapping can be described as a smooth multidimensional surface which plays the role of universal approximator [45]. The method uses a binary tree in which each non-leaf node (which is commonly a second-order polynomial) is modeled by a differentiable function. Additionally, each leaf node contains an input variable of the system. Fig. 1 shows a representation of the system.
Where the are the polynomial coefficients and the are the independent variables. The Kolmogorov-Gabor polynomial is a universal format for function modeling because it can be used to approximate any continuous function mapping to an arbitrary precision, if there are sufficient terms. For two independent variables, Eq. (8) with six terms of the Kolmogorov-Gabor polynomial, has been employed by various researchers [45].
( ) = 0 + 1 1 + 2 2 + 3 1 2 + 4 1 2 + 5 2 2 For the current optimization problem, the coefficients of the polynomials in non-leaf nodes of the binary tree ( Fig. 1) must be tuned up by fitting the output of function (̂) to the desired output ( ). Nikolay Y. Nikolaev and Hitoshi Iba [45] have proposed 16 quadratic forms presented in Table 1 to be applied on the non-leaf nodes of STROGANOFF algorithm instead of the complete bi-variable quadratic polynomial (Eq. 9).

B. CHOOSING THE INDICES OF CBC AS THE SYSTEM'S INPUTS
After reviewing CBC indices, the four indices Hb, RBC, MCV, and HCT are chosen as the system's inputs. Indices are chosen for minimum redundancy. To achieve this goal, we use the algorithm described in the two next subsections and some sample collected CBC specimens, drawn from those shown in Table 2. The chosen indices can populate the leaf nodes of a binary tree such as the one shown in Fig. 1. Additionally, these indices have been successfully used in various previous studies [2], [46].

1) DATA ANALYSIS
The dataset used in this work has been obtained from [47]. The total number of collected CBC tests is 750, obtained from the blood specimens of 390 males aged 20-35 and 360 females aged 17-32 years. The information is organized in three classes as follows [47]. Based on the data, the total numbers of β-TT and IDA subjects are 269 and 263, respectively. Table 2 shows some examples of different classes of collected samples. Abnormal values are shown in bold face ( Table 2). It is clear that the samples in class III are sufficiently similar and the relationship is sufficiently complex that a simple mathematical formula cannot distinguish reliably between IDA and β-TT subjects in this class. As the study's purpose is to differentiate between IDA and β-TT subjects, we excluded the normal samples and the remaining samples (532 samples including IDA and β-TT) were randomly divided into the training and the testing sets of 132 and 400 samples, respectively.

2) AN ALGORITHM FOR ELIMINATION REDUNDANT INDICES
The elimination of redundant CBC indices from the input set leads to a more efficient system. Fig. 2 shows the algorithm which is employed for choose the most reliable indices for discrimination. this algorithm is called Pattern Based Index Selection (PBIS). The name is chosen because the elimination is conducted based on the similarity between indices with respect to their abnormal patterns, as seen in the Table 2.
As seen in Fig. 2, for each iteration of the algorithm, the index such has the least similarity with the other indices, with respect to their abnormality pattern, must be found. Eq. (10) can be used to calculate the Coefficient of Similarity (COS) between any two CBC indices.
Where is the number of samples, is the value of the ' ℎ CBC index for the ' ℎ sample, and ( ) and ( ) are functions that denote whether is below the lower bound or above the upper bound of the normal range for the ' ℎ CBC index, respectively.
After applying the above algorithm, four indices, RBC, Hb, HCT, and MCV, are chosen as system inputs.

C. PROPOSED METHOD
The proposed method is based on harmony search and certain genetic programming principles. Before defining a harmony, the elements of the binary tree ( Fig. 1) must be defined. The four indices of CBC chosen in the previous section populate the leaf nodes. Non-leaf nodes take the complete bi-variable polynomial that comes from the Kolmogorov-Gabor polynomial with six coefficients (Eq. 9). There are three non-leaf nodes in the tree and the polynomial shown in (Eq. 9) is applied to each of them. Hence, there are 18 parameters that must be tuned up to reach the optimum value. Because is a symmetric function, we can generate 15 different schemes for the tree (Fig. 1) and the four CBC indices chosen as input variables. Table 3 lists the possible schemes for the tree.       A one-dimensional array is assigned to the harmony, containing 19 elements. Fig. 3 shows the fields of a harmony for the proposed method. Therefore, a harmony describes the structure of a tree, which can be one of the 15 different schemes in Table 3, and the coefficients of that tree. To determine the most suitable scheme for this problem, the harmony search algorithm can use a heuristic search to find the best structure for the system tree using the 'Id' field of the harmony and evaluating its corresponding tree to optimize the coefficients of the polynomials of the tree.

D. IMPLEMENTING THE PROPOSED METHOD
We use the programming language C# to implement the proposed method. A user interface is provided for quick testing and tuning the algorithm's parameters. The interface gives the user the ability to set various parameters, including the size of harmony memory (HMS), the initial values for HMCR and PAR, the number of improvisations (NI), and the choice of standard harmony search or dynamic harmony search (DHS). Additionally, the interface monitors the outputs of the program.

1) INITIALIZE THE HARMONY MEMORY
The harmony memory is initialized randomly by Eq. (2). The default value for HMS is 100, but it can be changed by the user. Hence, the harmony memory is a matrix of floating-point numbers of size 100×20. The elements of column 20 of the HM contain the fitness values for of each row/harmony. As shown in Fig. 3, a harmony is a one-dimensional array of size 19 where the first 18 elements contain the polynomials' coefficients and the last element contains the Id of the tree scheme. Hence, the values for the 18 first elements can be determined by trial and error, but the last element must be an integer in range of [0,14]. Therefore, we can calculate the values for the first 18 elements of harmony by Eq. (2), but the last element is obtained by Eq. (13) as follows: ( ) = (15 × ), = 1. . , = 19 (13) Where is a random number in range of [0, 1). The values for column 20 of the HM are determined by a fitness function described in the next section. Next, the HM rows are sorted in descending order with respect to their fitness (column 20). The harmonies in the HM are now organized best to worst from top to bottom.

2) CALCULATING THE FITNESS FOR EACH HARMONY
To determine the best harmony for the discrimination of IDA from thalassemia trait, suppose that and are the averages for the tree's outcome when the CBC indices of the samples with IDA and thalassemia trait are applied, and and are the variances of those values, respectively. It is clear that a function with a larger distance between averages and smaller variances is a better discriminator. Therefore, we use Eq. (14) as the fitness function for the harmony.

3) USING DYNAMIC HARMONY SEARCH (DHS) INSTEAD OF STANDARD HS
Dynamic harmony search combined with tuning the parameters leads to faster convergence [48]. In the standard HS, the parameters of the algorithm such as HMCR and PAR have fixed values, but in the DHS, the parameter HMCR changes with Eq. (15) as follows: where and are the normalized fitness/cost average of ten recently considered memory solutions and ten recently generated random solutions respectively in maximum finding. However, In case of minimum finding, the definitions of solutions are contrariwise. This improvement in the HS algorithm is tunned to be occurred every ten iterations. Likewise, variable (which is derived from memory consideration) is used to evaluate the adjustment of new solutions for improving the pitch adjustment rate (PAR) as stated in Fig. 4. Variable can be further used in Eq. (16) to improve the parameter PAR for the proposed harmony search algorithm.
= + ×(1− ) 100 (16) Fig. 5 illustrates the improving process of PAR for every 10 recently generated memory solutions. For comparison between standard HS and DHS, the user interface can perform the proposed method with both algorithms optionally. As can be seen in Fig. 4, the standard HS and DHS both have converged to the optimal value, but the standard HS needs more time than the dynamic harmony search.

E. TESTING THE PROPOSED METHOD
To test the proposed method, 750 CBC samples have been obtained from [47]. After eliminating unwanted samples, 532 CBC samples remain for further examination. Of the remaining subjects, there are 269 IDA and 263 β-TT subjects. Six popular medical metrics, sensitivity (SENS), specificity (SPEC), positive predictive value (PPV), negative predictive value (NPV), accuracy (ACC) and Youden's index (YI), have been used in all of the experiments. These metrics are obtained by Eqs. (17)(18)(19)(20)(21)(22) Where TP, TN, FP and FN are true positive, true negative, false positive and false negative, respectively. Table 4 shows the results of applying the method to the test samples.

V. COMPARISON
This section compares the best results of the proposed method against well-known works in the literature. This section is organized in two subsections. The first compares the proposed method with some traditional and some more recent methods on the same dataset. The second subsection compares the performance of the proposed method with the information reported by other researchers.

A. Traditional method for discrimination between IDA and β-TT
There are several mathematical methods that traditionally have been used for the IDA/ β-TT discrimination problem. Some of these methods are shown in Table 5. For better comparison, the traditional methods were tested on the same CBC samples used for the proposed method. The results for these methods are shown in Table 6. The results in Table 6 show that the proposed method is significantly more accurate than traditional methods. The range of accuracy for the traditional methods, which are based on a simple mathematical formula, is in (78, 93.2), while the proposed method has an accuracy of approximately 98%.

B.Artificial intelligence-based models for classification anemia
There are several AI-based models for the diagnosis and classification of anemia. The second set of experiments compares the performance of the proposed method with other, similar methods. As can be observed in Table 7, the proposed method is more accurate than these other similar methods.
Although the results for [50] seem to be better than this work, but the method which is used by [50] get thee image of blood cells as the system input instead of CBC, which needs additional laboratory test on patient. This Work 1 adaptive neuro-fuzzy inference system; 2 artificial neural network; 3 multi-layer perceptron; 4 support vector machine; 5 Knearest neighbor; 6 pattern-base input selection artificial neural network; 7 radial basis function; 8 probabilistic neural network; 9 dynamic harmony search ;

VI. CONCLUSIONS
In this study, a new method based on dynamic harmony search (DHS) was proposed for the discrimination of iron deficiency anemia (IDA) and β-thalassemia trait (β-TT). The method is implemented in C# and has been successfully tested using collected CBC sample data as input. Choosing the most suitable CBC indices to use as the system input is performed by a pattern-based index selection (PBIS) algorithm. The proposed method has been trained on 132 CBC samples. The results indicate that the proposed method, with an accuracy of approximately 98%, outperforms the other methods in the literature. The existing artificial neural network (ANN), MLP, and ANFIS methods, respectively, show the nearest performance on the anemia classification problem.