Modeling of acetylcholinesterase inhibition by tacrine analogues using Bayesian-regularized Genetic Neural Networks and ensemble averaging

Abstract

Acetylcholinesterase inhibition was modeled for a set of 136 tacrine analogues using Bayesian-regularized Genetic Neural Networks (BRGNNs). In the BRGNN approach the Bayesian-regularization avoids overtraining/overfitting and the genetic algorithm (GA) allows exploring an ample pool of 3D-descriptors. The predictive capacity of our selected model was evaluated by averaging multiple validation sets generated as members of diverse-training set neural network ensembles (NNEs). The ensemble averaging provides reliable statistics. When 40 members are assembled, the NNE provides a reliable measure of training and test set R values of 0.921 and 0.851 respectively. In other respects, the ability of the nonlinear selected GA space for differentiating the data was evidenced when the total data set was well distributed in a Kohonen Self-Organizing Map (SOM). The location of the inhibitors in the map facilitates the analysis of the connection between compounds and serves as a useful tool for qualitative predictions.

Keywords::

• Bayesian regularized-Genetic Neural Networks
• neural network ensemble
• self-organizing maps
• acetylcholinesterase inhibitors
• QSAR

Introduction

Tacrine (9-amino-1,2,3,4-tetrahydroacridine) was the first drug for the symptomatic treatment of Alzheimer's disease (AD) [1]. The rationale for its use is related to the elevation of the acetylcholine (ACh) levels by reversible inhibition of acetylcholinesterase (AChE), which can compensate for the cholinergic deficit associated to the brain lesions present in AD. However, there are several deficiencies of tacrine as a drug which are related to liver toxicity and peripheral cholinomimetic actions [2]. In order to reduce these undesirable side effects, many analogues of tacrine have been reported 3-11, most of these are structurally closely related to the parent compound and retain the aminopyridine or aminoquinoline moiety (Figure 1).

Figure 1 Structure of the tacrine analogues from modification of tacrine.

Computational-based rational design of drugs has increased in the last decade. Most of those approaches are focused on quantitative structure–activity relationship (QSAR) studies, using different kinds of molecular descriptors for encoding chemical information [12]. After computing a set of descriptors, multivariate linear or/and nonlinear relationships are established between a reduced subset of variables and the inhibitory activity, leading to a mathematical model. AChE inhibitors have been approached by many kinds of computational strategies. A serious source of information about this topic is offered by Dimoglo et al. [13]. In this paper, the authors made a comprehensive review of the ‘state of the art’ concerning AChE inhibition modeling.

Artificial neural networks (ANNs) arose from attempts to model the functioning of the human brain [14]. They usually overcome methods limited to linear regression models such as multiple linear analysis or partial least square 15-18. Contrary to these methods, ANNs can be used to model complex nonlinear relationships. Since biological phenomena are complex by nature, this ability has promoted the employment of ANNs in drug design studies. In the current study we model the AChE inhibitory activity of a set of tacrine analogues 3-11. The data set includes 136 compounds with the biological activity reported as IC₅₀ values. The characteristics of the inhibitors were represented by relevant 3D descriptors extracted by genetic algorithm (GA) feature selection using Bayesian-regularized ANNs (BRANNs) as predictors. The model was carried out by a neural network ensemble approach which provides reliable statistics. In addition, the capacity of the selected variables for differencing the data was evaluated by means of the unsupervised training of Kohonen Self-Organizing Maps (SOMs).

Materials and methods

Dataset and molecular descriptors

AChE inhibitory activities [log(10⁵/IC₅₀)] of 136 tacrine analogues were taken from the literature 3-11. The chemical structures and experimental activities are shown in Table I. IC₅₀ values represent the micromolar concentration that inhibits 50% of AChE activity. For inactive compounds the threshold value of log(10⁵/IC₅₀) = 3 (IC₅₀ = 100 μM) was assigned. Since the whole dataset was extracted from reports of several groups, there are some differences between the AChE sources for pharmacological assays. In order to evaluate the outcomes of these differences, we analyzed the correlation between inhibitory activities using bovine and human AChEs for 17 compounds reported by Camps et al. [7,8]. Table II. shows that bovine log(10⁵/IC₅₀) and human log(10⁵/IC₅₀) values correlates for these compounds (R² > 0.8), which was improved when one unfitted compound was extracted. In the equation that describes this relationship, the slope tends to 1 and the intercept tends to zero. Furthermore, the equation fixing the intercept also adjusts well this relationship. According to this analysis, we consider that the differences between AChE sources do not cause a high error to our work.

Table I.  Experimental and predicted AChE inhibitory activities (log(10⁵/IC₅₀)) of tacrine analogues^a.

			log(10^5/IC50)^c
comp		exp	training	test
1 (tacrine)		5.60	5.40	5.33
							log(10^5/IC50)^c
comp	R	R1	R2	R3	R4	R5	exp	training	test
2	H	H	H	H	H	H	5.17	4.63	4.55
3	H	NO2	H	H	H	H	< 3.00	3.60	3.71
4	H	NH2	H	H	H	H	4.23	4.41	4.41
5	H	H	NO2	H	H	H	3.17	4.63	5.02
6	H	H	NH2	H	H	H	3.54	4.28	4.40
7	H	H	Cl	H	H	H	4.19	4.78	4.90
8	H	H	F	H	H	H	4.92	4.54	4.48
9	H	H	H	OMe	H	H	4.80	4.69	4.16
10	H	H	H	H	OMe	H	4.19	3.72	3.42
11	H	H	H	H	H	OMe	4.37	4.06	3.94
12	H	H	H	Me	H	H	4.41	4.53	4.52
13	H	H	H	H	H	Me	4.34	4.37	4.42
14	H	H	H	H	F	H	5.37	4.51	4.34
15	H	H	H	H	Cl	H	4.27	4.16	4.09
16	Bn	H	H	H	H	H	4.15	4.31	4.33
17	Me	H	H	H	H	H	4.89	4.90	4.96
18	n-C7H15	H	H	H	H	H	4.37	4.72	4.83
				log(10^5/IC50)^c
comp	R	R1	R2	exp	training	test
19	H	Me	H	4.09	5.07	5.23
20	H	H	Me	6.00	5.19	4.92
21	H	Cl	H	5.26	5.31	5.34
22	H	H	Cl	7.00	6.59	6.44
23	H	NO2	H	4.52	5.04	5.22
24	H	H	NO2	6.55	5.43	5.21
25	H	H	OMe	5.46	5.16	5.07
26	H	NH2	H	4.42	5.38	5.54
27	H	H	F	6.06	5.60	5.56
28	H	Cl	Cl	5.33	5.76	6.12
29	H	OMe	OMe	4.28	4.36	4.42
30	Bn	Me	H	4.43	4.39	4.39
31	Bn	H	Me	5.12	4.61	4.54
32	Bn	H	Cl	5.77	4.68	4.54
33	Bn	NO2	H	4.80	4.57	4.51
34	Bn	H	NO2	4.32	4.64	4.73
35	n-C7H15	Me	H	5.41	5.63	5.68
36	n-C7H15	H	Me	5.89	5.55	5.51
37	n-C7H15	H	Cl	6.89	6.51	6.35
38	n-C7H15	H	NO2	5.54	5.82	5.87
39	n-C7H15	H	OMe	5.34	5.67	5.80
40	n-C7H15	H	F	6.35	6.09	6.02
41	H	H	Br	6.18	5.83	5.72
						log(10^5/IC50)^c
comp^b	X(C(Xanti)Xsyn)	R	R1	R2	R3	exp	training	test
42	CH2	H	H	H	H	4.31	5.04	5.11
43	C(Me)OMe	Me	H	H	H	4.87	5.05	5.13
44	C(OMe)Me	Me	H	H	H	4.17	4.53	4.62
45	CO	Me	H	H	H	4.68	5.56	5.71
46	(E)–C ( = CHMe)	Me	H	H	H	4.94	5.69	6.17
47	(Z)–C ( = CHMe)	Me	H	H	H	5.49	5.30	5.29
48	o-phenylene	Me	H	H	H	4.68	4.83	4.91
49	CH2	n-Pr	H	H	H	5.37	5.84	5.95
50	CH2	i-Pr	H	H	H	5.99	5.73	5.53
51	CH2	CH2CH = CH2	H	H	H	5.82	6.36	6.48
52	CH2	n-Bu	H	H	H	5.55	5.93	5.97
53	CH2	t-Bu	H	H	H	5.57	5.61	5.46
54	CH2	Ph	H	H	H	5.90	5.08	4.80
55	CH2	Me	H	H	H	6.19	5.80	5.77
56	CH2	Me	H	H	Me	6.91	6.14	5.93
57	CH2	Me	F	H	H	6.50	6.49	6.52
58	CH2	Me	H	H	F	7.07	6.39	6.25
59	CH2	Me	F	H	F	7.61	6.94	6.78
60	CH2	Me	H	H	Cl	7.37	7.13	7.10
61	CH2	Et	H	H	H	6.41	6.05	6.02
62	CH2	Et	Me	H	H	6.53	6.27	6.07
63	CH2	Et	H	H	Me	6.92	6.55	6.43
64	CH2	Et	Me	H	Me	7.44	7.30	7.30
65	CH2	Et	F	H	H	6.33	6.52	6.59
66	CH2	Et	H	H	F	7.13	6.76	6.68
67	CH2	Et	F	H	F	7.58	7.18	7.07
68	CH2	Et	Cl	H	H	6.79	7.17	7.34
69	CH2	Et	H	Cl	H	5.59	6.58	6.88
70	CH2	Et	H	H	Cl	7.56	7.33	7.27
71	CH2	Et	Cl	H	Cl	6.40	6.90	8.00
		log(10^5/IC50)^c
comp^b	R	exp	training	test
72	Me	5.29	5.40	5.39
73	Et	5.78	5.28	5.22
				log(10^5/IC50)^c
comp^b	X(C(Ranti)Rsyn)	R(Rexo,Rendo)	exp	training	test
74	CH2	H,H	3.38	4.46	4.65
75	CH2	OH,H	4.79	4.45	4.40
76	o-phenylene	H,H	4.59	4.30	4.21
77	CH2	H,OH	4.37	4.47	4.49
78	o-phenylene	H,OH	3.82	4.22	4.61
				log(10^5/IC50)^c
comp^b	X(C(Ranti)Rsyn)	R(Rexo,Rendo)	exp	training	test
79	CH2	Me	4.25	5.02	6.02
80	o-phenylene	O	4.53	3.70	3.39
81	CH2	H,OH	3.87	3.62	3.66
				log(10^5/IC50)^c
comp^b	R	R1	n	exp	training	test
82	H	Et	1	4.81	4.50	4.45
83	4-Me	Et	1	4.74	4.32	4.29
84	4-Cl	Et	1	4.71	4.78	4.96
85	4-CN	Et	1	4.46	4.58	4.58
86	4-OMe	Et	1	5.06	4.23	4.16
87	3-NO2	Et	1	4.72	4.58	4.54
88	H	Et	0	4.41	4.22	4.26
89	3-OMe	Et	0	< 3.00	3.80	3.98
90	3-OMe	Et	1	3.52	4.07	4.12
91	3-OMe	Et	2	3.89	3.83	3.82
92	2-OMe	Et	0	3.62	3.81	3.90
93	2-OMe	Et	1	4.11	4.01	4.00
94	2-OMe	Et	2	3.68	4.05	4.11
95	3,4-diOMe	Et	0	3.57	3.91	3.56
96	3,4-diOMe	Et	1	4.34	4.13	4.10
97	3,4-diOMe	Et	2	3.48	3.85	3.89
98	H	i-Pr	1	3.52	3.87	3.91
				log(10^5/IC50)^c
comp	R	R1	n	exp	training	test
99	H	Et	0	3.92	4.20	4.27
100	H	Et	1	5.09	4.56	4.49
101	H	Et	2	3.68	4.46	4.81
102	4-OMe	Et	0	3.86	3.87	3.91
103	4-OMe	Et	1	4.05	4.20	4.24
104	4-OMe	Et	2	4.87	3.99	3.94
105	4-OMe	Et	1	4.92	4.75	4.64
106	3,4-diOMe	Et	1	4.89	4.77	4.74
107	4-F	Et	1	5.15	4.69	4.64
108	4-Cl	Et	1	4.89	4.67	4.63
109	3-NO2	Et	1	5.03	4.97	4.96
110	4-Me	Et	2	3.89	4.20	4.30
111	3-OMe	Et	2	3.52	4.40	4.51
112	2-OMe	Et	2	3.66	3.96	3.97
113	3,4-diOMe	Et	2	4.19	4.00	3.95
114	4-F	Et	2	4.27	4.11	4.08
115	4-Cl	Et	2	3.72	4.18	4.24
116	3-NO2	Et	2	3.85	4.01	4.05
117	H	Et	3	3.77	3.70	3.66
118	4-OMe	Et	3	3.77	3.48	3.39
119	H	i-Pr	1	4.92	4.63	4.59
				log(10^5/IC50)^c
comp	X	R	n	exp	training	test
120	O	H	1	5.36	4.86	4.63
121	O	Me	1	5.42	5.24	5.23
122	O	OMe	1	5.34	5.43	5.51
123	O	H	2	5.50	5.27	5.19
124	S	H	1	< 3.00	2.75	2.14
125	S	H	2	< 3.00	3.67	4.59
			log(10^5/IC50)^c
comp	R	N	exp	training	test
126	Cl	0	< 3.00	3.12	3.15
127	OMe	0	< 3.00	3.12	3.18
128	Me	0	< 3.00	3.21	3.33
129	Cl	1	< 3.00	3.14	3.27
130	NO2	1	< 3.00	3.03	2.98
131	OMe	1	< 3.00	3.11	3.08
132	Me	1	< 3.00	2.68	2.58
133	Cl	2	< 3.00	2.90	2.82
134	NO2	2	< 3.00	2.99	2.98
135	OMe	2	< 3.00	2.99	3.02
136	Me	2	< 3.00	2.99	2.97

^a Compounds 2–18 are from Ref. 3; 19–41 are from Ref. 4; 42–48, 55, 61 and 74–81 are from Ref. 5; 49–54, 56–58 and 62–63 are from Ref. 6; 59–60 and 64–71 are from Ref. 7; 72–73 are from Ref. 8; 82–88 and 99–104 are from Ref. 9; 120–125 are from Ref. 10 and 89–98, 105–119 and 126–136 are from Ref. 11.

^b racemic mixtures.

^c IC₅₀: micromolar concentration that inhibits 50% of AchE activity.

Table II.  Correlation between log(10⁵/IC₅₀) using bovine AChE and log(10⁵/IC₅₀) using human AChE^a.

Features of the correlation.	n	slope	intercept	R^2
Adjusting intercept	17	1.108	− 0.497	0.814
Fixing intercept	17	1.034	0.000	0.810
Adjusting intercept	16^b	1.145	− 0.818	0.934
Fixing intercept	16^b	1.023	0.000	0.924

^a Equation: bovine-log(10⁵/IC₅₀) = slope × human-log(10⁵/IC₅₀)+intercept.

^b An unfitted compound was eliminated.

Prior to molecular descriptor calculations, 3D structures of the studied compounds were geometrically optimized using the semiempirical quantum-chemical method PM3 [19] implemented in the MOPAC 6.0 [20] computer software. The 3D descriptors from the Dragon software [21] were calculated for each compound: aromaticity indices [22,23], Randic molecular profiles [24], geometrical descriptors [25], RDF (radial distribution function) descriptors [26], 3D-MoRSE (molecule representation of structures based on electron diffraction) descriptors [27], WHIM (weighted holistic invariant molecular) descriptors [28] and GETAWAY (GEometry, Topology and Atom-Weight AssemblY) descriptors [29]. In all, 721 descriptors were calculated. Descriptors that stayed constant or almost constant were eliminated, and pairs of variables with a correlation coefficient greater than 0.95 were classified as intercorrelated, and only one of these was included in the model. Finally, 271 descriptors were obtained.

Artificial neural network regression procedure

Only a subset of the available pool of 271 descriptors is statistically significant in terms of correlation with the modeled AChE inhibitory activity. Our general procedure is described in Figure 2.

Figure 2 Flowchart of the modeling procedure.

Bayesian-Regularized Genetic Neural Network (BRGNN) is a framework that combines BRANNs with GA feature selection [30]. Our BRGNN approach is a version of the So and Karplus GA feature selection method [31] incorporating Bayesian regularization.

Bayesian networks are optimal devices for solving learning problems. They diminish the inherent complexity of ANNs, being governed by Occam's Razor, when complex models are automatically self-penalizing under Bayes's rule. The Bayesian approach to ANN modeling considers all possible values of network parameters weighted by the probability of each set of weights. The BRANN method was designed by Mackay [32,33] for overcoming the deficiencies of ANNs. Only a brief summary will be provided here. The Bayesian approach yields a posterior distribution of network parameters P(w|D,H) from a prior probability distribution P(w|H) according to updates provided by the training set D using the BRANN model H. Predictions are expressed in terms of expectations with respect to this posterior distribution. Bayesian methods can simultaneously optimize the regularization constants in ANNs, a process that is very laborious using cross validation. Instead of trying to find the global minimum, the Bayesian approach finds the (locally) most probable parameters.

Bayesian approach produces predictors that are robust and well matched to the data. These properties become BRANNs in accurate predictors for QSAR analysis [34,35]. They give models which are relatively independent of ANN architecture, above a minimum architecture, since the Bayesian regularization method estimates the number of effective parameters. The concerns about overfitting and overtraining are also eliminated by this method so that the production of a definitive and reproducible model is attained. The joining of BRANN and GA feature selection (BRGNN) increases the possibilities of BRANNs for modeling as we indicated in previous works [16,17,30]. This method is relatively fast and considers the whole data set in the training process. For other hybrids of ANN and GA the use of the MSE as fitness function could lead to undesirable well fitted but poor generalized networks as algorithm solutions. In this connection, BRGNN avoids such results by two aspects: 1) keeping network architectures as simple as possible inside the GA framework and 2) implementing Bayesian regulation in the network training function.

Fully connected, three-layer BRANNs with back-propagation training were implemented in the MATLAB environment [36]. In these nets, the transfer functions of input and output layers were linear and the hidden layer had neurons with a hyperbolic tangent transfer function. Inputs and targets took the values from independent variables selected by the GA and log(10⁵/IC₅₀) values, respectively; both were normalized prior to network training. BRANN training was carried out according to the Levenberg-Marquardt optimization [37]. The initial value for μ was 0.005 with decrease and increase factors of 0.1 and 10, respectively. The training was stopped when μ became larger than 10¹⁰.

The GA implemented in this paper keeps the same characteristics of the previously reported in earlier work [30]. Initially, a set of 50 chromosomes were randomly generated. The population fitness was then calculated and the members were rank ordered according to fitness. The 2 best scoring models were automatically retained as members for the next round of evolution. More progeny models were then created for the next generation by preferentially mating parent models with higher scores. Crossover operator and single-point mutations were used in the evolution process until the best MSE scoring model remains constant for at least 10 generations. Our GA was programmed within the MATLAB environment using the genetic algorithm and neural networks toolboxes [36]. The predictors are BRANNs with a simple architecture (two or three neurons in a sole hidden layer). We tried the MSE of data fitting for BRANN models, as the case may be, as the individual fitness function. The best models according to R value (R > 0.8) were selected, and they were tested in cross-validation experiments for avoiding chance correlations (Figure 2).

The predictive power of the model was measured by an external validation process that consists of predicting the activity of unknown compounds forming the test set. To avoid the influence of casual external sets, neural network ensembles (NNEs) were employed; building all members by the random partition of the whole data set into training (80%) and test sets (20%) following the Agrafiotis et al. proposition [38]. As a result, averaging external predictions were obtained (Figure 2). The quality and reliability was settled by examining the correlation coefficient R and the MSE of the test set fitting.

Assembling multiple versions of a predictor provides ‘smoother’ more stable predictions [39]. Ensemble averaging minimizes uncertainty and produces more stable and accurate predictors. Recently, our group demonstrated the advantages of the ensemble solution for QSAR validation [40]. The robustness of this method lies on the adequacy of many external predictions; therefore, it can replace the traditional validation processes.

Kohonen self-organizing maps

The Kohonen SOMs [41] are ANNs related to classic clustering algorithms, in that they generate groupings of data points taken to be described by a single vector of typical values. However, the SOMs are distinct from standard clustering methods in that they do not operate with separate clusters: rather, they allocate data points to groups which are related [42].

Kohonen SOMs are networks of spatially related nodes each of which represents a ‘prototype’ of a particular region of data (input) space. Each node comprises a set of weights corresponding to the dimensions of the data. Their characteristic feature is their ability to map nonlinear relations in multidimensional data sets into easily visualizable two-dimensional grids of neurons displaying the topology of a data set. Essentially, SOMs permit the perception of similarities in objects.

To settle structural similarities among the modeled AChE inhibitors, a Kohonen SOM was built. The 3D descriptors selected by GA were used for unsupervised training of 13 × 13 neuron maps. SOMs were implemented in the MATLAB environment [36]. Neurons were initially located in a grid topology. The ordering phase was developed in 1000 steps with a 0.9 learning rate until a tuning neighborhood distance (1.0) was achieved. The tuning-phase learning rate was 0.02. Training was performed for a period of 2000 epochs in an unsupervised manner.

Results and discussion

Molecular descriptors

BRGNN methodology was applied for evaluating nonlinear relationships between the inhibitory activity and molecular descriptors. We found a seven descriptor space able to explain the AChE inhibition; the descriptors that constitute this space are shown in Table III. The selected space includes the geometrical descriptor SPAM, the RDF descriptor RDF075m, four 3D-MoRSE descriptors and the GETAWAY descriptor HATS6u. It is noteworthy that there is no significant intercorrelation between these descriptors, as is seen in Table IV.

Table III.  Symbols of the 3D-descriptors selected by Genetic Algorithm and their definitions.

Variable	Definition
SPAM	Average span R.
RDF075m	Radial Distribution Function at 7.5Å weighted by atomic masses.
Mor32u	3D MoRSE signal 32 unweighted.
Mor03v	3D MoRSE signal 3 weighted by atomic van der Waals volumes.
Mor17v	3D MoRSE signal 17 weighted by atomic van der Waals volumes.
Mor32p	3D MoRSE signal 32 weighted by atomic polarizabilities.
HATS6u	Leverage-weighted autocorrelation of lag 6 unweighted.

Table IV.  Correlation matrix of the descriptors selected by Genetic Algorithm.

	SPAM	RDF075m	Mor32u	Mor03v	Mor17v	Mor32p	HATS6u
SPAM	1
RDF075m	0.368	1
Mor32u	0.098	0.023	1
Mor03v	0.054	0.280	0.044	1
Mor17v	0.553	0.432	0.093	0.232	1
Mor32p	0.000	0.004	0.582	0.106	0.010	1
HATS6u	0.080	0.224	0.203	0.000	0.252	0.196	1

Geometrical descriptor SPAM [43] describes the 3D structures of tacrine analogues based on the Volkenstein approach for estimation of the statistical properties of polymer chains considering short-range interactions.

RDF descriptors [26] are calculated from the radial distribution function of an ensemble of N atoms that can be interpreted as the probability distribution of finding an atom in a spherical volume of radius r. Equation (1) represents the radial distribution function code: where f is a scaling factor, N is the number of atoms, A_i and A_j are atomic properties of atoms i and j, r_ij represents the interatomic distances and B is an smoothing parameter, that defines the probability distribution of the individual distances. g(r) was calculated at a number of discrete points with defined intervals. Different atomic properties A_i were used, such as atomic mass, atomic van der Waals volumes, atomic Sanderson electronegativities, and atomic polarizabilities. The possibility for choosing an appropriate atomic property gives great flexibility to the RDF space for adapting it to the problem under investigation. RDF075m takes into account the atoms around 7.5Å in the atomic mass weighting scheme.

3D-MoRSE [27] code considers the molecular information derived from an equation used in electron diffraction studies. Electron diffraction does not directly yield atomic coordinates but provides diffraction patterns from which the atomic coordinates are derived by mathematical transformations. 3D-MoRSE code is applied by Equation (2):

In this equation, A_i and A_j are atomic properties of atoms i and j, r_ij represents the interatomic distances, and s measures the scattering angle. The value of s (0,…, 31.0Å^{− 1}) is considered only at discrete positions within a certain range. Values of I(s) are defined at 32 evenly distributed values of s in the range of 0–31.0Å^{− 1}. These 32 values constitute the 3D-MoRSE code of the three-dimensional structure of a molecule. Like in the RDF approach, atomic properties A_i were used (atomic masses, atomic van der Waals volumes, atomic Sanderson electronegativities and atomic polarizabilities). The possibility for choosing an appropriate atomic property gives great flexibility to the 3D-MoRSE code for adapting it to the problem under investigation. In this work, 3D-MoRSE-selected descriptors were the unweighted Mor32u, weighted by van der Waals volumes Mor03v and Mor17v and weighted by atomic polarizabilities Mor32p. This code represents a restricted 3D space which captures relevant molecular information, which is related to the modeled AChE inhibitory activity.

GETAWAY descriptors [29] are based on the molecular influence matrix (MIM). They match the 3D molecular geometry provided by the MIM and atom relatedness by molecular topology, with chemical information. The diagonal elements h_ii of the MIM, called leverages, encode atomic information and represent the ‘‘influence’’ of each molecule atom in determining the whole shape of the molecule; in fact, mantle atoms always have higher h_ii values than atoms near the molecule center. Each off-diagonal element h_ij represents the degree of accessibility of the jth atom to interactions with the ith atom. Specifically, HATS6u comes from H-GETAWAY descriptors based on Moreau-Brotto autocorrelation descriptors [44]. In such descriptors, geometrical information provided by leverage values is combined with atomic weightings, accounting for specific physicochemical properties of molecule atoms. HATS indices consider the MIM diagonal elements. Like in the Moreau-Broto autocorrelations, HATSk(w) descriptors are defined as:

Where k (1, 2, …, d) is the path length (lag) in the molecular graph, d_ij is the topological distance between atoms i and j, while w_i and w_j are the A-dimensional property vector of the atoms i and j. The function δ(k; d_ij) is a Dirac-delta function defined as

The information contained in RDF, 3D-MoRSE and GETAWAY selected descriptors is related to distributions of relevant atomic properties across the entire molecules. In this sense, the elucidation of the molecular key features of the modeled data set from our proposed space is a difficult task. However, it can be perceived according to the nonlinear structural information here obtained, that an adequate distribution of atomic masses, van der Waals atomic volumes and polarizabilities has a great influence on the AChE inhibitory activities of the tacrine analogues. It suggests that molecular size, shape and atomic constitution, instead of electronic properties, play an important role in the tacrine analogues' activity. These facts agree well with reports where the access of the inhibitors to the AChE buried active site are strongly limited by the presence of a long and narrow gorge or channel leading from the surface of the enzyme [45]; this aspect imposes shape and size-related constraints for AChE inhibitors. A striking feature of this gorge is that it is lined by 14 aromatic residues which make up approximately 40% of its surface and allow the AChE active site to be qualified as hydrophobic.

In order to gain a deeper inside on the relative effects of each 3D descriptor in our model, a recently reported weight-based input ranking scheme was carried out. The black-box nature of three layers ANNs has been “deciphered” in a recent report by Guha et al. [46]. Their method allows an understanding of how an input descriptor is correlated to the predicted output by the network and consists of two parts. First, the nonlinear transform for a given neuron is linearized. Afterward, the magnitude in which a given neuron affects the downstream output is determined. Next, a ranking scheme for neurons in the hidden layer is developed. The ranking scheme is carried out by determining the Square Contribution Values (SCV) for each hidden neuron (see Reference 46 for details). This method for ANN model interpretation is similar in manner to the partial least squares interpretation method for linear models described by Stanton [47].

Results of the ANN deciphering study appear in Table V. The reported effective weight matrix for our model shows that the third hidden neuron has the major contribution to the model with a SCV value 4-fold higher in comparison to the first hidden neuron. On this neuron, SPAM and Mor03v descriptors have the highest impacts equal to − 2.774 and 2.676, respectively. From this analysis we can also derive the approximate effect of the selected descriptors. The sign of the weights indicates the trend of the output value. According to the sign of the effective weights in each hidden neuron, neither descriptor shows a completely positive or negative trend which suggests complex nonlinear effects. The most relevant descriptors SPAM and Mor03v, which have opposed effects in the third hidden neuron (negative and positive respectively), change their effects in other neurons.

Table V.  Effective weight matrix for the optimum model for the AChE inhibitory activity of the studied tacrine analogues^a. Most relevant descriptors appear in bold letter.

	Hidden neurons
Network inputs	3	1	2
SPAM	− 2.774	1.119	2.850
RDF075m	0.846	− 0.037	0.021
Mor32u	− 0.071	0.549	− 1.384
Mor03v	2.676	− 0.681	− 1.412
Mor17v	− 0.801	0.371	− 0.061
Mor32p	− 0.652	− 0.170	2.472
HATS6u	0.924	− 0.522	0.297
SCV	0.777	0.204	0.019

^a The columns are ordered by the SCVs for the hidden neurons, shown in the last row

Ensemble averaging

The above selected descriptors constitute the vectorial space that describes the modeled activity in a better way. For validating our vectorial space we evaluated the correlation coefficient and standard deviation (Q² and S_CV values) in LOO cross-validation process. The network hidden layer's architecture was explored, but no relevant differences were found due to the Bayesian-regularization. Q² reached the value of 0.745 and S_CV = 0.611 when the 7-3-1 architecture was employed. As can be seen, the correlation and predictive power behaviours were satisfactory so that the ability of the GA-selected descriptors to act as relevant information for the modeled activity is also confirmed in this case.

In order to obtain a more reliable model, we evaluated the predictive power by creating several training sets and predicting activities of unknown compounds. Instead of the selection of a sole test set, we generated multiple ones by means of the NNE approach. Then, we averaged external predictions. Members of NNEs were randomly generated by dividing the whole data set into 109 inhibitors for the training sets (80%) and 27 inhibitors for the test sets (20%), keeping the previous BRANN architecture. Averaged multiple correlation coefficients (averaged R) of the test sets for 20 instances of NNEs with 1, 5, 10, 20 and 40 predictors were examined (Figure 3). All averaged R values stem from adding up 1, 5, 10, 20 and 40 external set R values containing 27 inhibitors each.

Figure 3 Multiple correlation coefficients (R) of 20 replies of neural network ensembles using 1, 5, 10, 20 and 40 members for test sets.

NNEs containing one member are cases when single training and test sets were selected without integrating them to any ensemble. It is obvious that random partitioning is highly unsatisfactory. Diverse partitions generate a broad scope of external sets, and some are far better predicted than others. Therefore, even when it is broadly employed, to assess the predictive power by means of a single external data set, random selection yields a rather fortuitous result. It is noteworthy that more reliable information can be acquired when the number of members in the ensemble is increased. Ensembles containing 40 members are similarly predictive according to averaged test set R values around 0.841. The accumulation of members leads to an averaged model that weights the contribution of each predictor; in this form, deceptive conclusions are suppressed. Since our test set consists of 20% of the whole data set, the probability that an inhibitor will be selected as part of a test set is low. When a new partition is carried out, the new test set can contain elements from the original test set and new elements. The assemblage of successive members allows storage of predictions for each compound. Whereas the higher is the number of members, most replies of predictions can be collected; therefore, we can establish a test set prediction for the ensemble by averaging them, which includes all compounds when the number of members is sufficiently large. The NNE provides a reliable measure of training and test set R values. According to this, we report training set R = 0.921 and test set R = 0.851 when 40 members are assembled. The plots of training and external predictions versus experimental activities employing 40-member NNEs are shown in Figure 4.

Figure 4 Plot of predicted versus experimental log(10⁵/IC₅₀) values for AChE inhibition by tacrine analogues using 40-member NNEs. () training predictions; (○) external predictions.

Since the ensemble approach leads to stable ANN models, it becomes easeful to inspect the influence of some aspects that may affect the quality of our modeling procedure. The number of hidden nodes defines the number of parameters, which must be carefully adjusted in traditional ANNs. This is less critical in BRANNs, however it is useful to study the statistical behaviour of our ensemble when the number of hidden nodes are increased. Figure 5 shows the consequence of varying the number of hidden nodes in test set R value for 40-member NNE. We conclude that there is not an appreciable effect due to the Bayesian regularization. Otherwise, we explored other partitions instead of 20% for generating test sets. This analysis allows evaluating which proportion of the dataset is foreseeable by the rest, which means the redundancy in the whole dataset. Figure 6 shows the consequence of varying dataset partition in test set R value for 40-member NNE. The statistics of test set fitting declines when the test set elements are above 50% which indicates that one half of the data contains the information of all the modeled tacrine analogues.

Figure 5 Plot of test set R values vs. number of hidden nodes in neural networks for 40-member NNE.

Figure 6 Plot of test set R values vs. % of compounds in test set in neural networks for 40-member NNE.

Kohonen self-organizing neural network analysis

In order to achieve data differentiation, a Kohonen SOM with 13 × 13 neurons was mapped with GA selected descriptors as input vectors. In a self-organizing neural network, if two input data vectors are similar, they will be mapped into the same neuron or into very close neurons in the two-dimensional map. Therefore, either group in the map can be interpreted as a set of analogues defined by the vectorial space.

Figure 7 depicts a Kohonen SOM for the 136 AChE modeled inhibitors. It is clearly seen that compounds are adequately distributed across the entire map: 77 out of a total of 169 neurons were occupied. As observed, compounds with a similar range of activities were grouped into neighboring areas. It is noteworthy that there is a kind of gradient from the less-active compounds across the upper zone to the most-active compounds at the lower-right zone. As a consequence, this map can be used to carry out qualitative predictions. The position in the map would be able to assign an approximate range of activity for unknown compounds.

Figure 7 Kohonen SOM for the data set using descriptors selected by genetic algorithm. Maps at bottom represent the ranges of AChE activities (log10⁵/IC₅₀).

The analysis of test set predictions of the 136 modeled AChE inhibitors reported in Table I, suggests that compounds 5, 24, 32, 46, 69, 71, 74, 79 and 125 can be considered as test set outliers (residuals > 2 S_CV) in our vectorial space. The reason for the wrong prediction of these compounds is related to inaccurate associations with similar molecules containing unequal activity. Some evidence can be derived from a closer inspection of the Kohonen SOM of Figure 7. Compounds 5, 46, 69, 71, 74, 79 and 125 were predicted as more active than they really are. This fact can be attributed to the high molecular similarity that these compounds share with more-active compounds, thus, avoiding a correct discrimination; it is hinted from the position they occupy in the map. Compounds with a range of log(10⁶/IC₅₀) values between 3 and 3.99 are in the upper-left, but compounds 5 and 74 are out of this zone. Meanwhile, the inactive compound 125 is located far from the other inactive ones. Besides, compounds 46, 69, 71 and 79 were located in the lower-right zone, related to the most-active inhibitors. On the other hand, compounds 24 and 32 were related to lesser-active inhibitors. This can be inferred by their surroundings in the SOMs: compound 32 (log(10⁶/IC₅₀) = 5.77) was located in the upper-left zone where the majority of compounds have log(10⁶/IC₅₀) values between 3 and 4.99, while compound 24 is located in a complex zone where the majority of compounds present an inferior AChE inhibitory activity.

Despite the help that an analysis of the SOMs may offer a simple comparison of the structures of the outliers with the associated compounds in the data set does not reveal why these compounds are poorly predicted. As a result of the complexity of the modeled activity and the structural diversity, any obvious explanation is not feasible.

Comparison with previous computational models

Depicting the interaction of tacrine analogues with the AChE active site is not trivial. The flexibility of the ligands in conjunction with the rather large volume available to them inside the gorge allow incredibly high numbers of binding modes. There is experimental evidence that tacrine is able to bind also at a peripheral site [48]. This knowledge originated from dimers of tacrine which can be more potent and selective than tacrine since their tacrine units simultaneous bind to the catalytic and peripheral sites of AChE [49]. Since some of the tacrine analogues can bind at the catalytic site, while others can bind at the peripheral site, a broad QSAR model would represent an all-embracing relationship for subsets of compounds acting at different or both sites.

In general, previous SAR and QSAR studies have identified that hydrophobicity and the presence of ionizable nitrogen are essential features for the inhibitors to interact with AChE [13]. In addition, docking and molecular dynamics approaches confirm that the 3D dimensional positioning of the inhibitor in the active site of the enzyme, i.e. the mode of interaction, varies among different chemical classes.

Some previous QSAR models for AChE inhibition modeling include tacrine analogues. Recanatini et al. carried out a comparative QSAR analysis aimed at individuating the physico-chemical properties governing the inhibitory activity of 13 series of compounds including benzylamines, physostigmine analogues and 2 series of tacrines [50]. The QSARs for tacrine series were bilinear models in steric effects. However, the collinearity between steric and hydrophobic parameters did not allow the authors to draw any final conclusion about this model. In other work, Recanatini et al. derived the Hansch approach and CoMFA analysis for a series of 23 tacrine analogues substituted in positions‐6 and ‐7 of the acridine nucleus and bearing selected groups on the 9-amino function [4]. Both methods provided two separate models that show a satisfactory consistency, pointing out the negative steric effect of substituents in position‐7 and the relative steric freedom of position‐6 as main SAR aspects of the tetrahydroacridine-based AChE inhibitors. The report of Martín-Santamaría et al. includes 21 tacrines, 7 huprines, and 7 dihydroquinazolines in a QSAR modeling study [51]. These authors identified the key residues that modulate the inhibitory potencies of these classes of AChE inhibitors using the comparative binding energy (COMBINE) methodology. These authors report an interpretable COMBINE model that was able to fit and predict the activities of the three series of inhibitors reasonably well (Q² = 0.76). They also found a more robust predictive model when the same chemometric analysis was applied to the huprines set alone (Q² = 0.81), but the method was unable to provide predictive models for the other two families when they were treated separately from the rest.

QSAR studies may be classified as interpretive and predictive depending on the purpose of the study [52]. Interpretive studies often use a relatively small number of compounds and molecular descriptors that can be easily related to structural characteristics. They attempt to illustrate how the descriptors found to be important in a model relate to the interactions between the ligand and target. Predictive QSAR studies usually use a large and diverse dataset and computationally efficient descriptors. In distinction of the above-mentioned earlier interpretive studies, our QSAR model is predictive, since it encompasses a diverse set of tacrine analogues, therefore, our model can be useful like a predictive tool more than for mechanistic evaluation of this information for molecular design. In consequence, we employed molecular descriptors chosen for their computational efficiency and information-rich character. The combination of BRANNs and GA leads to a highly satisfactory model taking into account its predictive ability.

Conclusions

In conclusion, a large and diverse data of tacrine analogues has been assembled for which AChE inhibitory activity has been assessed. A nonlinear QSAR model was developed with descriptors generated from 3D molecular structure using BRGNN methodology. The model successfully explains the modeled structure-activity relationship, based on both statistical significance and predictive ability. In analogy with previous studies, this work has generated a system capable of rapid virtual screening of tacrine-related compounds for AChE inhibition. Unlike previous studies, this model is derived from an ample and varied set of tacrine analogues which integrates the current trends of tacrine modification. In this sense, our model provides more robust extrapolation in chemical space compared to models created in previous reports.