ETM-ANN Approach Application for Thiobenzamide and Quinolizidine Derivatives

The structure anti-influenza activity relationships of thiobenzamide and quinolizidine derivatives, being influenza fusion inhibitors, have been investigated using the electronic-topological method (ETM) and artificial neural network (ANN) method. Molecular fragments specific for active compounds and breaks of activity were calculated for influenza fusion inhibitors by applying the ETM. QSAR descriptors such as molecular weight, EHOMO, ELUMO, ΔE, chemical potential, softness, electrophilicity index, dipole moment, and so forth were calculated, and it was found to give good statistical qualities (classified correctly 92%, or 48 compounds from 52 in training set, and 69% or 9 compounds from 13 in the external test set). By using multiple linear regression, several QSAR models were performed with the help of calculated descriptors and the compounds activity data. Among the obtained QSAR models, statistically the most significant one is the one of skeleton 1 with R2 = 0.999.


Introduction
Influenza (also known as flu) virus is an etiologic cause of widespread epidemics of acute respiratory tract infections that occur annually during the winter season [1,2]. Less frequent pandemic episodes of influenza are associated with greater mortality and morbidity, usually a consequence of genetic recombination of viral strains. Indeed, influenza is responsible for over 20000 deaths annually in the US, despite the annual seasonal vaccination campaigns that have become familiar [3,4]. It is particularly dangerous to the very young, elderly, and to those who have suppressed immune systems. Vaccination has been the main preventive measure and has provided limited control, but vaccines must be reformulated each year in response to antigenic drift and may be ineffective against new variants of influenza viruses [5].
There are three types of influenza viruses: A, B, and C. The first two types are responsible for annual epidemics and pandemic outbreaks; while influenza C is endermic vaccination provides limited protection based on the ability to predict the exact strains of influenza that will predominate a year in advance of the influenza season [4]. Millions of people each year become infected with the influenza A virus. Influenza A is characterized by the abrupt onset of constitutional and respiratory signs and symptoms (e.g., fever, myalgia, headache, severe malaise, nonproductive cough, sore throat, and rhinitis) [3]. Currently, control of influenza infection relies mainly on a preventive strategy dependent on the annual identification, production, and distribution of a multivalent vaccine designed to predict the predominant epidemic viral strains. Until recently, antiinfluenza drugs were restricted to the M2 channel blockers  amantadine and rimantadine, both of which are approved for the prophylactic and therapeutic treatment of the influenza viruses within the A family. However, the licensing of the neuraminidase inhibitors, zanamavir, and oseltamivir has provided the clinician with antiviral agents that effectively prevent and treat both type A and B influenza infections [1]. A new viral target has recently been exploited to develop clinically useful anti-influenza drugs [5].
All conformational and quantum-chemical data were obtained by means of the MM2 (molecular mechanic) method of the molecular mechanics and a semiempirical  quantum-chemistry method known as AM1. Activity features' selection has been carried out by means of the ETMsoftware. To have more stable activity features, every active molecule was used as a template for comparison with the rest of molecules. As a result of this comparison, activity features (pharmacophores) P1, P2, and P3 were revealed.
To decide which of pharmacophores is better, each inactive molecule was used as a template for comparison with the rest of molecules. So, inactivity features (antipharmacophores) AP1, AP2, and AP3 were revealed also to complete the system for the thiobenzamide and quinolizidine derivatives inhibitors activity prediction.

Materials and Methods
In this study, the structure anti-influenza relationships of thiobenzamide and quinolizidine derivatives being influenza fusion inhibitors have been studied by means of the ETM. An ETMC is an n×n matrix, where n is the number of atoms in the corresponding molecule. Diagonal elements are values of an atomic property such as atomic charge. An off-diagonal element can be of one of two kinds. It is bond property, if the bond exists, and distance for the corresponding pair of atoms, otherwise [19]. Detailed descriptions of the ETM can be found elsewhere [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20]. The matrices are easily understandable and extremely convenient for computer handling. In this way, the structures of the molecules under study get a unified description that is not bound to the atoms' identity. This circumstance is of primary importance when searching for pharmacophores in quite diverse structures of biologically active molecules.
The main steps of the ETM procedure are shown in Figure 2 [7,9,10,18]: The descriptors, such as E HOMO , E LUMO , ΔE, chemical potential, softness, electrophilicity index, dipole moment, and Mulliken charges of some atoms on studied molecules were calculated by the DFT method using B3LYP exchange correlation functionals [21,22] with the 6-31G(d,p) basis set by using Gaussian Program [23]. The Fukui functions

Results and Discussion
In this study, the optimized geometry data and electronic characteristics were used to form ETMCs for all molecules in a series of thiobenzamide and quinolizidine derivatives including 65 molecules [1,2]. The effective charge on atoms are shown on diagonal elements, while bond characteristics and optimized distances are represented on non diagonal elements.
The parameters responsible for the activity form a matrix called electron-topological submatrix of activity (ETSA), calculated from an ETMC that represents one of the most active molecules ("a template" for comparison).

Pharmacophores Identification.
For each template molecule, its ETMC was compared with the ETMCs of the rest of molecules in both series. The comparison resulted in a few common structural fragments for the two cases. The fragments were found as submatrices of the corresponding template ETMCs (i.e. electron-topological submatrices of contiguity (ETSCs).
Molecule 19 taken as template molecule from active ones was accepted to be the active one and from the template molecule 19, an activity feature 2 (or pharmacophore, P2) was found and presented in Figure 3, alone with the corresponding ETSCs, which describe electronic-topological characteristics of the fragments (see Figure 3).
In the matrices, the effective charges on atoms (local atomic characteristics, Q i ) were chosen for diagonal elements, while in place of nondiagonal elements representing interatomic characteristics there are either Wiberg's indices [21], (W i j , for bonds) or optimized distances (R i j , inÅ, for chemically non-bonded pairs of atoms). The distance between C 19 and O 22 atom is 8.80Å. Submatrices given in Figure 3 are found after setting some allowable limits for the matrix elements comparison. For both series, the limits are δ 1 = 0.04 for diagonal elements of the ETMCs and δ 2 = 0.14 for their off-diagonal elements comparison. The pharmacophores found from the ETM-calculations are realized in all active molecules (30). Statistical estimates for the pharmacophores are given in Table 2.
An activity feature P2 is shown along with its ETSA of the order 6x6 (see Figure 3). Only its upper triangle is given because of the symmetry of bounds. The pharmacophore was found in 21 of 30 active molecules, n A , and it was found in 3 of 28 inactive molecules, n IA . Thus, the probability P A of its realization in this class is about 0.84. As seen from the  pharmacophore structure, the active feature 2, P2, consists of the 6 atoms (C 3 , C 5 , H 11 , C 15 , C 19 , and O 22 ). Molecules 15 and 22 were chosen as templates and compared with the rest of the molecules in the series. Pharmacophore 1 (P1) was found in 20 of 30 active molecules having the probability P A = 0.84 and pharmacophore 3 (P3) were found in 21 of 30 active molecules having the probability P A = 0.85 as seen in Figure 3 and statistical estimates for the pharmacophores are given in Table 2.
Superimposition of three template compounds that correspond to the calculated pharmacophores P1-P3 is shown in Figure 4. In 3D space, three separate regions can be indicated for each active molecule, where atoms of three pharmacophores can be found.
These regions are shown in Figure 4 by dotted lines. One of them is formed by the atoms belonging to the phenyl ring, while the other two are formed by atoms that represent the heterocyclic rings. It is quite possible that the atoms from  these regions play an important role in the ligand-receptor interaction.

Antipharmacophores Calculation.
Antipharmacophores (alone with pharmacophores) are also of interest in this study, as those parts of molecules that are responsible for the considerable decrease or complete loss of the activity To find antipharmacophores, inactive molecules (39, 48, and 65) were selected as template molecules (their structures are given in Figure 5). Again, both protonated and unprotonated forms of molecules (39, 48 and 65) were studied. As an example, AP1, AP2, and AP3 antipharmacophores are shown in the figure by their numbers, while corresponding submatrices are given nearby.
The antipharmacophores consists of AP1 5 atoms, AP2 and AP3 4 atoms (see Figure 5). AP1 enters the structures of 15 inactive molecules and are found 2 active molecules (AP1, AP2, and AP3, see Table 2). The probability level of them are realized between 81% and 89%, because AP1, AP2, and AP3 are found in, respectively 15, 12, 15 from 28 inactive molecules in total. The distribution of the atoms are as shown in Figure 5.
The molecular orbitals (HOMO, LUMO) for template molecules taken as P1, P2, P3, AP1, AP2, and AP3 are shown in Figure 6, but we could not find any relationship between pharmacophore group and the molecular orbitals (HOMO and LUMO).
When comparing the structures of the pharmacophores and antipharmacophores, one can pay attention to the differences in their spatial and electron characteristics. Thus, pharmacophores and antipharmacophores can play their role in the activity prediction only if both types of fragments participate in the process of prognosis. Thus, the set of activity/inactivity fragments found as the result of this study forms a basis for a system of the analgesic activity prediction.
As seen from Figure 7, the pharmacophores P1-P3 and antipharmacophores AP1-AP3, found as a result of the ETM application, were used as a basis for a system formation that is capable of the prediction of the activity of the thiobenzamide and quinolizidine derivatives.

Combined ETM-ETM-ANNs Approach (ETM-ANNs)
The algorithm for the data resulting from the ETM calculations analysis (ETM-data) extends the ideas of volumelearning algorithm [26]. This method is also implemented as an iterative application of the Kohonen SOMs [27] and associative neural networks (ASNNs) [28]. The algorithm evaluates the "weights" of fragments represented by the ETMCs that have been obtained from the ETM calculations.
To do this, their projections on the Kohonen's maps corresponding to each ETMC are calculated and a weight of each fragment's presence in a molecule is determined [19,29]. The calculated weights are used for the ANN training. The original implementation of two last algorithms was used.

Selforganizing Map of Kohonen (SOM).
A SOM is neural network trained using unsupervised learning to create a nonlinear projection of high-dimensional input samples to a lower two-dimensional output space that is represented by a two-dimensional array of neurons, called a map [27]. The map seeks to preserve the topological properties of the input space. Like most artificial neural networks, SOMs operate in two modes: training and mapping. Training builds the map using input examples. It is a competitive process, also called vector quantization. During the SOM training, input vectors with similar properties in the high dimensional space are mapped to the same (or to the nearby) neurons on the two-dimensional space. Thus, by considering all input vectors projected to the same output neuron it is possible to determine clusters of vectors having similar properties in the high-dimensional space. Mapping automatically classifies a new input vector.  The avoidance of overfitting/overtraining has been shown to be an important factor for improvement of generalisation ability and correct selection of variables in ASNN. the current study used the early stopping over ensemble (ESE) technique to accomplish this. A detailed description of ESE can be found elsewhere [30]. In brief, each analysed artificial neural network ensemble (ANNE) was composed of M = 100 networks. The values calculated for analysed cases were averaged over all M neural networks, and their means were used for computing statistical coefficients with targets. We used a subdivision of the initial training set into two equal learning/validation subsets. The first set was used to train the neural network while the second one was used to control the training process. Two stopping points  were used to test network performance measured by root mean square error (RMSE). The first point (early stopping) determined a best fit of a network to the validation set while the second point corresponded to the error minimum for the learning set, and as a rule coincides with the end of the network training. The quality of each final model was assessed by the leave-one-out cross-validation method (LOO). By the method, each molecule was removed from the training set, and the remaining set was used to separate molecules into classes of activity, thereby predicting the activity of this molecule and evaluating the quality of the decision rule. The detailed description of the used methods can be found in [30,31].

Associative Neural Network (ASNN
Sensitivity analysis methods estimate the rate of change in the output of a model caused by the changes of the model inputs. It is mainly used to determine which input descriptor is more important or sensible to achieve accurate output values. It is also used to understand the behavior of the modeled system that to evaluate of the model applicability and to determine the stability of a model. To evaluate the importance of the original ETMC fragments, we have used sensitivity analysis methods named pruning algorithms [32,33]. The pruning algorithms introduce some measures of importance of weights matrix of ASNN by so called "sensitivities" (S). Since details of the sensitivity methods can be found in literature, we only give here the short description of the used method. The sensitivity of input neuron i was calculated as, where, S j was a sensitivity of the jth neuron in the upper layer and max a was taken over all weights ending at neuron j. This was a recurrent formula calculating the sensitivity of a unit on layer s via the neuron sensitivities on layer s+1. The sensitivities of output layer neurons were set to 1 [32]. The sensitivities of all neurons were calculated and the ETM descriptors associated with input neurons having smallest sensitivities was deleted.

Results of the ETM-ANN Approach Application.
A dataset containing 65 molecules was used. Fifty-two of the compounds were used for the model development, and 13 randomly selected compounds (5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, and 65) were used for the model validation. For the data, 210 fragments were selected. The importance of the detected fragments for the observed activity was evaluated by using pruning methods. The most part of the ETMC fragments (as descriptors) were detected as nonsignificant and removed by the pruning algorithms. As the result, only 13 ETMC-fragments were chosen from 210 fragments in total as the most important ones. By this, ASNNs classified correctly 92%, or 48 compounds from 52 in training set and 69%, or 9 compounds from external test set.

Model Interpretation.
The ETM-ASNNs models can be predominantly applied as virtual screening tools and also be useful for elucidating the essential structural and physicochemical requirements for activity. Usually, the main QSAR and pharmacophore elucidation approaches use three-dimensional steric and electrostatic potentials or chemical feature-based pharmacophoric ngerprints whereas the model reported in this paper uses 2D ngerprints (Mulliken atomic charges, Wiberg indices and bond length), in addition to several global physiochemical descriptors in the ASNN model development. In view of these differences, the current ETM-ASNN model reveals inherent features in terms of abstract substructural features (ETMS) that discriminate active and inactive molecules. This list of substructures preferentially present in actives and inactives (activity is specifically dependent on its conformational disposition on the cyclohexanol ring. In the trans position, it is dependent on the cyclohexanol group and the phenyl group and the substituents on the aromatic ring (especially the groups on C 3 and C 5 prevent the activity) can serve as a look-up guide to chemists during a hit-to-lead or lead optimization campaign. In general, substructures that have an incremental effect on activity (Figure 3) include the following: substituted (C 3 , C 5 , H11, C15, C 19 , and O 22 ) and substructures that have a decremental effect on activity include the following: substituted (C 3 , C 5 , H 14 , C 19 , and O 21 ).
The ASSN models based on descriptors set can be useful for understanding of the descriptor H1 influenza inhibitory activity of analysed molecules. The descriptors such as, chemical potential (μ), hardness (η), softness (s), electrophilicity index (ω), and condensed fukui functions, and quantum chemical descriptors such as highest occupied molecular orbital (HOMO) energy and lowest unoccupied molecular orbital (LUMO) energy, dipole moment, Mulliken charges are chemical reactivity descriptor [32][33][34][35][36]. For example, HOMO plays the role of the binding ability of the molecules. HOMO energy of molecules having high activity is lower than nonactive ones.

QSAR Study.
There are two main alternative approaches to random testing to find compounds with superior properties for correlation and quantitative prediction of chemical and physical properties from structure. One consists of theoretical calculations using quantum and statistical mechanics, and the other is QSAR/QSPR. The major goal of QSAR or QSPR studies is to find a mathematical relationship between activity or any other property under investigation (eg LD 50 , pKa etc.) and one or more descriptors related to the structure of the molecule. While some of the descriptors can themselves be experimental properties of the molecule, it is generally more useful to use descriptors derived mathematically from either 2D or 3D molecular structure, since this allows any relationship derived to be extended to the prediction of the property or activity for unavailable compounds. If an acceptable model of this type can be found, it can guide the synthetic chemist in the choice between alternative hypothetical structures. More fundamentally, such studies can illuminate, or even elucidate the mechanism by which the property or activity in question is related to the chemical structure.
The mathematical foundation of the QSAR is based on the principle of polylinearity. Multiple linear regression is a common method used in QSAR studies. The QSAR equations can be obtained by forward stepwise multiple regression shown below. Activity = f (physicochemical properties and/or structural properties) or f biological activity = f (electronic) + f (steric) or Biological activity = a 0 + a 1 D 1 + a 2 D 2 + a 3 D 3 + · · · a n D n , where D 1 , D 2 , D 3 and D n are descriptors, n is the number of descriptors. The intercept (a 0 ) and the regression coefficient of the descriptors are determined using the least squares method. In this study, by using the multiple regression, several QSAR models were performed with the help of some of the calculated descriptors and the compounds activity data from the different skeletons. Below are the regression equations. Skeleton From the three equations above, it shows that statistically the most significant one is the correlated parameters in Skeleton I with R 2 = 0.9999. There are some deviations for skeleton II and III. It also implies that the activity of the compounds in the skeletons does not only depend on the quantum chemical descriptors but may be influenced by other parameters.
Below are tables of theoretically calculated activities and experimental activities (Tables 3, 4, and 5). Skeleton I shows better correlation between the theoretically calculated activities and experimental activities while skeletons II and III do not show good agreement between them.
For each compound, 19 QSAR descriptors such as molecular weight, E HOMO , E LUMO , ΔE, chemical potential, softness, electrophilicity index, dipole moment, and others were calculated. Preliminary study by ASNNs for all descriptors resulted in model with good statistical qualities (classified correctly 87%, or 45 compounds from 52). Only 3 compounds in the external test set were predicted incorrectly (classified correctly 77%, or 10 compounds from 13).
Thus, the comparison of both models allowed us to conclude that both approaches provide practically similar results using the different set of descriptors. It should be noted that three identical molecules (5, 10, and 65) were classified incorrectly by both models. As it is seen, the approach presented in this paper has shown quite satisfactory results. This fact tells in favor of workability of both found models. These models can be applied to the design of new potent drugs.

Conclusion
A series of thiobenzamide and quinolizidine derivatives activity-binding affinity is studied by means of the ETM and ANN, which takes into account both structural and electronic characteristics of molecules. Based on pharmacophores and antipharmacophores calculated by the ETMsoftware (as submatrices containing important spatial and quantum chemistry characteristics), a system for the activity prognostication is developed. The system was tested on a few molecules with molecular skeletons other than those that were characteristic of the training sets. It allows for identifying the presence/absence of the analgesic activity (with probabilities 84%-85%) in molecules with diverse structures and predicting the level of the activity. QSAR models were performed with the help of some of the calculated descriptors and the compounds activity data from the different skeletons by using the multiple regression. Skeleton 1 shows better correlation between the theoretically calculated activities and experimental activities while skeletons 2 and 3 do not show good agreement between them.
QSAR descriptors such as molecular weight, as E HOMO , E LUMO , ΔE, chemical potential, softness, electrophilicity index, dipole moment, and others were calculated to compare classical and ETM model. Thus, the comparison of both models allowed us to conclude that both approaches provide practically similar results using the different set of descriptors. It should be noted that three identical molecules (5, 10, 65) were classified incorrectly by both models. As it is seen, the approach presented in this paper has shown quite satisfactory results. This fact tells in favor of workability of both found models. These models can be applied to the design of new potent drugs.