Molecular Modeling Studies on 11H-Dibenz[b,e]azepine and Dibenz[b,f][1,4]oxazepine Derivatives as Potent Agonists of the Human TRPA1 Receptor

A computational strategy based on comparative molecular fields analysis (CoMFA) and comparative molecular similarity indices analysis (CoMSIA) was performed on a series of the 11H-dibenz[b,e]azepine and dibenz[b,f][1,4]oxazepine derivatives as potent agonists of the human TRPA1 receptor. The CoMFA and CoMSIA models resulting from a 21 molecule training set gave r2cv values of 0.631 and 0.542 and r2 values of 0.986 and 0.981, respectively. The statistically significant models were validated by a test set of five compounds with predictive r2pred. values of 0.967 and 0.981 for CoMFA and CoMSIA, respectively. A systemic external validation was also performed on the established models. The information obtained from 3D counter maps could facilitate the design of more potent human TRPA1 receptor agonists.


Introduction
The transient receptor potential ankyrin 1 (TRPA1) receptor is a member of the transient receptor potential (TRP) family of the cation-selection channel and the only mammalian member of the TRPA subfamily [1][2]. It is expressed in the dorsal root ganglion, trigeminal ganglion (TG) neurons [3], and non-sensory tissue [4]. It plays an essential role as a biological sensor to irritant chemicals [5] and is implicated in a growing number of diseases, such as bladder disorders [6], inflammatory pain [7], and airway diseases [8]. The activation of TRPA1 by a diversity of chemical agents is widely accepted. It can be activated by many pungent chemicals, including methyl salicylate (from wintergreen oil) [3], isothiocyanates like allylisothiocyanate, the pungent compound in mustard oil (MO), wasabi, and horseradish [3], cinnamaldehyde [3], Δ9-tetrahydrocannabinol (Δ9-THC, the psychoactive compound in marijuana) [9], allicin and diallyl disulphide (from garlic) [10], acrolein (an irritant found in vehicle exhaust fumes and tear gas) [10], and the lacrimators 1-chloroacetophenone (CN), dibenz[b,f]- [1,4]oxazepine (CR) and 2-chlorobenzylidene malononitrile (CS) [11]. Most of the known activating compounds contain reactive, electrophilic chemical groups that react with cysteine residues in the active site of the TRPA1 channel [11]. A considerable amount of investigation has been carried out to develop reversible ligands that target TRPA1 receptor and several agonists have been under pharmacological evaluation. Recently, a series of compounds containing 11H-dibenz[b,e]azepines, and dibenz[b,f ] [1,4]oxazepines that function as extremely potent activators of the human TRPA1 receptor activities were reported by literature [5].
In this paper, molecular modeling studies of these 11H-dibenz[b,e]azepine and dibenz[b,f][1,4]oxazepine derivatives were performed using 3D-QSAR approaches. Thus, 3D-QSAR, including comparative molecular field analysis (CoMFA) and comparative molecular similarity indices analysis (CoMSIA) methods, were performed to predict the agonistic activities of these agonists and provide the regions in space where interactive fields may influence the activity. Furthermore, the 3D-QSAR models were validated by external validation. The constructed models can not only help in understanding the structure-activity relationship of these compounds but also may serve as a useful guide for the design of new agonists with much higher agonistic potencies.

Results and Discussion
A data set of 21 diverse analogues was selected as a training set to derive the conventional CoMFA and CoMSIA models and an additional five compounds (indicated with "*" in Table 2) were used to test the accuracy of these models. Selection of the training set and test set molecules was done randomly. The structures and associated biological activities are shown in Tables 1 and 2.

CoMFA PLS Analysis
The CoMFA model of 11H-dibenz[b,e]azepines and dibenz[b,f ] [1,4]oxazepines derivatives gave a cross-validated correlation coefficient (r 2 cv ) of 0.631 (>0.5) with an optimized component of 6, which indicated that this model may be considerably reliable to predict the pEC 50 values of the compouds of the test set.    The actual and predicted pEC 50 values of the training set and test set by the model were given in Table 2, and the graph of actual activity versus predicted pEC 50 of the training set and test set is illustrated in Figure 1.    Table 3. The actual and predicted pEC 50 values and residual values for training set and test set compounds are given in Table 2. The relationship between actual and predicted pEC 50 of the training set and test set compounds is illustrated in Figure 2.

External Validation Analysis for the CoMFA and CoMSIA
As can be seen from Table 4, the predictive correlation coefficient (r 2 pred ) values based on molecules of the test set were 0.967 and 0.981 for the CoMFA and CoMSIA models, respectively. were 0.002 and −0.007 (<0.1) or the CoMFA and CoMSIA models, respectively. All the data discussed above suggested the fact that both CoMFA and CoMSIA models have not only good estimation abilities, but also robust predictive powers.

CoMFA Contour Maps Analysis
To visualize the information of the derived 3D-QSAR models, CoMFA contour maps were generated by plotting the coefficients from the CoMFA model. These contour maps may help identify important regions where changes in the steric and electrostatic fields were predicted to increase or decrease the activity. As shown in Figure 3, the yellow contour around the C-8 position signified that a bulky substituent at this site would decrease potency. This is consistent with the fact that the analogues modified at C-8 with bulky substituents such as compounds 8 (8-COOMe, pEC 50 = 7.000), 11 (8-Br, pEC 50 = 7.482), 13 (8-CN, pEC 50 = 6.801) and 20 (8-COO i Pr, pEC 50 = 6.030, 10,000-fold decrease in potency compared to compound 10) showed relatively decreased activities. The yellow contour near the C-9 indicated that a small group at this site would be favorable.  3 OMe, pEC 50 = 8.114), and 23 (10-CONEt 2 , pEC 50 = 7.398) also resulted in decreased potencies when the 11-O was replaced by 11-CH 2 . In addition, two small yellow contours around the C-6 and C-4 respectively which were not observed in the CoMSIA steric field indicated that the small groups at these two sites may be beneficial for agonistic activity.

CoMSIA Contour Maps
The steric ( Figure 5) and electrostatic field ( Figure 6) contour maps of the CoMSIA model were almost the same as the CoMFA-steric and electrostatic contours (Figure 3 and Figure 4). However, in the CoMFA model, two more contours appeared: one yellow contour map near the C-4 and a yellow contour map around the C-6 position.
The hydrophobic field is presented in Figure 7, yellow and white contours highlighted areas where hydrophobic and hydrophilic properties were preferred, respectively. The presence of a large yellow contour surrounding the 10-COOMe of the template molecule (compound 10) indicated that hydrophobic substituents may be well tolerated in that region.   In the CoMSIA hydrogen bond fields, the purple contour near the 11 position (A) and CO of the template molecule revealed that hydrogen bond donor groups may decrease the potency. In fact, the -O and -CO at this position acted as hydrogen bond acceptors, and this may explain why compounds 10 (10-COOMe, pEC 50 = 10.301), and 18 (10-CONH 2 , pEC 50 = 10.097) showed relative better activities. A hydrogen bond donor contours features, displayed by cyan contours near the terminal of 10-COOMe revealed that this group acted as a hydrogen bond donor and would be favored over the groups found in compounds such as 18 (10-CONH 2 , pEC 50 = 10.097) and 19 (1-CONH 2 , pEC 50 = 9.959).
The hydrogen bond acceptor field contour map of CoMSIA is shown in Figure 9 using compound 10 as a reference molecule. The magenta and red contours represent favorable and unfavorable hydrogen bond acceptor groups. In the CoMSIA hydrogen bond fields, the magenta contour near the 10-CO and 11-O revealed that hydrogen bond acceptor groups may benefit the potency. These results are supported by the evidence of the high potency of 10 (10-COOMe, pEC 50 = 10.301), and 18 (10-CONH 2 , pEC 50 = 10.097). A red contour around the terminal of 10-COOMe suggested that hydrogen bond acceptor groups may increase the agonistic activity. This may be the reason why compounds 18 (10-CONH 2 , pEC 50 = 10.097) and 19 (1-CONH 2 , pEC 50 = 9.959) showed relative increased activities.

Summary of the Structure-Activity Relationship Based on CoMFA and CoMSIA Models
The structure-activity relationships revealed by 3D-QSAR studies are illustrated in Figure 10. In detail, bulky groups at the C-1 and C-2 positions are favorable; small groups at C-4 and C-6 may be essential for the agonistic potency; small, low electron density, and hydrophilic groups at C-8 position could increase the agonistic activity; small and electron-donating substituents at C-9 may benefit the potency; bulky and hydrogen bond donor groups at C-10 would increase the potency; high electron density, hydrophilic, and hydrogen bond acceptor substituents at the 11 position (A) would be essential for the potency of the agonists.

Data Sets
The twenty-six compounds involved in this study were taken from the literature [5]. The associated activities were reported as EC 50 values towards the human TRPA1 receptor. The EC 50 values were converted into pEC 50 by taking Log(1/EC 50 ). From Table 2, the pEC 50 values for the 26 studied agonists ranged from 6.030 to 10.301.

Molecular Modeling and Database Alignment
Molecular modeling and database alignment were performed by using the molecular modeling package SYBYL 8.1 (Tripos, Inc.) [12]. The 3D structures of all compounds were constructed by using the Sketch Molecule module. Energy minimization of each structure was performed using the SYBYL energy minimizer Tripos force field and Gasteiger-Hückel charge [13][14]. Molecular alignment was considered as one of the most sensitive parameters in 3D-QSAR analysis [15][16]. The quality and the predictive ability of the model are directly dependent on the alignment rule [17]. In this paper, all of the structures were aligned into a lattice box by fitting with (Z)-N-benzylidenebenzenamine group ( Figure 11) as a common structure using compound 10 as a template, which was the most active compound. The aligned superimposed molecules are shown in Figure 12.

CoMFA Modeling
CoMFA is a widely used 3D-QSAR method which relates the biological activity of a series of molecules with their steric and electrostatic fields. The CoMFA descriptor fields were calculated at each lattice with grid spacing of 1 Å and extending to 4 Å units in all three dimensions within defined region [18][19]. The Van Der Waals potentials and Coulombic terms, which represented steric and electrostatic fields, respectively, were calculated by using the standard Tripos force field. In CoMFA method, a sp 3 hybridized carbon atom with a charge of 1e was used as a probe atom, the energy values of the steric and electrostatic fields were truncated at 30 kcal/mol [19][20].

CoMSIA Modeling
The steric, electrostatic, hydrophobic, hydrogen bond donor and hydrogen bond acceptor CoMSIA potential fields were calculated at each lattice intersection of a regularly spaced grid of 1 Å and extending to 4 Å using a probe atom with radius 1.0 Å, +1.0 charge, and hydrophobic and hydrogen bond properties of +1. The attenuation factor was set to the default value of 0.3 [21][22].

PLS Analysis
The partial least-squares (PLS) approach, an extension of multiple regression analysis, was applied to linearly correlate the CoMFA and CoMSIA fields to the pEC 50 values. CoMFA and CoMSIA descriptors were used as the independent variables. Column filtering was used at the default value of 2.0 kcal/mol in the cross-validation part.
The cross-validation analysis was performed using the leave-one-out (LOO) method in which one molecule was omitted from the dataset. The activity of the omitted molecule was then predicted by using the model derived from the rest of the dataset [23]. The leave-one-out (LOO) cross-validation method could check the predictivity of the obtained model and identify the optimum number of components (ONC). Thus the optimum number of components (ONC) was the number of components lead to the highest cross-validated correlated correlation coefficient r 2 (r 2 cv ) [23]. Finally, the CoMFA and CoMSIA models were generated using non-cross-validated PLS analysis with the optimum number of components (ONC) determined by the cross-validation.

Models Validation-Predictive Correlation Co-Efficient (r 2 pred )
The predictive abilities of 3D-QSAR models were validated by predicting the activities of a test set of five compounds which were not included in the training set. These molecules were aligned to the template and their pEC 50 values were predicted by the produced models which were obtained using the training set. The predictive correlation coefficient (r 2 pred ), based on the molecules of test set, was calculated using the following Eq.: In this equation, SD is the sum of the squared deviations between the agonistic activities of the test set and the mean activity of the training molecules and PRESS is the sum of squared deviations between predicted and actual activity values for each molecule in the test set [13,[24][25][26].

Models Validation-External Validation
The previous researches provided the fact that a high cross-validated correlation coefficient, r 2 cv , was the necessary condition for a robust predictive power but not a sufficient condition. The only way to estimate the true predictive ability of the model was external validation. A reliable 3D-QSAR model should have a robust predictive ability, if it is close to the ideal one. 3D-QSAR models were considered acceptable if they satisfy all of the set of criteria for evaluation of predictive ability of QSAR models [27][28][29][30]: r 2 cv > 0.5, r 2 > 0.6, [(r 2 − r 0 2 )/ r 2 ] < 0.1, 0.85 ≤ k ≤ 1.15 and r 2 m > 0.5 In this paper, the CoMFA and CoMSIA models were subjected to rigorous external validation process. We calculated the statistical parameters of the test set according to references [27][28][29][30].

Conclusions
We have employed 3D-QSAR to explore the structure-activity relationship of a series of 11H-dibenz[b,e]azepines and dibenz[b,f] [1,4]oxazepines derivatives as potent TRPA1 agonists. CoMFA and CoMSIA analyses were used to build statistically significant models with good correlative and predictive capability for the activator of the human TRPA1 receptor. These models could be used to predict the potencies of related structures. Based on the excellent results of the external validation, the models established in the present study may be robust and reliable for predicting new derivatives. Moreover, the analysis of contours for the models has provided a clue for the design of new analogues with improved affinity and higher potency.