|Home | About | Journals | Submit | Contact Us | Français|
Back-propagation artificial neural networks (ANNs) were trained on a dataset of 104 VMAT2 ligands with experimentally measured log(1/Ki) values. A set of related descriptors, including topological, geometrical, GETAWAY, aromaticity, and WHIM descriptors was selected to build nonlinear quantitative structure-activity relationships. A partial least squares (PLS) regression model was also developed for comparison. The nonlinearity of the relationship between molecular descriptors and VMAT2 ligand activity was demonstrated. The obtained neural network model outperformed the PLS model in both the fitting and predictive ability. ANN analysis indicated that the computed activities were in excellent agreement with the experimentally observed values (r2 = 0.91, rmsd = 0.225; predictive q2 = 0.82, loormsd = 0.316). The generated models were further tested by use of an external prediction set of 15 molecules. The nonlinear ANN model has r2 = 0.93 and root-mean-square errors of 0.282 compared with the experimentally measured activity of the test set. The stability test of the model with regard to data division was found to be positive, indicating that the generated model is predictive. The modeling study also reflected the important role of atomic distribution in the molecules, size, and steric structure of the molecules when they interact with the target, VMAT2. The developed models are expected to be useful in the rational design of new chemical entities as ligands of VMAT2 and for directing synthesis of new molecules in the future.
Methamphetamine (METH), an amphetamine derivative, is an addictive psychostimulant drug and a significant health concern due to its abuse liability and potential neurotoxic effects.1 Chronic use of METH may cause long-term neural damage in humans, with concomitant deleterious effects on cognitive processes, such as memory and attention.2 Despite the serious consequences of METH abuse, currently there is no FDA approved clinical treatment for METH addiction. Thus, there is an increasing interest in identifying the underlying mechanisms of METH action, as well as the relevant pharmacological targets to promote the development of novel therapeutic agents as treatments for METH abuse.
The abuse liability of METH and structurally-related amphetamine compounds is thought to be due to alterations in dopaminergic neurotransmission.3,4 In this respect, the dopamine transporter (DAT) and the vesicular monoamine transporter-2 (VMAT2), presynaptic proteins critical for dopamine storage and release, are the primary targets for METH action.3,4,5 Specifically, METH interacts with VMAT2 to release dopamine from the synaptic vesicles into the cytosol of the presynaptic terminal,6,7 METH also inhibits monoamine oxidase and evokes the release of dopamine from the cytosol into the extracellular space via reverse transport of DAT, leading to an increase in dopamine concentration in the extracellular space.7,8,9
Lobeline, an alkaloidal constituent of Lobelia inflata LINN, is a nicotinic receptor ligand with high affinity for α4β2* nicotinic receptors.10 Lobeline was previously investigated as a therapeutic agent to treat tobacco dependence.10 Recent study indicated that lobeline has both temperature-dependent and temperature-independent neuroprotective effects against METH toxicity.11 Lobeline also inhibits dopamine uptake and promotes dopamine release from storage vesicles within the presynaptic terminal via an interaction with the tetrabenazine binding site on VMAT2.12 Lobeline attenuates d-amphetamine- and methamphetamine-induced hyperactivity, and inhibits the discriminative stimulus effects and self-administration of methamphetamine.13,14 However, lobeline does not support self-administration in rats,15 suggesting a lack of addiction liability. Thus, the development of lobeline and lobeline analogs with targeted selectivity at VMAT2, represents a novel approach for the treatment for psychostimulant abuse.10, 12
To date, very few VMAT2 ligands have been reported in the literature; these include low affinity ligands, such as 3-amino-2-phenylpropene derivatives,16 and high affinity tetrabenazine derivatives.17,18 Tetrabenazine was introduced in 1956 as an antipsychotic agent 19 and has currently been submitted for FDA approval as an anti-chorea drug.20 Recently, a small library of structural analogs of lobeline have been synthesized, and their activity and selectivity for VMAT2 have been evaluated.9,21,22,23
In the discovery of novel and more potent and selective lobeline analogs, we consider computational modeling as a valuable aid in drug design and optimization. In this respect, the nature of the interaction of these novel ligands with the binding site(s) on VMAT2 is not known due to the lack of crystal structure for this protein. Thus, a structure-based drug design approach is not available. On the other hand, neural network analysis approach, particular back-propagation network to data analysis, has received much attention over last decade. This artificial system emulates the function of the brain, in which a very high number of information-processing neurons are interconnected and are known for their ability to model a wide set of functions, including linear and non-linear functions, without knowing the analytic forms in advance.24 The rapid advancement of computing systems in the past 20 years is an important factor leading to the success of this approach in various engineering, business, and medical applications. So far, the neural network approach has been applied in a variety of biomedical areas, which includes analysis of appendicitis25 and cancer imaging extraction and classification,26 AIDS research and therapy.27 This approach has also been used in drug design and discovery;24 as well as in pharmaceutical applications such as pharmaceutical production development,28 pharmacodynamic modeling29 and mapping dose-effect relationships on pharmacological response.30
In this current study, the neural network analysis approach is used to build a quantitative structure-activity relationship (QSAR) model on a set of 104 tetrabenazine and lobeline analogs with known affinity for VMAT2. This is the first QSAR modeling study that addresses the interaction of a small library of ligands at the VMAT2 binding site. The goals of the current work are (i) to extract the relevant descriptors to establish the QSAR of the library of ligands, (ii) to establish the high predictive power of neural network modeling on this library of compounds, and (iii) to develop insights regarding the relationship between the descriptors of the compounds and their affinity for VMAT2. The developed models are expected to be valuable in the rational design of chemical modifications of first-generation VMAT2 ligands in order to identify the most likely candidates for synthesis and discovery of new lead compounds.
Molecular modeling was carried out with the aid of the Sybyl discovery software package.31a The software was used to construct the initial molecular structures utilized in the geometry optimization (energy minimization) for all molecules evaluated in this study. The geometry optimization was first performed by using the molecular mechanics (MM) method with the Tripos force field and the default convergence criterion. Since the pKa values of the basic nitrogen atoms in these compounds are between 8 and 10, and the synaptic vesicles are acidic (pH 5-5.6)32, most microspecies of the compounds are expected to be protonated when binding to VMAT233. Thus, in the construction of the initial molecular structures, the basic nitrogen atoms in these compounds were protonated with a formal charge of +1 assigned to the positively charged nitrogen atom. In this respect, the crystal structures of meso-transdiene (MTD, compound T78 in Table 1) shown in Figure 1, the structure of lobeline in crystal structure of acetylcholine-binding protein from aplysia californica in complex with lobeline (PDB code 2BYS), and the crystal structure of (−)-α-9-O-desmethyldihydro- tetrabenazine34a were utilized as a reference. It has been shown that only the (+)-isomer of dihydrotetrabenazine exhibited high affinity for VMAT2.34b Protonation of nitrogen atom in these molecules is also relevant, considering that some of the lobeline analogs in the dataset include quaternary ammonium nitrogen atoms, such as T6, P6 and P7 in Table 1. All of the obtained conformations optimized at the MM level were further refined to their lowest energy states with MOPAC, a semi-empirical molecular modeling routine, utilizing the PM3 Hamiltonian.
The 104 molecules listed in Table 1 constitute a database for the structure-activity relationship analysis. These data are a combination of 10 previously reported compounds with experimental Ki values,17,18 and 94 compounds with experimental Ki values from our own laboratory.9,21–23 A dataset of 89 molecules (T1–T89) was used for model training and leave-one-out (LOO) validation. A dataset of 15 molecules (P1–P15) from a different compound series was used for external testing. Table 1 also lists the experimental Ki values for each of these compounds, providing the pharmacological parameter which characterizes the interaction with the VMAT2 binding site. For the set of the 104 molecules utilized, the Ki values of 13 molecules were ≤ 0.01 μM, 10 molecules had Ki values in the range 0.1–1 μM, 63 molecules had Ki values in the range 1–10 μM, and 14 molecules had Ki values ≥ 10 μM.
The optimum three-dimensional conformations were used for generation of descriptors, some of which were geometry-dependent. A total of 807 molecular descriptors, consisting of zero-dimensional (constitutional descriptors), one-dimensional (functional groups, empirical descriptors, physical properties), two-dimensional (topological descriptors), as well as three-dimensional (geometrical, WHIM, GETAWAY, aromaticity descriptors) variables, were created by the DRAGON program.35a A reduced set of 149 descriptors was obtained after the constant and near constant descriptors and the highly inter-correlated (>0.95) descriptors were discarded.
Partial least squares (PLS) analysis was performed using the QSAR module of Sybyl version 7.0, with the NIPALS algorithm to extract the original variable into PLS components. All variables were initially auto-scaled to zero mean and unit variance. This method produces new variables by a linear combination of the original descriptors and uses them to predict the biological activities. The advantage of this method is that the method can be used for strongly correlated, noisy data with numerous independent variables (e.g. in modeling data sets where the number of descriptors greatly exceeds the number of observations).36
Experimental Ki values of lobeline and its analogs were measured according to the procedure described by Zheng et al.9 Experimental Ki values of tetrabenazine and its analogs, except for tetrabenazine and Ro4-1284 (P13 and P14 in Table 1 respectively) were taken from Lee et al.17 and Conney et al.18 The Ki values for tetrabenazine and Ro4-1284 were also measured in our laboratory to compare the Ki values with those from the other laboratories. The Ki values for P13 and P14 obtained from our own assays (0.013 μM and 0.028 μM, respectively) were used in this study. These Ki values are very similar to the literature values for these compounds (i.e., 0.0081 μM and 0.042 μM, respectively) reported by Lee et al.17 Conney et al. 18 also reported a Ki value of 0.0067 μM for tetrabenazine. The log(1/Ki), with Ki values expressed as molar, was used as the target pharmacological criterion to derive the QSARs.
Feed-forward, back-propagation-of-error networks were developed using a neural network C program.37 Network weights (Wji(s)) for a neuron “j” receiving output from neuron “i” in the layer “s” were initially assigned random values between −0.5 and +0.5. The sigmoidal function was chosen as the transfer function that generates the output of a neuron from the weighted sum of inputs from the preceding layer of units. Consecutive layers were fully interconnected; there were no connections within a layer or between the input and the output. A bias unit with a constant activation of unity was connected to each unit in the hidden and output layers.
The input vector was the set of descriptors for each molecule in the series, as generated by the previous steps. All descriptors and targets were normalized to the [0,1] interval using the following formula:
Where Xij and Xij′ represents the original value and the normalized value of the j-th (j=1,…k) descriptor for compound i (i=1,…n). Xmin and Xmax represent the minimum and maximum values for the j-th descriptor. The network was configured with one or more hidden layers. During the ANN learning process, each compound in the training set was iteratively presented to the network. That is, the input vector of the chosen descriptors in normalized form for each compound was fed to the input units, and the network’s output was compared with the experimental “target” value. During one “epoch”, all compounds in the training set were presented, and weights in the network were then adjusted on the basis of the discrepancy between network outputs and observed log(1/Ki) values by back-propagation using the generalized delta rule.
Models were cross-validated using the “leave-one-out” approach.38
Generated models were then tested using a subset of 15 compounds (P1–P15 in Table 1). They were randomly selected to cover the experimental activity range as uniformly as possible, and were not used in either the variable selection or the model building processes.
QSAR models were assessed by Pearson correlation coefficient r2, root mean square deviation (rmsd), and predictive q2, which is defined as
Where SD is the sum of squared deviations of each measured log(1/Ki) value from its mean, and PRESS is the sum of squared differences between actual and predicted values.
To provide a comparison with the neural network model, a linear PLS analysis was carried out on the affinity for the training set of 89 compounds. Variable selection from 149 descriptors generated as described in section 2.2 was used to build the PLS model. After stepwise exclusion of low contribution variables, 12 variables were included in the model, resulting in an equation in which only variables that significantly increased the predictability of the dependent variable were included. During the PLS analysis, five components (latent vectors) that explain the most covariance between the 12 descriptors and experimental log(1/Ki) values are obtained. The linear model built with these five components corresponds to the highest q2 (0.780) and lowest LOO rmsd (0.353). A final set of 12 descriptors which gave the best PLS model is listed in Table 2.
To identify the best neural network QSAR model, an exhaustive search of all possible models with various numbers of descriptors was carried out.24c The 12 descriptors used in the PLS model were first examined based on different neural network configurations. Of these descriptors, seven having small contributions were removed. A step forward descriptor selection procedure was used to select additional descriptors which are important to the neural network analysis. To find the optimal number of neurons in the hidden layer, neural network architectures with n: h: 1 where n = 8 to 12 and h = 2 to 7, were respectively trained. The most promising descriptor combinations and network configuration were selected by the LOO cross- validation procedure. The best descriptor combinations with different numbers of descriptors and configurations selected by the internal LOO validation processes were applied to the 15 external compounds, which were put aside at the beginning of the analysis. The same training and (internal and external) validation sets were always used for model quality comparison. The calculated internal q2 as well as the LOO root mean square derivation (loormsd), and the external r2 and rmsd values for these molecules were used to demonstrate the predictive ability of the selected descriptor combinations and network configuration. As seen in Table 3, the predictivity of the neural network models increases from 8 to 11, and drops thereafter with the 11 descriptors. A network with more weights was generally trained to their best configuration with a relatively lower number of epochs, with the exception of NN8-4-1, which has an irregular loormsd variation between 0 and 100000 training epochs. When the input vector of the neural network was fixed with the 11 descriptors, its fitting ability is good. Further examination shows that whereas r2 from the whole training set always increases with increasing numbers of hidden neurons, Figure 2 and Table 3 shows that the 11:3:1 neural network architecture exhibited a good internal validation and externally predictive performance, when the optimal training epochs was set to 200000 (Figure 3).
A common belief about neural network systems is that the number of parameters in the network should be related to the number of data points in the dataset, and the expressive power of the network.39 For the current dataset, the results from the search for various configurations (Table 3) shows that the weights used in the neural network configuration to obtain better internal and external prediction are between 33 and 40. Comparing the NN11-3-1 model with the NN8-4-1 model, both configurations include 36 weights. While the numbers of descriptors used in both configurations are not overwhelmed (i.e., the ratio of the number of compounds in the training set to the number of descriptors is 8 and 11, respectively for NN11-3-1 and NN8-4-1), the former model includes more information from compounds themselves (i.e., a greater number of descriptors was used), and provides higher internal and external prediction ability (Table 3). Therefore, the best ANN model is an 11-descriptor model with 3 hidden neurons consisting of the following descriptors: Ram, PW5, LP1, SEige, VEp2, DISPm, G(N..N), H7m, RARS, R1p+, HOMT with 3 hidden nodes. Brief summaries of these descriptors can be found in Table 2.
The best PLS model for calculated descriptors is a 5-component model of 12 descriptors (Ram, PW5, DISPm, G(N..N), RCON, DISPp, HATS5u, R5u, R4u+, R1v+, R5v+, R5e+). The r2 values for the training and for the LOO cross-validation runs are 0.83 and 0.78, respectively; the corresponding rmsd values are 0.303 and 0.353, respectively. For the 15 test compounds, the r2 value is 0.87 and the rmsd is 0.415 (Table 4). Figure 4 shows the relationships of the trained, LOO and external predicted log(1/Ki) values versus the experimental log(1/Ki) values for the PLS model. The calculated log(1/Ki) values by the PLS model for the 104 molecules are shown in Table 1.
The best 11-descriptor neural network model consists of the following descriptors: Ram, PW5, LP1, SEige, VEp2, DISPm, G(N..N), H7m, RARS, R1p+, HOMT. The statistical results for this model are as follows: r2 =0.91, rmsd =0.225, q2 = 0.82, and loormsd = 0.316. The QSAR model demonstrated good predictivity in the test set: r2 =0.93, and rmsd = 0.282 (Table 4). Comparing the values of correlation coefficients and root mean square deviations listed in Table 4, the ANN model is better than the PLS model. The predictive ability of NN11-3-1 can also be judged from the plots of the trained, LOO and external test predicted versus experimental log(1/Ki) values shown in Figure 5. The computed log(1/Ki) values of the model for the 104 molecules in the database are listed in Table 1.
The statistics of the PLS model are summarized in Table 5. The X weights obtained from the PLS analysis are displayed in Table 6; the definition of the descriptors used in the model can be found in Table 2. Since q2 beyond the 5th component basically remains a constant for the PLS model, the most important descriptors may be located from the first 4 components only. In component 1, it is clear that PW5 and R5v+ are the most heavily weighted descriptors. PW5 is topological descriptors related to molecular shape;35h R5v+ is one of the GETAWAY descriptors. According to the definition, GETAWAY descriptors are a type of descriptor encoding both geometrical information given by the influence molecular matrix and the topological information given by the molecular graph, weighted by chemical information encoded in selected atomic weightings.35f GETAWAY descriptors are related to a molecular 3D structure. In the 2nd component, the most weighted descriptors are two 3D descriptors - RCON and G(N..N), among which G(N..N) is a 3D geometrical descriptor and characterizes the geometrical distance between nitrogen and nitrogen in those molecules that contain more than one N-atom; RCON is referred to R-connectivity index, which is sensitive to the molecular size, conformational changes and cyclicity.35f Accordingly, other relatively important descriptors identified by the component analysis were HATS5u, DISPp, R5e+ and DISPm. HATS5u and R5e+ are GETAWAY descriptors; DISPp and DISPm are geometrical descriptors. The details of these descriptors can be found in the previously published study.35(b)–(h)
The descriptors important for the NN11-3-1 model are plotted in Figure 6 according to a previous described procedure.24c, 40 The three most important descriptors are H7m, G(N..N) and DISPm. These descriptors indicate that molecular size, shape and atomic distribution in the molecules are important. The most important descriptor is the GETAWAY descriptor H7m. The second most significant descriptor, G(N..N), which is also an important descriptor in the PLS model. DISPm is related to the molecular geometry as well as molecular size; the other two important topological descriptors, PW5 and VEp2, also emphasize these important features. Among these important descriptors, PW5, G(N..N), and DISPm are common to both the PLS and the NN11-3-1 models. Analysis of the relationship between the values of the descriptors and the experimental endpoints of the compounds, several features were observed. To obtain an VMAT2 binding affinity Ki < 0.1 μM, it is shown that PW5 needs to be 0.12; H7m needs to be > 0.05 (most are between 0.05 and 0.15); VEp2 needs to be between 0.18 and 0.20; DISPm needs to be between 7.5 and 10.85; G(N..N) needs to be 0 or 4.34; and R1p+ needs to be between 0.03 and 0.05.
The linear PLS model and the NN11-3-1 model have five descriptors in common (Ram, PW5, DISPm, G(N..N), RCON). These descriptors are expected to significantly encode the linear relationship between the variables and the target bioactivity values. The neural network model contains five topological (Ram, PW5, LP1, SEige, VEp2) and six 3D descriptors (DISPm, G(N..N), H7m, RCON, R1p+, HOMT). These descriptors contain chemical information concerning size, symmetry, shape and distribution of the molecular atoms in the molecule. When compared to the linear model, the data show that the neural network model has been able to capture a more detailed analysis of the structure-activity relationships, and affords a high correlation with low root-mean-square-deviation.
To test the stability of the best neural network and PLS model, the dataset of 104 molecules was randomly divided into two sets; one set of 89 molecules for training and LOO validation, and the other set of 15 molecules for external testing. This process was randomly repeated five times. Among each division, the dissimilarity of one test set from the test set of another division was greater than 90%, i.e., less than one compound in 15 in one test set was the same as the compounds in another test set. The same eleven descriptors and the same neural network architecture, i.e., eleven input neurons, three hidden neurons and one output neuron, were used to train and test the neural network models. Learning epochs was set to 200,000. For the PLS model, the same twelve descriptors listed in Table 2 were used. Interestingly, during the five times training with the different set of eighty-nine molecules to build the PLS models, it was observed that the best PLS models are always those utilizing only five components of the twelve descriptors models for the five sets of data. The calculated results were listed in Table 7.
The data in Table 7 indicate that the generated models are stable with regard to the data division. The neural network model has a higher predictive power than the PLS model. These results are consistent with the conclusions drawn previously (see section on descriptor contribution analysis).
In this study, PLS and neural network approaches were used to build linear and nonlinear QSAR models for a set of tetrabenazine and lobeline analogs that are ligands at the VMAT2. These are the first models to predict Ki values for the interaction of these two families of compounds with VMAT2. While the linear PLS model is predictive, it was demonstrated that the fully interconnected three-layer neural network model trained with the back-propagation procedure was to be superior in learning the correct association between a set of relevant descriptors of compounds and these log(1/Ki) for VMAT2. The trained neural network model (NN11-3-1) including eleven-input and three-hidden neurons exhibited a high predictive power (r2=0.91, rmsd=0.225 and q2=0.82, loormsd=0.316). This model succeeded in predicting the Ki values of an additional set of 15 tetrabenazine and lobeline analogs which were not included in the model training (external r2=0.93, rmsd=0.282). The stability test of the model with regard to data division was found to be positive. Evaluation of the contributions of the descriptors to the QSAR reflected the importance of atomic distribution in the molecules, molecular size and steric effects of the ligand molecules, when interacting with their target binding site on the VMAT2. The nonlinear relationship between these factors and the endpoint bioactivity values has been clearly demonstrated. These results indicate that the generated neural network model is reliable and predictive. Thus, this new neural network model, reported herein, will be valuable for future rational design of novel second generation ligands targeted to VMAT2, aimed at developing novel therapeutics for the treatment of methamphetamine abuse.
This work was supported by NIH Grant No. DA013519.
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.