Profiling the structural determinants of aminoketone derivatives as hNET and hDAT reuptake inhibitors by field-based QSAR based on molecular docking

2.1Datasets preparation

This work collected 60 and 63 aminoketone derivatives with active data of hNET and hDAT from CHEMBL database [25] for field-based QSAR study as shown in Table 1. The datasets were randomly assigned 70% into training set of 42 molecules for hNET and 44 molecules for hDAT to generate QSAR models, and the rest of them as test sets were used to validate. The process of dividing datasets for both models should be considered the uniform distribution of compounds with as enough structural diversity and wide range activities as possible. The training and test sets of both models are demonstrated in Table 1.

2.2Homology modeling and molecular docking

The potent poses of compounds for each target were identified based on the complexes of target-ligand achieved by molecular docking and prime minimization. Since there is no crystal structural data reported in Research Collaboratory for Structural Bioinformatics Protein Data Bank (RCSB PDB) [26], the 3D structures of hNET and hDAT need to be constructed by homology modeling performed in SWISS-MODEL workspace [27]. The template used to construct models was confirmed according to sequence identity with targets by sequence alignment using ClustalW2 [23, 28, 29], and the stereochemical quality of built structures were checked by Ramachandran plot referenced our previous work [30, 31, 32].

The process of molecular docking was carried out in Maestro [33]. The first detailed work, active compounds of hNET and hDAT were collected from CHEMBL, and then their structures were preprocessed by the LigPrep [34] using OPLS-2005 force field [35]. The ionized state was assigned by Epik [36] at pH value of 7.0 ± 2.0. Secondly, to prepare hNET and hDAT structures for docking, the Protein Preparation Wizard module in Maestro [33] was used to add hydrogen atoms, assign partial charges using OPLS-2005 force field, assign protonation states and minimize the structure. The minimization was terminated when the root mean square deviation reached the maximum value of 0.30 Å. Finally, docking grid boxes were defined using the Receptor Grid Generation tool in Glide [37] by centroid of corresponding ligands from the template structures, and the prepared structures were docked in binding site of each target via standard precision (SP) in Schrödinger Glide module [37]. The preferred ligand was further rescored to calculate relative binding free energies using prime [38] MM-GBSA method at default.

2.3Field-based QSAR models construction and validation

Field-based QSAR studies were performed on Field-based QSAR module [39] in Maestro [33]. The potent poses were superimposed to bupropion by SMARTS method based on common scaffold in Flexible ligand alignment module. The field-based QSAR models were constructed by training sets at the extended Gaussian field, and the maximum PLS factors was set 5 for hNET and hDAT. The best models for each target were identified by statistical robustness, and their stabilities and predictive abilities were validated by the testing sets with leave-one-out (LOO) cross-validation methods. The most important statistics are the test set statistics including RMSE, Q2, and Pearson-r, which indicate how good the predictions are. If the predictions are not improving (much) as the number of PLS factors increases, the extra factors are not adding to the model and the model is probably over-fit. Scatter plots illustrated the correlation between the observed and predicted activity of all compounds for selected models. Moreover, the 3D contour maps for Gaussian steric, electrostatic, hydrophobic, hydrogen bond donor/acceptor, and aromatic ring fields were analyzed in QSAR Visualization module.

2.4ADME/T and chemical synthetize properties predicted for candidates

In order to obtain effective drugs, this work evaluated the ADME/T properties of the designed NDRIs by the pkCSM (http://biosig.unimelb.edu.au/pkcsm/prediction) and ADMETlab online server (http://admet.scbdd.com/calcpre/index/). On the SwissADME online server (http://www.swissadme.ch/index.php), the problematic fragments in novel candidates were identified by the complementary pattern recognition method of PAINS (pan assay interference compounds) and Brenk filter. Moreover, the synthetic accessibilities (SA) were evaluated to demonstrate the synthetic complexities for molecules. A molecule with SA score closer to 10 was considered difficult to synthetic, while was considered as easy with closer to 1.

3.1Homology modeling and validation for hNET and hDAT

The X-ray crystal structures of Drosophila melanogaster dopamine transporter (dDAT, PDB ID: 4XPH) showed high sequence identity with hNET (60.98%) and hDAT (55.83%), and showed higher sequence identity in binding site as the same as our previous work (69.05% and 78.57% for hNET and hDAT, respectively) [31]. Homology models of hNET and hDAT were constructed by dDAT as template, and the stereochemical quality and accuracy of models were illustrated by Ramachandran plot. Residues of models in allowed regions were 99.8% and 100% for hNET and hDAT respectively, indicating the reasonability of obtained models.

3.2The selected the potent poses by molecular docking and MM-GBSA prime

In this work, re-docking of 4XPH was first performed to validate the credibility of docking protocol. The binding pose of DCP in dDAT achieved by molecular docking was close to the original conformation with 0.474 of root mean square deviation (RMSD), indicating the validity of the docking protocol. Then, all active molecules were docked in each target using the same parameter settings as re-docking, and the obtained all complexes were minimized by MM-GBSA prime. Finally, the potent pose of bupropion was identified as the same as our previous work [31], including the electrostatic interaction between ammonium moiety and Asp75/Asp79 in hNET/hDAT recognition, and hydrophobic interactions formed in corresponding sub-sites for the other two groups. For other molecules, the docking poses were selected based on the docking score and the orientation of bupropion. Figure 1A and B show all selected poses gathered in the binding site surrounded by TM1, 3, 6, 8 and 10.

Figure 1.

Docking poses (A and B) and superimposed poses (C and D) template as bupropion of NDRIs compounds in the binding pocket of hNET and hDAT. The hNET (palecyan) and hDAT (lightblue) were shown in ribbon representation.

3.3Constructing and assessing of filed-based QSAR models for hNET and hDAT

All selected poses for hNET and hDAT were superimposed template as bupropion as demonstrated in Fig. 1C and D. In this work, 10 filed-based QSAR models were constructed using superimposed training sets compounds for hNET and hDAT, and the reliabilities of models were assessed by testing sets. All QSAR statistics are collected in Table 2.

For QSAR models of hNET, the values of R2 for the partial least square regression rose as increased factors, while the last 2 models should be discarded due to the lower stability value than the R2, which indicating the over-fitting of them. The best model of 3 factor performed good predictive and descriptive power based on QSAR statistics with 0.8479 for R2 and 0.5873 for Rcv2. Moreover, the 1.35e-15 for P value further indicated great degree of confidence. The reliability of model 3 was assessed by independent testing set of 18 compounds, and the value of Q2 for the predicted activities was 0.7352 suggesting its potent capability in predicting extend molecules. Besides, 0.8776 for Pearson-r showed a good correlation between predicted and observed activity for the test set, which further confirmed models’ reliability (robustness). As the similar analysis of QSAR statistics among hDAT models, the model of 3 factor was the best selection owing to the high linear correlation for training set (R2= 0.8040) and the potent capability in predicting testing set (Q2= 0.6266) without over-fitting. The experimental and predicted activities in data sets for selected models collected in Table 1 and evenly distributed around the diagonal as illustrated in Fig. 2, which intuitively verified the effectiveness of models.

Figure 2.

Correlation plots of the predicted and experimental PIC50 values for training and test sets of hNET (A) and hDAT (B).

Figure 3.

3D contour maps of 6 extended Gaussian fields referenced bupropion for hNET and hDAT field-based QSAR models. A-F and G-L were steric, hydrophobic, electrostatic, hydrogen bond acceptor, hydrogen bond donor and aromatic ring field contour maps based on bupropion in hNET (cyan) and hDAT (pink) models, respectively. The colors of favorable and unfavorable regions were green and yellow for steric fields, yellow and white for hydrophobic field, blue and red for electrostatic fields, red and magenta for hydrogen bond acceptor, purple and cyan for hydrogen bond donor and orange and gray for aromatic ring field, respectively.

3.4Contour maps analyses of QSAR models for hNET and hDAT

3D contour maps of 6 extended Gaussian fields referenced bupropion (compound 1) are illustrated in Fig. 3, which provide information about the contribution of substituents on phenyl-aminoketone scaffold to the activity for each target. As shown in Fig. 3, steric and hydrophobic fields performed more potent contribution to activity than other 4 fields with 0.539 and 0.225 respectively for hNET, and with 0.484 and 0.323 respectively for hDAT. Thus, the volume and hydrophobicity of substituents on phenyl-aminoketone played major roles in regulating activity.

From the steric contour maps (Fig. 3A), there were 2 green and 1 yellow contour regions around aminoketone moiety, with green meaning bulky groups are favored and yellow meaning bulky groups are disfavored. The green contour maps at R1 (substituted-carbonyl) were incarnated by the increased activities of compounds 1 (5.73), 40 (6.22) and 42 (6.43) with methyl, ethyl and propyl group respectively. However, compounds 46 (6.40), 26 (5.70) and 29 (5.31) with butyl, amyl and hexyl are gradually weaker than 42 for hNET, which demonstrated the larger stereoscopic space is not meaning blindly the longer carbon chain. This opinion was effectively proved by compounds 30 (6.44) and 46 (6.40) with the same carbon atoms for isobutyl and butyl group, the same condition was found in compounds 7 (7.06) and 21 (7.95) with propenyl and propyl. Moreover, the green and yellow contour maps at R3 (substituted-ammonium) signified the large group with short carbon chain was beneficial to increase the activity. It was reflected by the weak activity of compound 19 (5.60) with butyl rather than bulky tertiary butyl in compound 1, and by compound 53 (5.09) with cyclopentyl-ammonium comparing with the larger piperidyl group in 52 (6.14). The similarities of steric contour maps between hDAT (Fig. 3G) and hNET at R1 and R3 meant that the large moieties in these regions are favor of increasing activities for hDAT.

The yellow contour map of hydrophobicity was mainly sited at R2 (substituted aromatics) as illustrated in Fig. 3B and H, which demonstrated introducing the higher hydrophobic group would increase the activities against targets. The result was confirmed by 1 (7.91 for hNET, 7.35 for hDAT) with 4-methanol group versus 14 (7.06, 7.30) replaced by 4-phenolic hydroxyl, 9 (7.33, 7.50) comparing with 15 (6.77, 6.73) and 21 (7.95, 7.79) in contrast to 2 (7.64, 6.81). Moreover, the larger yellow polyhedron covering meta-position of R2 for hDAT (Fig. 3H) comparing with hNET (Fig. 3B), which supplied guidance to regulate the balance of activity for two targets. Conversely, the white contour at R3 (Fig. 3B and H) indicated that introducing weaker hydrophobic group would improve the inhibition of targets. For example, the replacement of cyclopentyl in this region by piperidyl with lower hydrophobicity could increase the activity. This regulation could be seen from compound 52 (6.14, 6.34) versus 53 (5.09, 5.64), 54 (6.24, 6.38) in contrast to 55 (5.06, 6.08), respectively.

The electrostatic field contour maps are illustrated in Fig. 3C and I. The regions in red surrounded at para-position of R2 implied where negative charge groups enhanced bioactivity, whereas the regions in blue sited at meta-position of R2 purported positive charge groups improved the activity. The orderliness for para-position of R2 could be adequately reflected by the decrease activities of compounds 21 (7.95, 7.79), 2 (7.64, 6.81) and 22 (6.26, 5.96), which is due to the gradual reduction of electronegativities for methyl, methyl formate and nitro, respectively. Meanwhile, it also appeared in compounds 9 (7.33, 7.50), 8 (7.08, 7.4), 15 (6.77, 6.73) with increasing of electronic deficiency for iodine, bromine and fluorine atom. At meta-position of R2, because of the more electropositive of carbon atom in naphthalene at blue region, compound 18 (7.93, 7.40) performed higher activity contrasting to compound 2 (7.64, 6.81) with methyl or 23 (6.5, 7.17) with oxygen atom in carbonyl. Additionally, the oxygen atom with negative charge in meta-nitro of compound 58 intersected the blue area, therefore, it performed the weakest active among similar compounds 54, 56, 57, 59, 60.

The H-bond acceptor contour maps (Fig. 3D and J) closed to R2, which indicated that introduction of H-bond acceptor atom in red region is favorable while magenta areas stood for the opposite. This might be one reason for compound 13 with hydroxy as H-bond acceptor atom performed potent activity for hNET and hDAT. The purple contours in Fig. 3E and K signify H-bond donor was beneficial, while cyan contours were unfavored to increase activity. This might could be explained by the hydrogen bond between protonated nitrogen (NH+) in R3 and negative charged oxygen (O-) of Asp75/Asp79 were essential for drugs recognition of hNET and hDAT in previous work [31, 40]. Finally, from gaussian aromatic ring contours as displayed in Fig. 3F and L, the gray polyhedron surrounding the para-position of R2suggesting that it was not prominently effectual to increase the activity by introducing the aromatic group.

Considering all discussed above, this work supplied the graphical description of the structural features for designing novel NDRIs. In order to improve the activity, researchers could replace the methyl by the bulky groups with suitable carbon chain at R1, and/or with additional weak hydrophobicity at R3. Moreover, hydrophobic and electronegative substituted at R2 should be helpful for enhancing the effectiveness. Meanwhile, designer could regulate the active balance between hNET and hDAT by changing the hydrophobicity with obvious contribution to activity of meta-position of R2. At last but the most important, it was necessary to retain the protonizable nitrogen in ammonium group to form hydrogen bond for recognizing by targets.

Table 3

The predicted activities of designed novel derivatives by selected hNET and hDAT QSAR models

‡Molecules with active as the high active NDRI compound 10. §Molecules with more active than compound 10.

3.5Design of novel derivatives

According to the information based on compound 1 provided by Field-based QSAR models, firstly this work replaced methyl of R1 by cyclopentyl, cyclopentadienyl, cyclohexyl and phenyl, meanwhile introduced cyclopentyl, cyclohexyl, pyrrolyl, piperidyl at R3. Secondly, the designed compounds were optimized by bring methyl, trifluoromethyl, methoxy, trifluoromethoxy and formate group in R2, and regulated the active balance for hNET and hDAT by changing the hydrophobicity of substituent between meta- and para-position. Additionally, Other useful proposals could be performed to import the electronegative group in para-position of R2 to yield goo candidates, such as methoxy, sulfo group and so on. Finally, 29 novel derivatives were designed and their activities were predicted by selected models of hNET and hDAT as collected in Table 3. All designed candidates performed with hNET > hDAT potency order, and their inhibitive ability were no less than the highest active NDRI compound 10.

3.6ADME/T profiles and medicinal chemistry friendliness properties of candidates

The ADME/T properties and medicinal chemistry of 29 high active candidates were predicted using pre-ADMET web servers, and the results are shown in Table 4. As results, the rang of HIA was from 89.056% to 98.477%, indicating that all candidates performed high absorption properties from intestine. The VDss of most molecules was between 0.71 and 2.81 L/kg, which meant that they were evenly distributed in the plasma and in tissues, while compounds D21, D25, D27 and D29 with low 0.71 L/kg were distributed more in the plasma rather than in tissues. There was no compound to be considered as poorly distributed to the brain and unable to penetrate the CNS according BBB and CNS permeability. For metabolism, all designed molecules were predicted as substrates for the CYP450 2D6/3A4/1A2 subtype, which indicated they could be metabolized by P450 isoforms. In addition, compound D6, D16, D25, D26, D27 and D29 could not inhibit the CYP450 family, whereas the other compounds might inhibit the CYP450 1A2/2C19/2C9/2D6/3A4 subtype. The predicted total clearance of all compounds demonstrated they could be cleared by combing hepatic and renal tissue. Based on the predicted the toxicity, some candidates would be harmful for the liver, and only a few compounds were likely to be associated with skin sensitization, and none of them performed mutagenic potency. Finally, as shown the medicinal chemistry friendliness properties in Table 4, there was no problematic fragment in novel compounds assessing by PAINS and Brenk filter. The SA of all designed candidates ranged from 2.85 to 3.58, indicating the low synthetic complexity of structure and high synthetic feasibility. Taking all predicted ADMET properties and medicinal chemistry profiles into consideration, the designed 29 novel derivatives could be served as lead compounds for further development of NDRIs.