Chromatography in itself is not an accurate analytical technique, but rather a separation one. The identification of oxygencontaining organic compounds can be made with the method of gas chromatographic peak in comparison with that of a standard sample of each compound. Because samples of pure compounds are not always available, it is important to develop QSRR that can efficiently predict retention parameters by using theoretical descriptors computed from chemical structure.

Quantitative Structure-Retention Relationships (QSRR) [1] establish the relationship between a chemical structure and its chromatographic retention value, which has been demonstrated to be a powerful tool for the investigation of chromatographic parameters. The main advantage of QSRR is the ability to distinguish in quantitative theoretical terms, packing materials of different chemical nature of the organic ligand and/or organic or inorganic support [2], Furthermore, it can be of valuable assistance in the prognosis of the behavior of new molecules, even before they are actually synthesized [3].

An important property that has been extensively studied in QSRR is the chromatographic retention index. The retention index is a generally accepted type of data used for the identification of chemical compounds by gas chromatography. A retention index is a continuous quantitative variable that relates the retention of a solute to the retention of a set of standard compounds. Retention indices are much less dependent on experimental factors (e.g., Temperature, flow, column, length etc.) than retention times. While Kovats retention indice [4] have linear collerations with column temuprature. And they were obtained by the logarithmic interpolation method.

QSRR on the Kovats retention indices have been reported for different types of organic compounds. The Kovats retention index is the most popular dependent variable in QSRR studies because of its reproducibility and accuracy. In many cases, the precision and accuracy of the QSRR models are not sufficient for identification purposes; still the models are useful to elucidate retention mechanisms, to optimize the separation of complex mixtures or to prepare experimental designs.

Topological descriptors computed on the basis of molecular graph are easy to be calculated with present computing facilities. Due to the simplicity and efficiency of graph-theoretical approaches, we take novel polarizability effect index (PEI), odd-even index (OEI), the sum eight values X_1CH of every C-H bond adjacency matrix S_x1CH.

An interesting and increasing application of QSRR is to test various chemometric methods from multiple linear regression (MLR) methods to Artifical neural network (ANN) methods. Multiple linear regression (MLR) is without doubt the most frequently applied technique in building QSRR models.

Gene Expression Programming (GEP) is a new evolutionary algorithm that evolves from computer programs (they can take many forms: mathematical expressions, neural networks, decision trees, polynomial constructs, logical expressions, and so on). The computer programs of GEP, irrespective of their complexity, are all encoded in liner chromosomes. Then the liner chromosomes are expressed or translated into expression trees (the branched structures). Thus, in GEP, the genotype (the liner chromosomes) and the phenotype (the expression trees) are different entities (both structurally and functionally), and because of this apparently trivial fact, this new evolutionary system can finally make a difference, successfully assisting researchers in the design of robust and accurate computer models [5].

Gene Expression Programming (GEP) is a new evolutionary algorithm that evolves from computer programs (they can take many forms: mathematical expressions, neural networks, decision trees, polynomial constructs, logical expressions, and so on). The computer programs of GEP, irrespective of their complexity, are all encoded in liner chromosomes. Then the liner chromosomes are expressed or translated into expression trees (the branched structures). Thus, in GEP, the genotype (the liner chromosomes) and the phenotype (the expression trees) are different entities (both structurally and functionally), and because of this apparently trivial fact, this new evolutionary system can finally make a difference, successfully assisting researchers in the design of robust and accurate computer models [5].

Data set

The Kovats retention index of 91 molecules (include esters, ketones, and alcohols) taken from reference [6] were presented in Table 1. Kovat’s retention index of all compounds was obtained under the same conditions on two stationary phases: OV-1 (dimethylpolysiloxane) and SE-54(5% phenyl -95% dimethylpolysiloxane) 74 molecules were used as training set for model generation and 17 molecules were used as test set for model prediction. The corresponding experimental and predicted values of the RI for all the molecules studied in this work are shown in Table 1.

Table 1: Data set and corresponding experimental (exp.) and predicted (cal.) values of RI.

GEP theory

Gene Expression Programming (GEP) was first proposed formally by Candida Ferreira in 2001. It was an elegant and efficient solution to expression-mutation problems. GEP, which is an extraordinarily powerful tool, is a subset of Genetic Algorithms, except it uses genomes whose strings of numbers represent symbols. GEP-an evolutionary algorithm inherits both the evolutionary simplicity of Genetic Algorithms (GA) and the expressional power in Genetic Programming (GP) by utilizing a genotype/phenotype representation system. The string of symbols can further represent equations, grammars, or logical mappings.

Ferreira [5] proposes the use of a set of genetic operators: Replication, Mutation, IS Transposition, RIS Transposition, Gene Transposition, 1-Point Recombination, 2-Point Recombination, Gene Recombination. As Ferreira comments, the advantages of a Genetic Representation like the one in GEP are simple entities: linear, compact, relatively small, easy to manipulate genetically. The genetic operators applied to them are less restricted than those used in GP [5].

Fortunately for us, in GEP, thanks to the simple rules that determine the structure of expression trees and their interactions, it is possible to infer immediately the phenotype given the sequence of a gene. It is easy for a computer program to follow these three rules while performing mutations, and it never has to check whether the resulting expression has valid syntax. By allowing a broad range of mutations, the process can efficiently explore a high dimensional space, and the expressions can change in size as functions are replaced by terminals and terminals by functions.

GEP are evolutionary tools inspired in the Darwinian principle of natural selection and survival of the fittest individual and uses populations of candidate solutions to a given problem in order to evolve new ones. These methods use an initial random population and apply genetic operators to this population until the algorithm finds an individual that satisfies some termination criteria. The evolving populations undergo selective pressure and their individuals are submitted to genetic operators.

Gene representation: GEP genes are composed of a head and a tail. The head contains symbols that represent both functions (elements from the function set F) and terminals (elements from the terminal set T), whereas the tail contains only terminals. Therefore, two different alphabets occur at different regions within a gene. For each problem, the length of the head h is chosen, whereas the length of the tail t is a function of h and the number of arguments of the function with the most arguments n, and is evaluated by the equation:

t=h (n-1) + 1 (1)

Consider a gene composed of {Q, *, /, -, +, a, b}. In this case n=2.

For instance, for h=15 and t=16, the length of the gene is 10+11=21. One such gene is shown below (the tail is shown in bold):

0123456789012345678901234567890

*b+a-aQab+//+b+babbabbbababbaaa (2)

It codes for the following ET:

A K-expression can be mapped into an expression tree (ET) following a first-order procedure. A branch of the ET stops growing when the last node in this branch is a terminal. For example, the ET shown above corresponds to chromosome (2). In this case, the open reading frames end at position 7, whereas the gene ends at position 30.

Chromosomes in GEP are usually composed of more than one gene of equal length. For each problem or run, the number of genes, as well as the length of the head, is chosen. Each gene codes for a sub-ET and the sub-ETs can be linked by pre-defined rules forming a more complex multi-subunit ET.

Selection method and genetic operators: References [3] suggests there is no difference between different selection methods. It is strongly advised to use a simple elitism in any GEP implementation. The elitism means copying the best (or few best) individual to the offspring population without modifying them. GEP uses the well-known roulette-wheel method for selecting individuals.

In GEP, individuals were selected according to fitness by roulette wheel sampling coupled with the cloning of the best individual. The fitter the individual is, the higher the probability of leaving more offspring. Thus, during replication the genomes of the selected individuals are copied as many times as the outcome of the roulette. The roulette is spun as many times as there are individuals in the population, always maintaining the same population size.

GEP uses simple elitism of the best individual of a generation, preserving it for the next one. Replication is an operation that aims to preserve several good individuals of the current generation for the next one. In fact, this is a do-nothing probabilistic operation that takes place during selection (using the roulette-wheel method), and replicated individuals will be subjected to the action of the genetic operators. The mutation operator aims to introduce random modifications into a given chromosome. A particularity of this operator is that some integrity rules must be obeyed so as to avoid syntactically invalid individuals. In the head of a gene, both terminals and functions are permitted (except for the first position, where only functions are allowed). However, in the tail of a gene only terminal is allowed.

Mutation, Inversion, Transposition and Recombination

Mutation: Mutations can occur anywhere in the chromosome. Simple mutation just replaces symbols in genes with replacement symbols. However, the structural organization of chromosomes must remain intact. Symbols in the heads of genes can be replaced by functions or terminals (variables and constants). Symbols in the tail sections can be replaced only by terminals. Randomly change symbols in a chromosome. Symbols in the tail of a gene may not operate on any arguments. Typically two-point mutation per chromosome is used. It is worth noticing that in GEP there are no constraints neither in the kind of mutation nor the number of mutations in a chromosome: in all cases the newly created individuals are syntactically correct programs.

Inversion: Inversion reverses the order of symbols in a section of a gene. A portion of a chromosome is chosen to be inserted in the head of a gene. The tail of the gene is unchanged. Thus symbols are removed from the end of the head to make room for the inserted string. Typically a probability of 0.1 of insertion is used.

Transposition: Transposition selects a group of symbols and moves the symbols to a different position within the same gene. Gene transposition moves entire genes around in the chromosome. One gene in a chromosome is randomly chosen to be the first gene. All other genes in the chromosome are shifted downwards in the chromosome to make place for the first gene.

An IS element is a variable-size sequence of elements extracted from a random starting point within the genome (even if the genome was composed of several chromosomes). Another position within the genome is chosen as the insertion point.

This target site must be within the head part of a gene and cannot be the first element (gene root). The IS element is sequentially inserted in the target site, shifting all elements from this point onwards and a sequence with the same number of elements is deleted from the end of the head, so that the structural organization is maintained. This operator simulates the transposition found in the evolution of biological genomes. RIS is similar to the IS transposition, except that the insertion sequence must have a function as the first element and the target point must be also the first element of a gene (root).

The transposable elements of GEP are fragments of the genome that can be activated and jump to another place in the chromosome: (1) Short fragments with a function or terminal in the first position that transpose to the head of genes, except to the root (insertion sequence elements or IS elements); (2) Short fragments with a function in the first position that transpose to the root of genes (root IS elements or RIS elements); (3) Entire genes that transpose to the beginning of chromosomes.

Recombination: During recombination, two chromosomes are randomly selected, and genetic material is exchanged between them to produce two new chromosomes.

The cross over operation this can be one point (the chromosomes are split in two and corresponding sections are swapped), two point (chromosomes are split in three and the middle portion is swapped) or gene (one entire gene is swapped between chromosomes) recombination. Typically the sum of the probabilities of recombination is 0.7.

In GEP there are three kinds of recombination: one-point, twopoint, and gene recombination. (1) One-point: During one-point recombination, the chromosomes crossover a randomly chosen point to form two daughter chromosomes; (2) Two-point: In twopoint recombination the chromosomes are paired and the two points of recombination are randomly chosen. The material between the recombination points is afterwards exchanged between the two chromosomes, forming two new daughter chromosomes; (3) Gene recombination: recombines entire genes. This operator randomly chooses genes in the same position in two parent chromosomes to form two new off springs. In gene recombination an entire gene is exchanged during crossover. The exchanged genes are randomly chosen and occupy the same position in the parent chromosomes. It is worth noting that this operator is unable to create new genes: the individuals created are different arrangements of existing genes.

Fitness function: A fitness function is the most important part of any EA application. Fitness function given with above equations allows for fulfilling all of the set conditions. In GEP, fitness is based on how well an individual model the data. If the target variable has continuous values, the fitness can be based on the difference between predicted values and actual values. Evolution stops when the fitness of the best individual in the population reaches some limit that is specified for the analysis or when a specified number of generations have been created or a maximum execution time limit is reached.

All of the fitness functions produce fitness scores in the range 0.0 to 1.0 with 1.0 being ideal fitness – that is, the individual exactly fits the data. If a function is unviable – for example, it takes the square root of a negative number or divides by zero – then its fitness score is 0.0.

GEP evolution process: The GEP evolution begins with the random generation of linear fixed-length chromosomes for individuals of the initial population. The chromosomes are translated into ETs and subsequently into mathematical expressions, and the fitness of each individual is evaluated based on a pre-defined fitness function. The individuals are then selected by fitness to reproduce with modification. The individuals of this new generation are, in their run, subject to the same developmental process. The selection and reproduction is accomplished by roulette-wheel sampling with elitism, which guarantees the survival and cloning of the best individual to the next generation. Variation in the population is introduced by applying one or more genetic operators to selected chromosomes, including crossover, mutation and insertion.

Models

GEP model: The GEP program was coded by the combination of MATLAB and VC++. The MATLAB software has the advantage of computing matrix conveniently and programming efficiently, but its operating efficiency is relatively low. So VC++ was combined for its powerful function and the characteristics of higher operating efficiency with MATLAB. In this paper, MATLAB engine was used to achieve the combination with VC++ programming. There are two steps: (1) Add MATLAB engine library header files and library functions of the path. (2) Add libmx.lib libeng.lib libmex.lib to complete the import of the corresponding MATLAB engine static link library.

From the data in Table 1, GEP method was used that 6 topological index as input, output for its retention index. During the run, parameter values were needed to adjust constantly in order to achieve the optimal results. The set of optimal parameter values were listed in Table 2 and the predicting results of test set on OV-1 and SE-54 were listed in Figures 1 and 2. It can be seen from the figures that the predictive values of gas chromatography retention index of oxygen-organic compounds were in good agreement with the experimental data.

Table 2: Parameters of GEP models.

Figure 1: Plot of the predicted RI against the experimental values on OV-1 column for test set based on GEP.

Figure 2: Plot of the predicted RI against the experimental values on SE-54 column for test set based on GEP.

ANN model: Non-linear statistical treatment of QSRR data is expected to provide models with better predictive quality as compared with related MLR models. In this perspective, functioning and applications of ANN have been adequately described elsewhere [7-10]. Extensive use of ANN, which has inherent ability to incorporate nonlinear and cross-product terms into the model and does not require prior knowledge of the mathematical function as well, largely rests on its flexibility and less sensitivity to collinearity among variables. The theory behind ANN and their use in chromatography have been reported elsewhere [11-13].

Multi-layer feed forward networks, with good self-learing ability and adaptability is widely used in the field of QSRR modeling [14]. Commonly, they consist of three layers: one input layer formed by a number of neurons that equal to the number of descriptors, one out neuron (providing the model response) and a number of hidden neurons fully connected to both input and out neurons. Among the available learning algorithms, back-propagation of errors is one of the most widely used [8,15].

Usually, there are four steps involved in ANN modeling: (1) assembling the training data of input (independent variables) and output (dependent variables), (2) deciding the network architecture, (3) training the network, and (4) simulating the network response to new inputs. The training process is simply an optimization process which aims at finding the set of weight and biases associated with each layer that will minimize the error objective function related to the deviations of the network predictions from the true response output data of the training set.

Before data set was used for the training of ANN, it was normalized separately. Its minimum value was set to zero and maximum to one. The proper number of nodes in the hidden layer was determined by training the network with different number of nodes in the hidden layer. The root-mean-square error (RMSE) value measures how good the outputs are in comparison with the target values. In this paper, following a troubleshooting study to investigate the effects of the number of hidden layers and the number of neurons involved in these hidden layers, a 2-3-1 network, with tansig-logsig transfer functions, was found to be the most optimum in terms of the root mean squared errors (RMSE) obtained.

ANN with basic back-propagation of errors learning algorithm was used in this study to predict oxygen-containing retention index. A three-layer network with a sigmoid transfer function was designed for ANN. The ANN program was coded in MATLAB 7 for windows [15].

The MLR: For regression analysis, data set was randomly divided into two groups: training and test sets. The training set, composed of 74 molecules, was used for the model generation. The test set, composed of 17 molecules, was used to evaluate the generated model. The program used for MLR analysis was compiled in Statistical Product and Service Solutions (SPSS version 19.0 IBM) software. In MLR analysis, in order to minimize the information overlap in descriptors and to reduce the number of descriptors required in regression equation, the concept of non-redundant descriptors was used in this study. The best equation was selected on the basis of the highest multiple correlation coefficients (R) and the lowest root mean squared error (RMS). The linear equation between these descriptors and the retention parameters of fluid catalytic cracking (FCC) gasoline was:

(1)

Where b₀ is the intercept and b_j is the regression coefficient for descriptor j. The statistical results obtained by using the two molecular descriptors based on MLR are listed in Table 3 and plotted against the experimental values in Figures 3 and 4.

Table 3: Model parameters value and coefficients for MLR model.

Figure 3: Plot of the predicted RI against the experimental values on OV-1 column for test set based on MLR.

Figure 4: Plot of the predicted RI against the experimental values on SE-54 column for test set based on MLR.

It is common to consider four statistical parameters for regression equation. These parameters are the number of descriptors, correlation coefficient (R) for training and test sets, root mean squared error (RMS) for training and test sets, and F statistic. A reliable MLR model is one that has high R and F values, low RMS and least number of descriptors. In addition to these, the model should have a high predictive ability. Consequently, among different models, the best model was chosen, whose specifications are presented in Table 3. Here the corresponding descriptors used in MLR were applied as inputs for ANN in order to compare the performance of the two models.

The main aim of the present work was developing a QSRR model to predict the retention parameter (RI) of oxygen-containing compounds appeared in Table 1. A linear model of MLR was developed, whose specifications are given in Table 3. All statistic tests were performed at a significance level of 5%. MLR model performance was measured by three metrics: (1) R, which gives the fraction explained variance for the analyzed set, was used to measure the model’s fit performance. (2) Root Mean Squared error (RMS), which can give the bias in the prediction, was used to evaluate the model’s predictive precision: the lower the RMS, the better the prediction precision. It can be calculated as below:

(3)

where d_i is the target value, o_i is the experimental value and n is the number of compounds in analyzed set. (3) The variance ratio of calculated and observed activities F.

After the linear model was gained, non-linear characteristics of the retention parameter were also performed using ANN. Here a feedforward neural network with basic error back-propagation algorithm was constructed to model the nonlinear QSRR models. Therefore, a 2-3-1 BP-ANN, with tansig-logsig transfer functions, was developed. Figure 5 demonstrates the plot of the ANN predicted versus the experimental values of the RI for the data set. A correlation coefficient of this plot indicates the reliability of the model. As can been seen in Table 4, the correlation coefficient R on OV-1 and SE-54 for the ANN models are larger than that of MLR models respectively, which indicates that the ANN models are slightly improved to MLR models. The residuals of calculated values of RI by ANN are plotted against the experimental values in Figure 6. The propagation of the residuals on both sides of zero line indicates that no symmetric error exists in the development of ANN model.

Figure 5: Plot of the predicted RI against the experimental values on OV-1 column for test set based on ANN.

Figure 6: Plot of the predicted RI against the experimental values on SE-54 column for test set based on ANN.

Table 4: Result of correlation coefficient (R) with MLR, ANN and GEP for the test set.

Correspondence author is grateful to the financial support from National Natural Science Foundation of China, Lanzhou University of China and Liaoning Province Education Department of China, which made this work possible.