Convergent QSAR studies on a series of NK₃ receptor antagonists for schizophrenia treatment

Abstract

The dopamine hypothesis states that decreased dopaminergic neurotransmission reduces schizophrenia symptoms. Neurokinin-3 receptor (NK₃) antagonists reduce dopamine release and have shown positive effects in pre-clinical and clinical trials. We employed 2D and 3D-QSAR analysis on a series of 40 non-peptide NK₃ antagonists. Multivariate statistical analysis, PCA and HCA, were performed to rational training/test set splitting and PLS regression was employed to construct all QSAR models. We constructed one highly predictive CoMFA model (q²= 0.810 and r²= 0.929) and acceptable HQSAR and CoMSIA models (HQSAR q²= 0.644 and r²= 0.910; CoMSIA q²= 0.691, r²= 0.911). The three different techniques provided convergent physicochemical results. All models indicate cyclopropane, piperidine and di-chloro-phenyl ring attached to cyclopropane ring and also the amide group attached to the piperidine ring could play an important role in ligand–receptor interactions. These findings may contribute to develop potential NK₃ receptor antagonists for schizophrenia.

• CoMFA
• CoMSIA
• HQSAR
• NK₃ receptor antagonists
• schizophrenia

Introduction

Schizophrenia is a psychotic disorder that affects 24 million people around the world, with a major peak of appearance between 15 and 35 years of age1. The most accepted theory for schizophrenia is the dopamine hypothesis, which states that the disorder results from increased dopaminergic neurotransmission. According to this view, pharmacological approaches that decrease dopaminergic neurotransmission in the brain should alleviate psychotic symptoms2^,3. Therefore, all drugs currently employed to treat schizophrenia act, to some extent, on dopamine receptors, mainly as antagonists of D2-type receptors. However, antagonism of D2 receptors in the nigrostriatal pathway is responsible for severe side effects of antipsychotics, such as Parkinsonism, acute dystonias, akathisia, and tardive dyskinesia3.

In the search for molecules with better side-effect profiles, neurokinin-3 receptor (NK₃) has been investigated as a new target for schizophrenia. Stimulation of this receptor is known to increase dopamine release and, therefore, binding of antagonists may decrease dopaminergic neurotransmission without interfering with dopamine receptors2^,3. Neurokinin receptors are G protein-coupled receptors (GPCRs) activated by neurokinins, peptides from the tachykinin family, which includes substance P, neurokinin A and neurokinin B. NK₃ receptor is mainly activated by neurokinin B. The exact distribution of this receptor in the brain is still unknown, but it seems to be found predominantly in cortical regions such as the frontal and parietal cortex, and also in the hippocampus3–5. Data from pre-clinical and clinical trials of non-peptide NK₃ antagonists used to treat schizophrenia have shown positive effects, with significant reduction of symptoms and lack of major side effects3^,6–8.

Of the many available drug design strategies, quantitative structure–activity relationship (QSAR) models are widely employed due their advantages, such as the ability of mathematical models to predict the biological activity of untested compounds, providing physicochemical insights into molecular optimization (for example, the influence of certain substituent or molecular properties on the biological activity, ligand–target interactions, etc.)8–14. To this date, there are few quantitative structure and activity relationship studies for NK₃ antagonists15^,16; therefore, there is a lack of information in this research field.

Here, we constructed QSAR models employing molecular hologram quantitative structure activity relationship (HQSAR) as a 2D-QSAR approach, comparative molecular fields analysis (CoMFA) and comparative molecular similarity indices analysis (CoMSIA) as 3D-QSAR approaches. The models presented good convergence towards physichochemical properties influencing activity. Significant information was extracted regarding structure–activity relationships, contributing to define important structural features for NK₃ antagonism to be used in the discovery of novel antipsychotic drugs. Furthermore, the compounds studied in this work showed high affinity to NK₃ receptors at nanomolar range and some compounds also presented good balance between binding affinity and pharmacokinetics profile related to CNS effects4. Therefore, the constructed QSAR models highlight the main physicochemical properties related to NK₃ antagonism of a promising drug candidates series.

Material and methods

Biological data, molecular alignment and descriptors calculations

The data set contains 40 NK₃ receptor inhibitors (Table 1) synthesized and in vitro evaluated (K_i) by Juhl et al.4. The pK_i values (−log K_i) were used as the dependent variable in the QSAR models, spanning approximately three logarithmic unities. Visual inspection demonstrated that pK_i values are bimodally distributed, with an adequate number of compounds with low, moderate, and high pK_i values.

Table 1. 2D structure of data set compounds with its K_i and pK_i values.

Training set
N-benzyl amide analogues
Cpd	Stereochemistry	Ar	Ki (nM)	pKi
1	(1R,2S)	C6H5	83	7.081
2	(1S,2R)	4-F-C6H4	5.8	8.237
Aromatic ring on cyclopropane core
Compound	Stereochemistry	Ar	Ki (nM)	pKi
3	(1S,2R)	4-F-C6H4	67	7.174
4	(1S,2R)	4-Cl-C6H4	36	7.444
5	(1S,2R)	4-CF3-C6H4	1000	6.000
6	(1S,2R)	4-CN-C6H4	270	6.569
7	(1S,2R)	3-F-C6H4	140	6.854
8	(1S,2R)	3-Cl-C6H4	32	7.495
9	(1S,2R)	3-Me-C6H4	40	7.398
10	(1S,2R)	3-OMe-C6H4	56	7.252
11	(1S,2R)	3,4-F2-C6H3	47	7.328
12	(1S,2R)	3,4-Cl2-C6H3	5.8	8.237
13	(1S,2R)	3,4-Me-C6H3	7.3	8.137
14	(1S,2R)	3,4-OCH2O-C6H3	140	6.854
15	(1S,2R)	2-tiofenil	93	7.032
16	(1S,2R)	3-tiofenil	110	6.959
Phenyl piperidine analogues
Cpd	Stereochemistry	R	Ki (nM)	pKi
17	(1S,2R)		4.2	8.377
18	(1S,2R)		6.1	8.215
19	(1S,2R)		1.5	8.824
20	(1S,2R)		5.6	8.252
21	(1S,2R)		12.0	7.921
22	(1S,2R)		3.6	8.444
23	(1S,2R)		7.8	8.108
Piperazine analogues
Cpd	Stereochemistry	NR^1R^2	Ki (nM)	pKi
24	(1S,2R)		3.3	8.481
25	(1S,2R)		3.7	8.432
26	(1S,2R)		1.8	8.745
27	(1S,2R)		17	7.770
28	(1S,2R)		7.3	8.137
29	(1S,2R)		7.6	8.119
30	(1S,2R)		5.8	8.237
31	(1S,2R)		1.9	8.721
32	(1S,2R)		150	6.824
33	(1S,2R)		140	6.854
Test set
N-benzyl amide analogues
Cpd	Stereochemistry	Ar	Ki (nM)	pKi
34	(1S,2R)	C6H5	5.9	8.229
Aromatic ring on cyclopropane core
Compound	Stereochemistry	Ar	Ki (nM)	pKi
35	(1S,2R)	4-Me-C6H4	37	7.432
36	(1S,2R)	4-OMe-C6H4	350	6.456
37	(1S,2R)	3-CF3-C6H4	210	6.678
38	(1S,2R)	3,5-Cl2-C6H3	90	7.046
Phenyl piperidine analogues
Cpd	Stereochemistry	R	Ki (nM)	pKi
39	(1S,2R)		2.8	8.533
Piperazine analogues
Cpd	Stereochemistry	NR^1R^2	Ki (nM)	pKi
40	(1S,2R)		12	7.921

Ar, aromatic rings.

The 3D molecular structures were constructed with Marvin Sketch17, which was also used to calculate molecular and structural descriptors. Each structure was further refined using the OMEGA 2 software18^,19, which employs the Merck molecular force field20 for conformational analysis by systematic search in order to generate 50 minimum energy conformers. The molecular alignment, a necessary step for 3D-QSAR calculations of molecular interaction fields (MIFs) and similarity indices, was performed using ROCs software (OpenEye Scientific Software, Inc., Santa Fe, NM) (Figure 1)21^,22. In this step, we considered the three most active compounds as a template to the molecular alignment which was performed using a chemical similarity score (linear combination of Tanimoto indices for 3D shape and chemical features, Figure 1A).

Figure 1. (A) Molecular alignment of the three most active compounds. Rings 1, 2, 3 and the carbonyl oxygen (acceptor) were used as superimposition sites as well as its 3D shapes. (B) Molecular alignment of all 40 studied compounds. (C) Chemical feature sites were employed as the template for the molecular alignment (distances are shown in Å).

After the generation of conformers by molecular alignment, partial atomic charges were calculated employing the PM3 semi-empiric method23^,24 with MOPAC software. Polar surface area (PSA) and charged partial surface area (CPSA) descriptors, atoms, fragment and ring counts, ClogP, molecular volume, and molecular surface area were calculated with SYBYL-X 2.1.125; molecular orbital energies (HOMO and LUMO) were calculated using the Hartree–Fock/6-31G(d,p) method26 with Gaussian 03W27; ClogD at pH 1.5, 5.0, 6.5, and 7.4, molecular refractivity, molecular polarizability, and topological indices were calculated with Marvin Beans17. These descriptors made up the matrix X and were used in the hierarchical cluster analysis (HCA) and principal component analysis (PCA).

Data set division

The test set should be representative of the total data set with respect to the molecular structure and biological activity. To this end, HCA and PCA were performed using the Pirouette 3.11 software (Infometrix Inc., Bothell, WA)28 in order to explore similarities between the samples (compounds) and guide to the selection of molecules for the test set. HCA considers the complete-linkage clustering method and Euclidean distance to calculate similarities between samples and variables. Equation (1)(1) was used to calculate the distances between the samples: (1) The multivariate distance d_kl between the sample vectors, k and l, is calculated taking into account the numeric differences in all the m variables. M represents the order of the Euclidian distance (M = 2). Once all intersample distances are calculated, the values are organized in a matrix made of similarity indices (similarity _kl = 1 – d_kl/d_max, where d_max is the maximum distance in the data set). The similarity indices span from zero (dissimilar samples) to one (identical samples), and the larger the similarity index value is, the more similar the species are29. The results can be seen as dendrograms (see “Results and discussion” section) and highlights samples agglomeration into clusters. Ideally, one should choose few compounds from each cluster to comprise the test set.

In PCA, the variables are rearranged in principal components (PCs) in such a way that each PC accounts for most of the variability of the data set. The first PC (PC1) is defined taking into account the maximum variability of the data set. Each succeeding PC is taken orthogonally to PC1 and describes the highest possible remaining variability.

After data transformation, just few PCs are necessary to explain the most relevant information contained in the whole data set. The original data matrix X (I × J, where I corresponds to the molecules and J is the number of descriptors), is decomposed in two new matrices, T and L, in such a way that X = TL^T. T is the scores matrix and represents the position (classification) of the compounds on system defined by the PCs. The scores matrix holds important information about the relationships between samples, their similarities and differences, as well as how samples are distributed across the PCs. L is the loadings matrix and describes the weight by which variable should be multiplied to obtain the each component score. In other words, L specifies how the PCs are obtained on the basis of the original variables30^,31. In this study, PCA was run up to eight factors.

QSAR models

All QSAR models were constructed with SYBYL-X 2.1.1 package (Tripos, St. Louis, MO)24 employing partial least squares (PLS) regression technique. Various HQSAR, CoMFA, and CoMSIA models were initially constructed and selected according to their robustness, measured as the leave-one-out (LOO) internal validation regression coefficient (q²) and predictive power, measured through external validation coefficients and its error ( and RMSEP). Next, leave-N-out tests were conducted in triplicate in each of the best HQSAR, CoMFA, and CoMSIA models, leaving 5–30 cross-validation groups out of training set in internal validation. This approach was adopted to verify how the robustness of each model varied according to the size and the composition of the training set. Finally, external validation was carried out using a test set of seven compounds, which were not included in QSAR model development. The predictive ability of the models was evaluated by the residual of the predicted pK_i values and the external validation coefficient ()32. The best-validated QSAR models were used to generate contribution and contour maps for the data set compounds, as well as mutually corroborating physicochemical information. All validation coefficients were evaluated according to acceptable reference values9^,11^,32. Our strategy was based on many works which combine different QSAR methods in the drug design33–38.

HQSAR

HQSAR is a 2D-QSAR method that employs molecular fingerprints to generate models correlating activity with structural features39. Molecules are broken down into fragments which are used to generate molecular holograms related to biological data. First, we fixed the fragment size (4–7 atoms), used the 12 default series of hologram length and evaluated the influence of fragment distinction [atomic type (A), bonds (B), connections (C), hydrogen atoms (H), chirality (Ch), and H-bond donor/acceptor atoms (DA)]. Then, we used the best obtained models, fixed the fragment distinction and evaluated the influence of fragment size in models. First, HQSAR models were generated with the training set ranging the parameters that affect hologram generation. The fragment distinction combinations used were the following: A/B, A/C, A/B/C, A/B/H, A/C/DA, A/C/Ch, A/B/C/H, A/C/H/Ch, A/C/Ch/DA, and A/B/C/H/Ch. In all combinations, we tested all options of hologram lengths (53, 59, 61, 71, 83, 97, 151, 199, 257, 307, 353, and 401). For the best models, the fragment size was varied in length (2–5, 3–6, 5–8, and 6–9) to improve statistical significance.

CoMFA and CoMSIA

CoMFA and CoMSIA descriptors were calculated after molecular alignment, described in Biological data, molecular alignment and descriptors calculations section. In CoMFA, electrostatic and steric energy field descriptors were calculated in each grid point (2 Å) of a cubic box using a sp³ hybridized carbon probe with a charge of + 1.0 and 1.52 Å van der Waals radius. The interactions were calculated with the Tripos force field40 with an energy cut-off of 30 kcal/mol. CoMSIA similarity indices descriptors (steric, electrostatic, hydrophobic, H-bond donor, and H-bond acceptors) were obtained for each point of the grid box using the same probe as in the CoMFA analysis. In order to obtain better CoMFA and CoMSIA models, we employed focus strategy which consists of combining the two follow variations: (i) to focus of fields varying the distance of each point at 3D grid from 1 to 4 Å multiplying the original value by a distance factor (d) and (ii) to apply a standard deviation weight factor (w) from 0.3 to 0.7.

Results and discussion

Variable selection and data set division

After descriptor calculation, the initial X matrix contained 63 descriptors (independent variables) plus the pK_i values (dependent variable) and 40 compounds. The list of all 63 calculated descriptors is available at Supplementary Table 1. Visual inspection of the plots of each independent variable versus pK_i was used as a first filter to exclude invariant descriptors, since they do not contribute to explain the variability of the biological activity. Next, we employed HCA recurrently to identify pairs of descriptors with over 85% of correlation between themselves. The decision of what descriptor should be included in the final model was made assisted by the loading values. We decided to exclude variables with low absolute loading values, since they contribute to a lesser degree to the PCA axis formation. Moreover, variables with similar absolute loading values but opposite signals were not considered because they do not significantly contribute to the new PCA axes. We selected the X variables in order to represent the descriptors with the same nature of the employed ones in HQSAR (constitutional), CoMFA (steric and electrostatic), and CoMSIA (electrostatic and hydrophobic) model construction.

After the elimination of redundant information, the final matrix contained eight variables. The first two PCs were sufficient to explain almost 88% of the original data variance. The scores plot shows samples (compounds) disposition and values on PC1 and PC2, separated in dashed circles. The loading chart shows the contribution of each descriptor for each PC (Figure 2). Steric and topological descriptors (Harary_index, van der Waals volume, and Surface_Area) have greater absolute loading values on PC1. Descriptors related to partial positive charge distribution (WPSA1 and WPSA3) have greater contribution on PC1, while descriptors related to partial negative charge distribution (PNSA3, FNSA3 and RNCS) have greater contribution on PC2.

Figure 2. PCA findings for the studied data set. (A) Factor selection, (B) scores plot, and (C) loadings values.

HCA findings further clarify and support PCA pattern separation, seen as a dendrogram. At 50% of similarity index (solid line in Supplementary Figure 1), three clusters were formed. Molecules from each group were selected based on the shared similarity index (S_i) with other molecules. This way, we ensured that the test set comprised representatives of the total data set. The following compounds were chosen: 34 (94% of S_i with 21); 35 (97% of S_i with 4); 36 (98% S_i with 10); 37 (95% S_i with 5); 38 (98% S_i with 12 and 2); 39 (90% S_i with 18); and 40 (94% S_i with 27). Moreover, the test set pK_i values also represent the distribution of the original data, ranging from 5.4 (6) to 7.5 (25).

All test set has a similar distribution according to its pK_i range values in comparison with training and complete set (Figure 3). The PCA and HCA findings support that all test set compounds have physicochemical similarity to train set ones at the two main found clusters. These results also corroborates the PC1 versus PC2 plots of HQSAR, CoMFA and CoMSIA models (Supplementary Figure 2) which indicates that all test set compounds were inside the chemical spaces employed to QSAR models construction and, therefore, test set was suitable to external validation performing.

Figure 3. Distribution of test, training, and complete sets compounds according to pK_i range.

QSAR models

From the molecular alignment (Figure 1B), we can note that all compounds are within the spatial restrictions of the chemical features. In particular, the cyclopropane ring at the center of all compounds is well aligned. Due to the high flexibility and stereochemistry of some compounds, it is expected that some groups seem to be unaligned. Figure 1(C) displays the distances in Å between each chemical feature site.

HQSAR model

After screening all constructed models with different fragment distinction, the three best models were 29, 41, and 42, according to the statistical results obtained from PLS analysis (Table 2). Model 29 resulted from the combination of A, B, and H fragment distinctions; models 41 and 42 resulted from the combination of A, B, C, and H. We observed that H atoms are important to describe NK3 inhibition since the addition this parameter increases q², as it is observed in models 1 (A/B, q² = 0.596) and 10 (A/B/H, q² = 0.610); 2 (A/B/C, q² = 0.566), and 3 (A/B/C/H, q² = 0.587); and 8 (A/C/Ch, q² = 0.594) and 7 (A/C/H/Ch, q² = 0.611).

Table 2. Statistical parameters analysis of the three bests HQSAR models using fragment size default (4–7).

Model	Distinction	Fragment size	HL	N	q^2	SEP	r^2	SEE		RMSEP
29	A/B/H	1–4	61	6	0.638	0.487	0.898	0.259	0.416	0.353
41	A/B/C/H	1–4	71	6	0.659	0.472	0.905	0.250	0.500	0.437
42	A/B/C/H	2–5	151	6	0.644	0.483	0.910	0.243	0.592	0.416

HL, hologram length; N, optimal number of components; q², LOO cross-validated correlation coefficient; SEP, standard error of prediction; r², non-cross-validated correlation coefficient; SEE, standard error of estimate; , external validation correlation coefficient; RMSEP, root mean squared error of external prediction.

The fragment size controls the minimum and maximum lengths of each fragments encoded in the hologram fingerprint39. We varied the size of the three best models to improve statistical parameters (Supplementary Table 3). The best model is shaded in light gray. Model 42 used fragment size 2–5 and showed acceptable statistical parameters (q²= 0.644 and r²= 0.910), indicating internal consistency. Moreover, the difference between r² and q² is less than 0.30 suggesting that the data are not overfitted. Our model shows that atom types, bonds, connections, hydrogen atoms, and chirality are essential molecular features for NK₃ inhibition.

For model 42, we carried out LNO tests on six different cross-validation groups, in triplicate (Supplementary Figure 3). The results indicate that all LNO cross-validation runs achieved q² values superior to 0.60, which indicates the robustness of our model.

CoMFA model

CoMFA models were generated for the training set based on the molecular alignment. After varying the grid distance and the weight factor, the model with the highest q² value was selected. The best model was generated using a grid distance of 2 Å and a weight factor of 0.50. A statistical comparison between the best CoMFA models is shown in Table 3. Other generated CoMFA models are shown in Supplementary Table 4. Field contribution percentages of the best model show that steric and electrostatic parameters are important to define molecular activity, the latter being more significant (64.1%).

Table 3. Statistical data of the best CoMFA models.

	No focus	Focus1	Focus2	Focus3
w	–	0.30	0.30	0.50
d	–	1.00	1.50	1.00
q^2	0.632	0.758	0.695	0.810
SEV	0.491	0.397	0.442	0.351
N	6	6	6	6
r^2	0.946	0.944	0.925	0.929
SEE	0.188	0.192	0.219	0.216
S	0.366	0.411	0.434	0.359
E	0.634	0.589	0.566	0.641
	–	0.961	0.919	0.969
RMSEP	–	0.567	0.463	0.569

w, standard deviation weight factor; d, distance factor; q², LOO cross-validated correlation coefficient; SEV, standard error of validation; N, optimal number of components; r², non-cross-validated correlation coefficient; SEE, standard error of estimate; Field contribution: S, steric; E, electrostatic; , external validation correlation coefficient; RMSEP, root mean squared error of external prediction.

LOO validation showed that our CoMFA model has internal consistency, with q² higher than 0.70 (0.810) and SEV lower than 0.50 (0.351). LNO internal validation was performed with six different group sizes, varying from 5 to 30. All the q² values were higher than 0.60 (Supplementary Figure 3), emphasizing the consistency of the chosen CoMFA model.

CoMSIA model

CoMSIA model was also generated after molecular alignment using training set for all possible combinations of the five similarity fields (steric, electrostatic, hydrophobic, H-bond donor, and H-bond acceptor). After varying the grid distance and the application of weight factors of the two best non-focused models, the model with the highest q², and lowest RMSEP was selected. The best final model had contributions of steric, electrostatic, hydrophobic, H-bond donor, and H-bond acceptor similarity fields and was generated with a grid distance of 2 Å (d = 1.0) and a standard error weight factor of 0.3. It is important to note that the model generated with a grid distance of 1 Å is statistically equivalent to the chosen model; nonetheless our choice was based on the lowest number of employed grid points. The statistical results of the two best models and the final CoMSIA model are shown in Table 4. The best CoMSIA models obtained with the best two field distinctions and also all models constructed with focus strategy are shown in Supplementary Tables 5 and 6.

Table 4. Statistical data of the best CoMSIA models.

	E/H/D/A No focus	E/H/D/A w = 0.3; d = 1.0	S/E/H/D/A w = 0.3; d = 1.0	S/E/H/D/A w = 0.5; d = 1.0
q^2	0.639	0.679	0.691	0.684
SEV	0.486	0.442	0.426	0.439
N	6	4	3	4
R^2	0.903	0.907	0.911	0.908
SEE	0.252	0.246	0.242	0.246
S	–	–	0.089	0.090
E	0.349	0.340	0.317	0.314
H	0.248	0.222	0.211	0.200
D	0.222	0.214	0.174	0.177
A	0.180	0.223	0.210	0.219
	0.868	0.808	0.781	0.775
RMSEP	0.413	0.643	0.406	0.403

w, standard deviation weight factor; d, distance factor; q², cross-validation correlation coefficient; SEV, standard error of validation; N, optimal number of components; r², noncross-validation correlation coefficient; SEE, standard error of estimation. Field contribution: E, electrostatic; H, hydrophobic; , external validation correlation coefficient; RMSEP, root mean squared error of external prediction.

On one hand, field contribution percentages of the selected models demonstrate that the steric differences between the molecules are less important to describe the activity. On the other hand, electrostatic characteristics are the most significant feature to describe the molecular activity of our CoMSIA model, since this field weight is 31.7%.

The CoMSIA model LOO validation showed q² of 0.691 and SEV of 0.426, indicating its acceptancy. The LNO internal validation was performed with six different group sizes, varying from 5 to 30. All the q² values were also higher than 0.60 (Supplementary Figure 3), reflecting the consistency of the obtained CoMSIA model.

QSAR models external validation

The predictive power of the three best QSAR models was tested by predicting the pK_i values of the test set compounds, which were not included during the construction of the models (Table 5). The plot in Figure 4 shows experimental versus predicted activities.

Figure 4. Plots of predicted versus experimental pK_i values HQSAR (A), CoMFA (B), and CoMSIA (C) best models.

Table 5. Experimental (Exp.) and predicted (Pred.) biological activities (pK_i) for the test set (Cpd).

		HQSAR		CoMFA		CoMSIA
Cpd	Exp. pKi	Pred. pKi	Residual	Pred. pKi	Residual	Pred. pKi	Residual
34	8.229	7.777	−0.452	8.094	0.135	8.266	−0.037
35	7.432	7.561	0.129	7.583	0.151	7.353	0.079
36	6.456	7.052	0.596	6.997	0.541	7.199	−0.743
37	6.678	5.985	−0.693	7.137	0.459	7.455	−0.777
38	7.046	7.955	0.909	7.254	0.208	7.911	−0.865
39	8.553	8.591	0.038	8.516	0.037	8.483	0.070
40	7.921	8.495	0.574	8.149	0.228	8.194	−0.273

Residual values are expressed in the last column for each model.

The external validation of the HQSAR model showed good agreement between experimental and predicted pK_i values (Table 5 and Figure 4A). The HQSAR model obtained is robust and reliable, thus can be used to predict biological activity of novel compounds within this structural class. Both CoMFA and CoMSIA models were able to predict the activity of the test compounds with accuracy, as can be seen on Figure 4(B) and (C). Furthermore, the residual values were all lower than 1.00 (Table 5), reinforcing the model reliability. All constructed QSAR models showed an external validation coefficient () values higher than 0.50, which is acceptable by literature32. The calculated values of HQSAR, CoMFA, and CoMSIA models were equals to 0.592, 0.969, and 0.781, respectively.

In order to verify the applicability domain of our models, we calculated the missing information and also plotted the scores values of first and second principal components of all models. From the first analysis, the missing information of test set compounds did not exceed 0%, 7%, and 3% for HQSAR, CoMFA, and CoMSIA models, respectively. These results indicate that the test set compounds were able to generate around 100% of descriptors used to construct the models. The principal component plots showed that all predicted activity of the test set is inside the range of the training set, suggesting that our models respect the applicability domain.

QSAR models physicochemical interpretation

HQSAR also provides information about what molecular fragments are related to the studied biological activity. This is shown as contribution maps that use color coding to represent the atomic contributions to the activity. Colors reflecting low contribution to the activity are at the red end of the spectrum, while colors reflecting favorable contribution are at the green end41.

The contribution map for compound 19, a phenyl piperidine analogue which is the most potent antagonist, is shown in Figure 5. The piperidine ring and its nitrogen atom, and both cyclopropane and phenyl ring attached to it contribute positively to the activity. Furthermore, other atoms such as the methylenic hydrogen adjacent to the carbamate group and the carbon of the aromatic ring connected to N-benzyl amide group also contribute positively to the activity.

Figure 5. HQSAR contribution maps for the most potent antagonist (A) and the least active antagonist (B). Green and yellow correspond to positive contributions to activity. Red and orange correspond to negative contributions to activity. The colored maps are available in online version.

The least active antagonist, compound 6, presents an aromatic ring on the cyclopropane core (Figure 5). In this case, the aromatic ring connected to the piperidine ring has negative contribution to activity, which was not observed in compound 19. We concluded that substituents on the aromatic ring connected to the cyclopropane ring influence the contribution to the activity of this ring connected to the piperidine moiety.

The results of CoMFA model can be visualized as stereochemical and electrostatic contour maps. For better interpretation, these maps are showed for the most (19) and the least active (6) compounds (Figure 6). Stereochemical contours are depicted in green (favored) and yellow (disfavored), while electrostatic fields are shown in red (negative charges enhancing the activity) and blue (positive charges enhancing the activity).

Figure 6. Energy contribution maps for the most (A) and the least active antagonists (B) generated for the best CoMFA model.

As the steric contour maps show, bulky substituents in the aromatic ring adjacent to the cyclopropane core favor the activity. This property explains the great antagonist activity of compound 19, since it has two chlorine atoms. On the other hand, the least active compound has a trifluormethyl substituent, which has smaller volume. Moreover, compound 19 has a chlorine atom as substituent in the meta position, which increases the volume in a neighboring region, probably favoring interaction with the receptor. Moreover, the methyl group attached to the carbamate is located near a green polyhedral which is not present in compound 6, indicating a steric interaction in this region.

We observed that the positive charge adjacent to the amide nitrogen favors the activity, which indicates the importance of a positively charged group in this position. These maps also show that the negatively charged area near the carbonyl has a great significance to the activity. The comparison between the two compounds highlights the advantage of a carbamate group instead of an amide moiety to increase the electron density in that region, favoring the activity.

Similar to the CoMFA analysis, CoMSIA results can also be interpreted with the aid of contour maps generated for the best model. The steric and electrostatic similarity fields are shown in Figure 7 and the hydrophobic, H-bond donor and H-bond acceptor fields are shown in Figure 8. The steric contour map generated with the CoMSIA model showed that both carbamate group and para substituents of the phenyl attached to the cyclopropane ring have important role in the antagonist activity.

Figure 7. Steric and electrostatic contribution maps for the most (A) and the least active antagonists (B) generated for the best CoMSIA model.

Figure 8. Hydrophobic, H-bond donor, and H-bond acceptor contribution maps for the most (A) and the least active antagonists (B) generated for the best CoMSIA model.

CoMSIA electrostatic contour maps are in agreement CoMFA results. Besides the importance of the negative partial charge around the carbonyl and the positive charge on the nitrogen of the carbamate group, this model also shows the role of a substituted aromatic ring. As said before, great electron density in the aromatic ring can improve compounds’ antagonism.

The hydrophobic contour map calculated for the CoMSIA best model shows a polar favorable region indicating that highly polar groups, such as the carbamate in compound 19, favor the activity. The aromatic substituent also influences the hydrophobicity and, therefore, the molecule efficacy. This result shows that, besides electrostatic contribution, the two chlorines can also favor the antagonism, since they have greater hydrophobic contribution in comparison with a trifluormethyl substituent42^,43. Juhl et al.4 have also observed that molecules with hydrophobic substituents are potent antagonists, while compounds with polar moieties are less active.

H-bond acceptor and donor contour maps show three favorable regions (Figure 8). An H-bond acceptor in the vicinity of the carbamate methoxy region, existent only in compound 19, contributes positively to the activity, as well as the carbonyl of the amide attached to the cyclopropane core. Similarly, the amide nitrogen might play an important role on the interaction with the receptor, since it favors activity.

All QSAR models showed convergent physicochemical evidences of the possible interactions between NK₃ antagonists and its biological targets. Due the statistical consistency of constructed models, this information (summarized in Figure 9) could be employed in the design of new NK₃ antagonists.

Figure 9. Structure–activity relationship obtained from all QSAR models.

Conclusions

HQSAR, CoMFA, and CoMSIA models to predict NK₃ receptor antagonism were developed and validated. All constructed models were robust (q²= 0.644, 0.810, and 0.691, respectively) and externally predictive ( values were higher than 0.50 and all residual values of pK_i prediction were lower than 1 log unity). Further physicochemical interpretation of models was based on contribution and contour maps. From these maps, we could infer that the central cyclopropane ring and the piperidine ring contribute positively to the inhibitory activity of the compounds. The di-chloro-phenyl ring attached to cyclopropane ring also contributes positively to biological activity due its high volume and hydrophobicity. The carbamate attached to piperidine ring has electrostatic contributions (CoMFA) and also contributes to H-bond donor and acceptor favored CoMSIA map. Last but not least, the carbonyl of the amide attached to the cyclopropane core also provides an H-bond favorable site to the activity. These insights show that some cores and atoms have different effects on the activity of the inhibitors and are useful for the design of new NK₃ receptor antagonists.

Supplementary material available online Supplementary Figures 1-3 and Tables 1-6.

Acknowledgements

The authors thank NAP – Univ São Paulo Center for Innovation in Diagnostics & Therapeutics (IDx&T). We also thank OpenEye Scientific Software Inc. for providing Academic License of the software (ROCs and OMEGA 2). We thank alNAP – Univ São Paulo Center for Innovation in Diagnostics & Therapeutics (IDx&T). The authors also acknowledge Prof. Dr. Marcelo S. Castilho (Federal University of Bahia/Brazil) for having provided access to Pirouette 3.11 software.

Declaration of interest

The authors would like to thank FAPESP (M. C. P. 2013/15947-3), CAPES, and CNPq (G. H. G. T. 477002/2013-7) for financial funding.