Three-Dimensional Quantitative Structure-Activity Relationships (3D-QSAR) on a Series of Piperazine-Carboxamides Fatty Acid Amide Hydrolase (FAAH) Inhibitors as a Useful Tool for the Design of New Cannabinoid Ligands

Fatty Acid Amide Hydrolase (FAAH) is one of the main enzymes responsible for endocannabinoid metabolism. Inhibition of FAAH increases endogenous levels of fatty acid ethanolamides such as anandamide (AEA) and thus consitutes an indirect strategy that can be used to modulate endocannabinoid tone. In the present work, we present a three-dimensional quantitative structure-activity relationships/comparative molecular similarity indices analysis (3D-QSAR/CoMSIA) study on a series of 90 reported irreversible inhibitors of FAAH sharing a piperazine-carboxamide scaffold. The model obtained was extensively validated (q2 = 0.734; r2 = 0.966; r2m = 0.723). Finally, based on the information derived from the contour maps we designed a series of 10 new compounds with high predicted FAAH inhibition (predicted pIC50 of the best-proposed compounds = 12.196; 12.416).


Introduction
The endocannabinoid system (ECS) remains a highly relevant topic in the scientific community as it is involved in several regulatory actions and pathophysiological conditions [1]. Current available knowledge suggests that the ECS is constituted by the cannabinoid receptors, type 1 and 2, the main endogenous ligands anandamide (AEA) and 2-arachidonyl glycerol (2-AG) as well as the enzymes that participate in their biosynthesis (N-acyl phosphatidylethanolamine phospholipase D, or NAPE-PLD) and degradation (Fatty Acid Amide Hydrolase or FAAH and Monoacylglycerol Lipase, or MAGL) [2].
FAAH is an integral membrane protein of~60 kDa (579 amino acids) that belongs to the amidase family of enzymes [3]. It exists as a dimer in its membrane-associated form [3] and is highly expressed Some quantitative structure-activity relationships (QSAR) studies have been performed on the cannabinoid system (CB1 and CB2 receptors) [8][9][10][11][12][13], however, few structure-activity relationship studies have been performed on FAAH inhibitors, and mostly with carbamate-type structures. Dainese et al. calculated theoretical molecular descriptors in a series of naturally occurring FAAH inhibitors [14]. Käsnänen et al. reported the synthesis and 3D-QSAR studies of carbamate inhibitors [15]. Mor et al. constructed 2D-QSAR equations that could explain the inhibition activity of biphenyl-alkylcarbamates. [6]. Vacondio et al. developed structure-property relationships to explain the hydrolytic stability of carbamate inhibitors [16]. Han et al. reported a comparative molecular field analysis (CoMFA) study on a series of oleoylethanolamide structure inhibitors [17]. To date, there are no 3D-QSAR studies of irreversible inhibitors with the piperazine-carboxamides structure. This type of general structure was shown to have good physical and pharmacokinetic properties and has been reported to be capable of elevating plasma concentrations of AEA, AEP, and AEO in rats [18]. For this reason, the formulation of a QSAR model for the design and prediction of FAAH inhibitor activity based on this structural moiety is significant from a pharmacological point of view.
In the present work, three-dimensional quantitative structure-activity relationships (3D-QSAR) studies based on comparative molecular similarity indices analysis (CoMSIA) were carried out on a set of various reported urea-based FAAH inhibitors. The aim of our 3D-QSAR is to derive useful binding information in order to guide the design of future FAAH inhibitors. The importance of steric, electrostatic and hydrogen-bond characteristics can be analyzed by aligning similar analogues based on key pharmacophoric features [19]. Knowledge of binding requirements can then be used to derive predictive 3D-QSAR models that can, in turn, aid in the design of new inhibitors.

Statistical Results
The statistical results for CoMSIA are presented in Table 1. All possible field combinations were tested for both CoMFA and CoMSIA. In the case of CoMFA, no combination was statistically significant. The CoMSIA models with the highest q 2 values were those that considered the field combinations SEDA, EDA, EHDA, and SEHDA. The SEDA and EHDA models presented a donor hydrogen-bond contribution of 0.099 and 0.093 respectively. While in the EDA model, the H-bond donor contribution was 0.111 versus 0.889 of the Electrostatic and H-bond Acceptor contributions. The imbalance in the field contribution of these models made us discard them. The final selected model SEHDA, presents a good balance between the field contributions, a high value of q 2 (0.734) and r 2 ncv (0.937) and higher value of F (138.36). Extensive additional validation was carried out to test the predictive quality of this model.  Table 2 presents a summary of the external validation parameters of the CoMSIA model. The model has a high value of r 2 (0.966), which indicates satisfactory external predictive power. Nevertheless, as described by Golbraikh and Tropsha, high values of q 2 and r 2 (conditions 1 and 2) are not sufficient for model validation. For a reliable predictive capacity, the line of experimental versus predicted activity should be as close as possible to the line y = x. This is observed when the conditions (3a or 3b), (4a or 4b) and (5a or 5b) shown in Table 2 are fulfilled. Condition 6 or r 2 m metrics, measures the agreement between the observed and predicted activity. Our derived CoMSIA model satisfied all the conditions for internal and external validation. In addition, we performed the calculation of several external validation descriptors (conditions 7 to 12) [20][21][22]. In all cases, our model passed the validation tests.
Furthermore, the Y-randomization test [23] (Table 3) was applied to assess the robustness of the model (see Table S1 of the Supplementary Material) as previously described by Lorca et al. [24]. The obtained QSAR models show low q 2 and r 2 ncv values ( Table 3). Close to value of ∆r 2 m <0.2 0.056 q 2 = the square of the LOO cross-validation (CV) coefficient; r 2 is the regression coefficient for the test set exclusively; r 0 2 and k are the correlation coefficient between the actual and predicted activities for test set and the respective slope of regression; and r 0 ' 2 and k' are the correlation coefficient between the predicted and actual activities for test set and the respective slope of regression. r 2 m was defined in equation 5. Parameters 8-12 are defined in the methods section. The values of experimental activity, predicted activity, and the residual values for the best CoMSIA model are shown in Table 4. All the compounds showed low residual values and deviations of the predicted activity greater than a logarithmic unit were not observed. Figure 2A shows a plot of experimental versus predicted activity and the data distribution is close to the y = x line. The model shows a good balance in terms of predictive capacity. Forty-two compounds showed negative residual values and 48 presented positive deviations ( Figure 2B). The residual range was from −0.82 to 0.89. As shown in Figure 2C the CoMSIA model shows a satisfactory predictive capability throughout the whole data set (training and test set) as well as a good predictive power for both, less active (1, 6 and 7) and most active compounds (65, 66, and 67).

Applicability Domain
In this work we used the method developed by Roy et al. [25] for determination of applicability domain (AD) as previously described by Lorca et al. [24].
The calculation was carried out using the free application available on the author's page and all compounds were found to be within the domain of applicability, except compounds 83 and 87. These two compounds are the only ones bearing an imidazopyridine or imidazopyrimidine ring and the only difference between them is the position in which the heterocycle is connected to the urea moiety. For this reason, compounds with these heterocyclic systems connected to the urea moiety were not proposed as new molecules.
In summary, the CoMSIA model generated here presents good internal and external validation parameters (q 2 = 0.734; r 2 = 0.966), and meets the validation criteria of Tropsha and Roy (r 2 m = 0.723). All the molecules studied are within the applicability domain (except compounds 83 and 87) and the model was validated by the Y-randomization test. Therefore, reliable information can be extracted from the contour maps as discussed in the next section.

Contour Maps Analysis
The result of a 3D-QSAR study can be visualized graphically unlike a traditional 2D-QSAR equation. Contour maps represented by colored polyhedrons can be seen around the molecules. The maps obtained in our study correspond to the steric, electrostatic, hydrophobic, H-bond donor and H-bond acceptor contour maps. Regions where a molecular property is favorable or unfavorable are indicated by different colored polyhedrons. Figure 3 presents the different maps around the most active compound (66, pIC 50 = 10.602; on the left) and the least active compound (6, pIC 50 = 5.176; on the right).

Steric Contour Map
The steric contour map shows a large yellow polyhedron close to the pyridazine ring of the most active compound indicating that bulky substituents in this region should be avoided in order to favor biological activity ( Figure 3A,B). Alternatively, a smaller five-membered ring could replace the pyridazine ring following the same steric requirement. This can be seen in the proposed molecules (Table 5) where the analog 2x shows a considerable increase in predicted activity when the pyridazine ring is replaced by a pyrazole ring. This relation is further supported by compounds 54 and 55 (Table 6) which bear phenyl substituents in the corresponding position and show low activity consistent with their bulkier nature.       On the other side, the yellow polyhedron that surrounds the pyrimidine ring indicates that reducing the size of this ring or replacing it with a smaller linker while maintaining the electronic properties would be beneficial for activity. Additionally, the green polyhedrons around the ortho and meta positions of the disubstituted benzene ring indicate that bulky substituents in these positions can be favorable for activity. This can also be seen in Table 5 where substitution with a methyl group considerably increases the predicted activity in all the proposed molecules. For the least active compounds ( Figure 3B), the steric factor by itself does not seem to explain the lower activity values.

Electrostatic Contour Map
Regarding the electrostatic contour maps, the red polyhedrons around the pyridazine ring ( Figure 3C) highlight the importance of nitrogenated heterocycles that can confer electronegative areas in this region. This may explain why molecules 1, 6 and 7 (Table 6) bearing benzene rings with more homogeneous charge distribution show lower activity. The blue polyhedron inside the pyridazine ring shows that an electropositive center is beneficial for activity, therefore, incorporating electro withdrawing substituents in positions 5 and 6 of the pyridazine ring could increase activity. This electron distribution with an electron rich edge and electron deficient center suggests possible pi-stacking or pi-cation interactions with the target enzyme. Similarly, the expansion of the blue contour at position 4 of the pyridazine ring indicates that electron withdrawing substituents particularly at this position are favorable for activity. In agreement with this, proposed analogues 1x, 9x and 10x ( Table 5) that follow this substitution pattern display high predicted activity. The polarization of the carbon atom directly attached to the electroattractive or electronegative groups nitrile and fluorine lower the electron density right where the blue polyhedron lies.
On the other side, the red contour over the pyrimidinic nitrogen ( Figure 3C) shows the importance of an electronegative atom at this position as is present in the most active molecules (65, 66, 67 and 68 from Table 6) and absent in the least active one (compound 6). Likewise, an electron rich benzene ring is favorable and thus replacing the fluorine atom for an electrodonating group would be recommended. Accordingly, the proposed molecules 9x and 10x substituted with electrodonating methylene groups show the highest predicted activities. In Figure 3D the red contour over the electron deficient nitrile carbon atom may explain the lower activity of this analog.

Hydrophobic Contour Map
The hydrophobic contour maps ( Figure 3E,F) show a gray polyhedron over the pyridazine nitrogen atoms of the most active molecules, indicating that incorporating hydrophilic atoms in this region can favor activity. The yellow polyhedron over the pyrimidine ring shows that a hydrophobic linker region is important. Additional yellow polyhedrons surrounding the disubstituted benzene ring suggest hydrophobic substituents in the meta and ortho positions can also increase activity. Following both the hydrophobic and the previously mentioned steric requirement all proposed molecules (Table 5) possess an ortho-methyl substituent in the benzene ring.

Donor and Acceptor Contour Maps
In the hydrogen bond donor map ( Figure 3G,H) cyan polyhedrons surrounding the pyridazine ring suggest that incorporation of hydrogen bond donor groups can favor activity. For this reason, the proposed molecule 3x (Table 5) was designed with a hydroxyl group able to form hydrogen bonds. Furthermore, cyan polyhedrons around the urea linker suggest that the urea moiety is involved in a hydrogen bond interaction with the target enzyme.
Finally, the hydrogen bond acceptor map shows red polyhedrons next to position 4 of the pyridazine ring ( Figure 3I,J), position 5 of the pyrimidine ring and over the urea carbonyl. Indicating that hydrogen bond acceptor groups in these positions is unfavorable for activity. Therefore, using different linker groups without hydrogen bond acceptor groups could be advisable in order to design new inhibitors.

Design of New FAAH Inhibitors
Based on the information obtained from the contour maps, we have designed a series of compounds evaluating multiple combinations of fragments. Substituents and functional groups were proposed, taking into consideration the electronic, steric, hydrophobic and hydrogen bonding properties suggested by the contour maps. Table 5 shows the compounds that presented the best predicted inhibitory activity. All proposed molecules have better predicted activity than the most active compound in the series (66, pIC 50 = 10.602). In general, the presence of 6-member rings or fused systems on the left side did not greatly increase activity (1x pIC 50 = 10.889; 4x pIC 50 = 10.990). However, the insertion of a pyrazole ring generated derivatives with a significant increase in the pIC 50 value (best compounds: 9x, pIC 50 = 12.416; 10x, pIC 50 = 12.196). This is because the pyrazole system meets the electronic, hydrophilic and hydrogen bonding requirements suggested by the contour maps. On the right side, the insertion of halogens and alkyl groups slightly increased the activity.
Due to the reported potential toxicity of heterocyclic compounds similar to those presented in Table 5 [26], we conducted a predictive toxicity study using the PreAdmet online platform (https://preadmet.bmdrc.kr/adme/). The calculation showed that only compound 1x presents a high risk of toxicity mediated by affinity to the anti-target hERG.

Molecular Alignment
CoMSIA studies were performed with Sybyl X-1.2 software (1.2, Tripos International, St. Louis, MS, USA) [27] installed in a Windows 10 environment on a PC with an Intel core i7 CPU. Each compound was drawn in ChemDraw and then geometry optimized using MM2 molecular mechanics as implemented in ChemBio3D software (15.1.0, PerkinElmer, Waltham, MA, USA). The structures were further minimized by Tripos force field implemented in Sybyl. MMFF94 charges were assigned to each atom. Among the techniques to perform molecular alignment are: (1) atom-by-atom alignments using a common fragment, (2) rigid alignment that minimizes rms distance, (3) flexible alignments [28] and (4) receptor guided alignments. In the present study, the first two options were carried out (Figure 4). In the first case, the piperazinyl-urea nucleus was chosen as the common scaffold for alignment. In the second case, the alignment was carried out using the Distill rigid protocol, as implemented in Sybyl. The best results in terms of statistical validation (q 2 , N) were for the atom-by-atom alignment, which was chosen as the basis for the development of the models.

CoMSIA Field Calculation
To derive the CoMSIA descriptor fields, the aligned training set molecules were placed in a 3D cubic lattice with grid spacing of 2Å in x, y, and z directions such that the entire set was included in it. The CoMSIA analysis, the standard settings (probe with charge +1.0, radius 1Å, hydrophobicity +1.0, hydrogen-bond donating +1.0, hydrogen bond accepting +1.0 [29]) were used to calculate five different fields: steric, electrostatic, hydrophobic, donor and acceptor. Gaussian-type distance dependence was used to measure the relative attenuation of the field position of each atom in the lattice, and led to much smoother sampling of the fields around the molecules when compared to CoMFA. The default value of 0.3 was set for the attenuation factor α.

Data Set Selection and Inhibitory Activity
CoMSIA studies were performed on a set of 90 piperazinyl urea derivatives reported in the literature [18,[30][31][32][33][34][35][36] (Table 6). The derivatives displayed potent fatty acid amide hydrolase (FAAH) inhibitors activity. The IC 50 values were converted to pIC 50 (-logIC 50 ). The compounds were divided into training (73 compounds, 81%) and test sets (17 compounds, 19%), ensuring that both sets contained structurally diverse compounds with high, medium and low activity, and a uniform distribution to avoid possible problems during the external validation. The distribution of pIC 50 values for the whole set, the training set and the test set is shown in Figure 5. In all three cases the biological activity follows a Gaussian distribution.

Internal Validation and Partial Least Squares (PLS) Analysis
PLS analysis [37] was used to construct a linear correlation between the CoMFA and CoMSIA descriptors (independent variables) and the activity values (dependent variables) as previously described by Lorca et al. [24].

External Validation
The models were subjected to external validation criteria according to the proposed test by Golbraikh and Tropsha [38,39], which considers a QSAR model predictive, if the following conditions are satisfied: 0.85 ≤ 1.15 or 0.85 ≤ k ≤ 1.15 (4) It has been demonstrated [38] that to adequately assess the predictive ability of a QSAR model the above criteria are necessary.
Furthermore, the external predictive power of the developed 3D-QSAR models using the test set was examined by considering r 2 m metrics as shown below [40]: where r 2 and r 2 0 are squared correlation coefficients between the observed and predicted activities of the test set with and without intercept, respectively. For a significant external model validation, the value of r 2 m should be greater than 0.5.
Additionally, the following descriptors were calculated: where, TR = training set, EXT = external prediction set, y i = experimental data values,ŷ i = predicted data values, y = average of the experimental data values,ŷ = average of the predicted data values. Finally, r 2 m is calculated using the experimental values on the ordinate axis, while r 2 m using them on the abscissa.

Applicability Domain (AD) Calculation
The AD was evaluated based on the simple standardization method reported by Roy et al. [25] and as described by Lorca et al. [24].

Conclusions
In this contribution, a 3D-QSAR CoMSIA study was carried out on an extensive database of 90 irreversible inhibitors of the enzyme FAAH with a pyrimidinyl-piperazine-carboxamide general structure. The best model obtained considered all the field contributions, being the electrostatic and hydrogen-bond acceptor properties the ones that contributed most to the activity (30.4% and 33.0% respectively). The model was validated internally (q 2 = 0.734) and externally (r 2 test = 0.966) and was also submitted to Tropsha validation criteria, r 2 m calculation (0.723) and Y-randomization test, passing all tests. The information derived from the contour maps was used to design a series of new compounds that showed promising predicted activities (pIC 50 of the most active compounds = 12.196 and 12.416). The main structure-activity relationships found in this study and summarized in Figure 6 are a useful tool to guide the future design of new FAAH inhibitors. The extensive database used in this study could motivate future work complementing the information obtained from contour maps with QSAR studies by Dragon-based descriptors [41].

Conflicts of Interest:
The authors declare no conflict of interest.

3D-QSAR
Three-dimensional Quantitative Structure-Activity Relationship CoMSIA Comparative Molecular Similarity Indices Analysis FAAH Fatty Acid Amide Hydrolase CB Cannabinoid