Modelling the effect of structural QSAR parameters on skin penetration using genetic programming

The data set used in this work was taken from a study of the effect of structural QSAR parameters on skin penetration. In [3], Ghafourian and Fooladi analysed 20 variables of 17 inputs and 3 outputs—permeability (^log  K ^_p ), ^skin/^water partition (^log  K ^_m ) and diffusion coefficient (^log  D/h), as listed in table 1. Ghafourian and Fooladi analysed 39 compounds. However, in this study, for ease of comparison with the statistical results in the literature, only 22 compounds (compound 1–14 and 26–33) were used for evaluating the predictability of GP for all responses. From a total of 22 compounds used for observation, 20 records are used as training data and 2 records are used as test data selected by the Smart Selection function in GP.

GP has been known primarily as a learning technique producing mathematical equations [10–12]. It is founded on the basic principle of Darwin's theory of evolution. In nature, biological structures that are more successful at operating in their given environment have a higher probability of surviving and reproducing. These structures are considered to be the result of Darwinian natural selection operating in an environment over a period of time. With a similar principle of natural selection, Koza introduced GP that attempts and succeeds in applying this evolutionary theory in order to find the best (fittest) or the most appropriate equation (solution) for a problem domain [10].

2.2.1. Representation scheme

GP modifies individuals of the population after a number of runs in every cycle, called a generation, to find the best solution. The structure of individuals in GP is a tree (see figure 1(a)). Possible structures in genetic programming are the set of all compositions of functions and the set of terminals.

Zoom In Zoom Out Reset image size

The elements in the function set may include:

• arithmetic operations (+, -, *, etc),
• mathematical functions (such as sin, cos, exp and log),
• Boolean operations (such as AND, OR and NOT),
• conditional operators (such as If-Then-Else),
• functions causing iteration (such as Do-Until),
• functions causing recursion,
• any other domain-specific functions that may be defined.

The terminals are typically either a variable or a constant value.

2.2.2. Genetic operations

2.2.3. Fitness function

2.2.4. Implementation scheme

The software used was a modified form of that described previously [14], but with additional exponential and other mathematical functions as well as the fitness functions listed. Table 2 lists the values of the control parameters used in the present study.

In order to evaluate the quality of a model generated by GP, the correlation coefficient R-squared (R²) was computed, with higher values of R² indicating improved quality of the model,

where is the mean value of the dependent variable y, is the predicted value from the model and n is the number of records.

The permeability coefficient has the following expression:

From equation (1), the R² value of the predictive model is considered as a satisfactory result and moreover the predictive equation is extremely simple for the permeability coefficient (^log  K ^_p ) response. In equation (1), only three parameters—the number of hydrogen bonding heteroatoms [OT], the sum of atomic charges on the hydrogen bonding heteroatoms [q ^{_s −b} ], and the second-order molecular shape index [²κ ]—from a total of 17 input variables influence the permeability coefficient ^log  K ^_p .

The values of the training parameters, which GP used to produce equation (1), are shown in table 3. For this response, the simplicity of the predictive equation resulted from a maximum of 10 nodes in the training trees. However, because of this small number of tree nodes, the population size required was 10 000 with 500 generations. In addition, maximum values of crossover and mutation in the training models were needed to obtain an acceptable model. For this response, the AIC fitness function was used to generate a high-quality model. The details of the equations from the statistical method are shown below:

As analysed in the literature using MINITAB statistical software [3], equation (2) shows that is correlated with the logarithm of total solvent accessible surface area calculated using a probe of 1.4 Å radius (log [TA]), the sum of atomic charge on the hydrogen bonding heteroatoms [q ^{_s −b} ]and the number of hydrogen bonding heteroatoms [OT] in equation (2). In equation (3), this response is also correlated with the number of hydrogen bonding heteroatoms [OT], the second-order molecular shape index [²κ] and the sum of atomic charge on the hydrogen bonding heteroatoms calculated by the AM1 method [q ^{_s +a} ].

When comparing the predictive models between GP and statistical methods, the results of regression analysis for both methods using whole data of 22 structures in figure 2 indicated that there are only minor differences between these techniques. While producing a slightly lower linear R² value, GP generated a better model with more satisfactory values of slope and intercept coefficients for this response.

Zoom In Zoom Out Reset image size

The ^skin/^water partition coefficient is determined by the following formula:

For this response from GP given in equation (4), the parameters of the octanol/water partition coefficient molecular weight [^log  P ^_oct ], the second-order molecular shape index [²κ] and the energy of the highest occupied molecular orbital E ^_HOMO correlate with the skin/water partition coefficient (^log  K ^_m ) response with a high R² value of 0.96. This result is similar to the prediction for the first response of these QSAR data in that the equation predicted is simple with only three parameters.

As seen in table 4, the AIC fitness function also gave a high quality model with a population size of 10 000 500 generations, as well as the same values of crossover and mutation as in the case for the first response. However, in order to achieve a satisfactory result for (^log  K ^_m ), a larger number of tree nodes equal to 15 are required, and the value of the reproduction parameter also increased from 0.05 to 0.2.

For the skin/water partition coefficient ^log  K ^_m response, the models predicted by the statistical method are

As reported in the literature, the [E ^_HOMO ] parameter is a hydrogen bonding acceptor ability parameter, and according to equation (5) it has a positive effect on (^log  K ^_m ). Also with regard to the literature reported correlation of [^log  P ^_oct ] with the two parameters [OT] and [log TA] in equation (7) [3], it can be seen that the parameters of [OT] and [log TA] influence the octanol/water partition coefficient [^log  P ^_oct ]. The GP (4) for this response shows that the predictive model also displays the correlation between (^log  K ^_m ) and two parameters of [E ^_HOMO ] and [^log  P ^_oct ].

When equations (6) and (4) are considered in more detail, the results from the two methods are consistent in that the skin/water partition coefficient is correlated with the size of the compounds, as presented by [log TA] (in equation (6)) and [²κ] (in equation (4)). Thus GP also produced a high- quality model for this response data with similar influential parameters as reported in the literature [6].

In order to examine the predictability of GP for the skin/water partition coefficient (^log  K ^_m ), regression analysis was applied, and the results of both methods for this response were compared. From figure 3, the results of regression analysis for both methods using whole data of 22 structures show that there are no major differences between the two methods in this case. The linear R² values of three equations and regression coefficients observed are also similar.

Zoom In Zoom Out Reset image size

The diffusion coefficient ^log  D/h is predicted by GP as follows:

It can be seen from equation (8) that the predictive model for this response has extremely high train and test R² values with four input variables influencing the final results of this response. The diffusion coefficient (^log  D/h) response is correlated with the number of hydrogen bonding hydrogens [HT], the sum of atomic charges on the hydrogen bonding heteroatoms calculated by both methods of the AM1 method [q ^{_s −a} ]and the classical method [q ^{_s -b} ] [3], and the most negative electrostatic potential on the solvent accessible surface of molecule [ESP] parameters.

GP required the same population size and generations as in the predictions for the previous responses of 10 000 and 500, respectively (see table 5). Apart from the training model using the MDL fitness function and a larger number of nodes set to 15, all values of the remaining parameters in the training model are the same for the GP prediction for the first response of this data set.

In the original paper, Ghafourian and Fooladi [3] showed that this response could also be given by

Based on equation (9), an alternative equation for this response predicted by GP using the predictions from the ^log  K ^_m and ^log  K ^_p responses was

As reported in the literature [3], the response of the diffusion coefficient ^log  D/h was modelled successfully using the statistical method with very high R² values, as shown below:

The model predicted by GP in equation (8) demonstrates that there are four parameters from a total of 17 influencing the diffusion coefficient (^log  D/h).

When the predictive model of GP (8) is compared with the statistical results (equations (11) and (12)), both methods indicate that the diffusion coefficient (^log  D/h) is negatively correlated with the number of hydrogen bonding hydrogens [HT] (those connected to the oxygen or nitrogen atoms), and positively correlated with the most negative electrostatic potential on the solvent accessible surface of molecules [ESP], which can be considered as the hydrogen bonding basicity parameter [3]. However, in the statistical result, the molecular weight [MW] influences (^log  D/h), while in the GP result, ^log  D/h is correlated with the sum of atomic charges on the hydrogen bonding heteroatoms [q ^{_s −a} ] and [q ^{_s −b} ]parameters.

From figure 4 comparing equation (8) from GP to equations (11) and (12), it is apparent that GP delivers a slightly higher quality model for parameter (^log  D/h) in comparison with the statistical results based on the regression analysis for whole data of 22 records. Although the linear R² values of the three models are equal to 0.98, the slope and intercept coefficients in the linear equations of GP are closer to 1.00, while those from the statistical equations are 0.00.

Zoom In Zoom Out Reset image size

From figure 4, it can also be seen that equation (10) is an improved model with more satisfactory values of slope and intercept coefficients in comparison with the statistical results.

In general, for the diffusion coefficient (^log  D/h) response, while not generating the same influential parameters as reported in the literature, the quality of prediction of GP for this response is high in both observations on the R² value as well as on the regression analysis.