3 D-QSAR on a series of pyrimidinyl-piperazine-carboxamides based Fatty Acid Amide Hydrolase ( FAAH ) inhibitors as a useful tool to obtain novel endocannabinoid enhancers

Fatty Acid Amide Hydrolase (FAAH) is one of the enzymes responsible of endocannabinoids metabolism. The inhibition of FAAH is a useful and indirect strategy to raise endogenous cannabinoid concentrations, which would be useful for the treatment of various pathological processes in which cannabinoid concentrations are lowered. In the present work, we present an extensive 3D-QSAR/CoMSIA study on a series of 90 irreversible inhibitors of FAAH of pyrimidinyl-piperazine-carboxamide structure. The final model obtained was extensively validated (q2 = 0.734; r2test = 0.966; r2m = 0.723), and based on the information derived from the contour maps we reported a series of 10 new compounds designed as powerful FAAH inhibitors (pIC50 of the bestproposed compounds = 12.196; 12.416).


Introduction
The endocannabinoid system remains a highly relevant topic in the scientific community as it is involved in several pathophysiological conditions and has been associated to many different regulatory actions.[1] The knowledge of this system, constituted by the cannabinoid receptors (CB1 and CB2), the endogenous ligands Anandamide (AEA) and 2-Arachidonyl glycerol (2-AG) and the enzymes associated with its biosynthesis (N-acyl phosphatidylethanolamine phospholipase D, or NAPE-PLD) and degradation (Fatty Acid Amide Hydrolase or FAAH and Monoacylglycerol Lipase, or MAGL), have provided neurochemical tools for the design of new drugs.[2] FAAH was the first enzyme responsible for endocannabinoid hydrolysis to be purified and characterized.FAAH is an integral membrane protein of ~60kDa (579 amino acids) that belongs to the amidase family of enzymes.[3] It is highly expressed in the brain, liver, kidney, and testis, but not in the heart or skeletal muscle.[4] FAAH exists as a dimer in its membrane-associated form, [3] and unlike most other mammalian serine hydrolases, which use a Ser-His-Asp triad for catalysis, FAAH possesses an unusual Ser-Ser-Lys triad.[4] FAAH cleaves and inactivates a broad range of lipid amides in vitro and can tolerate variations in the amide headgroup and lipid acyl chain.Bioactive lipid amide substrates for FAAH include, in addition to anandamide, related N-acylethanolamine (NAE) congeners, the sleep-inducing factor oleamide (cis-9,10-octadecanoamide), and the transient receptor potential (TRP) ion channel agonists the N-acyl taurines (NATs) [4].
Endogenous FAAH substrates such as anandamide serve key regulatory functions in the body and have been implicated in a variety of pathological conditions, including pain, anxiety, depression and vascular hypertension.[5] Therefore, inhibition of the endocannabinoid hydrolases such as FAAH offers a rational therapeutic approach to treat conditions in which a higher endocannabinoid activity would be beneficial.An advantage of this indirect cannabinoid activation over the direct action of agonists could be higher selectivity, as the enzyme inhibition would increase the activity of the endocannabinoid system only at the sites where the endocannabinoids are produced.This hypothesis has been supported by animal studies in which the FAAH inhibitor URB597 elevated endocannabinoid tone and, unlike the nonselective cannabinoid agonists, did not produce any motor side effects, such as catalepsy.[6] As a result, pharmacological inhibition of intracellular FAAH activity has been the focus of intense drug-discovery efforts and diverse FAAH inhibitors have been developed.[5] The first-generation of FAAH inhibitors were substrate-derived in their structure and designed with reactive groups able to form a covalent bond with the catalytic residue Ser241.[3,4] Although these first-generation inhibitors were capable of blocking FAAH in in vitro pharmacological assays, their lack of selectivity made them poor candidates for preclinical applications.[4] Subsequently, FAAH inhibitors were developed with significantly improved selectivity, including carbamates (ORG-231295), α-ketoheterocycles (OL-135) carbamoyl tetrazoles (LY-2183240), benzothiazole derivatives and piperidine/piperazine ureas [4] (PF-3845, PF-04457845, Figure 1).In our present work, three-dimensional quantitative structure-activity relationships (3D-QSAR) studies based on comparative molecular similarity indices analysis (CoMSIA) studies were carried out on a set of different 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 is revealed by aligning structurally similar analogues using pharmacophoric features as structural superimposition guides.[11] This will allow us to derive predictive 3D-QSAR models for the design of new potent inhibitors and prediction of the activities of new derivatives for this class of compounds.

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.This imbalance in the field contribution of these models made us discard them.The final model selected 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).In order to test the predictive quality of this model, extensive additional validation was carried out.a q 2 = the square of the LOO cross-validation (CV) coefficient; N = the optimum number of components; SEP = standard error of prediction; SEE is the standard error of estimation of non CV analysis; r 2 ncv is the square of the non CV coefficient; F is the F-test value; S, E, H, D and A are the steric, electrostatic, hydrophobic, hydrogenbond donor, and hydrogen-bond acceptor contributions respectively.
Table 2 presents the summary of the external validation of the CoMSIA model.The model has a high value for r 2 test (0.966), which is indicative of an adequate external predictive capacity.However, according to Golbraikh and Tropsha, high values of q 2 and r 2 test (conditions 1 and 2) are necessary but not sufficient conditions for the validation of a model.For a QSAR model to have a reliable predictive capability, the line of experimental versus predicted activity should be as close as possible to the line y = x.This is observed in the fulfillment of conditions [3a or 3b], [4a or 4b] and [5a or 5b] listed in table 2. Finally, condition 6 known as r 2 m metrics, is a quantitative measure to determine the proximity between the observed and the predicted activity for the test set.The CoMSIA model reported here fulfilled all the conditions for internal and external validation.>0.5 0.723 q 2 and r 2 test are the same parameters as listed in table 1; r0 2 and k are the correlation coefficient between the actual and predicted activities for test set and the respective slope of regression; and r0' 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 6.
Furthermore, the Y-randomization test [12] was applied to assess the robustness of the model (see Table S1 of the Supplementary Material for randomizations).The dependent variable (biological activity) was randomly shuffled and a new QSAR model was developed using the original independent variable matrix.If after multiple randomizations the new values of q2 and r2ncv are negative or below the limit of acceptability (q 2 <0.5, r 2 ncv<0.6) it is corroborated that the results obtained in the formulation of the final models are not by chance.In our case, the new QSAR models (after several repetitions) have low q 2 and r 2 ncv values (Table 3).4. All the compounds showed low residual values and deviations in the predicted activity over a logarithmic unit were not observed.Figure 2 shows the graphs of experimental versus predicted activity and it can be seen that the data distribution is close to the y=x line.The model shows a good balance in terms of predictive capacity.The model presents 42 compounds with negative residual values and 48 with positive deviations (Figure 2B).The residual range was -0.82 to 0.89.As shown in Figure 2C 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
The applicability domain (AD) is a theoretical region in chemical space encompassing both the model descriptors and modeled response, which allows one to estimate the uncertainty in the prediction of a compound based on how similar it is to the training compounds employed in the model development.In this work, we used the method developed by Roy et al. [13] for the determination of AD.This method is based on the basic theory of the standardization approach.
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-87.These compounds are the only ones bearing an imidazopyridine or imidazopyrimidine ring and the only difference between them is the position through which the heterocycle is connected to the urea moiety.For this reason, compounds with these heterocyclic systems connected to urea 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 test = 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-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 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, pIC50 = 10.602; on the left) and the least active compound (6, pIC50 = 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 same steric requirement.This can be seen in the proposed molecules (Table 5) where the analog 2x shows a considerable increase in 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 substituentes 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 characteristics 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 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 (Table 1).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.Inside the pyridazine ring the blue polyhedron 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.Proposed analogues 1x, 9x and 10x (Table 5) that follow this substitution pattern display high 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, which 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 can explain the lower activity of this analog.

Hydrophobic Contour Map
The hydrophobic contour maps (Figure 3E and 3F) 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.A yellow polyhedron over the pyrimidine ring shows that a hydrophobic linker region is important and additional yellow polyhedrons surrounding the disubstituted benzene ring suggest hydrophobic substituents in the meta and ortho positions to 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 and 3H) cyan polyhedrons surrounding the pyridazine ring suggest incorporating hydrogen bond donor groups to favor activity.For this reason, the proposed molecule 3x (Table 5) was design 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 and 3J), position 5 of the pyrimidine ring and over the urea carbonyl, which indicates that incorporating hydrogen bond acceptor groups in these positions is unfavorable for activity.Therefore, using different linker groups without hydrogen bond acceptor groups could be advisable to in order to design new 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 bestpredicted inhibitory activity.All proposed molecules have better activity than the most active compound in the series (66, pIC50 = 10.602).In general, the presence of 6-member rings or fused systems on the left side did not greatly increase activity (1x pIC50 = 10.889;4x pIC50 = 10.990).However, the insertion of a pyrazole ring generated derivatives with a significant increase in the pIC50 value (best compounds: 9x, pIC50 = 12.416; 10x, pIC50 = 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.

Selection of Conformers and Molecular Alignment
CoMSIA studies were performed with Sybyl X-1.2 software [14] installed in a Windows 10 environment on a PC with an Intel core i7 CPU.In order to acquire the best conformers for each molecule, every compound was drawing in ChemDraw and then were subjected to a preliminary geometry optimization using MM2 molecular mechanics as is implemented in ChemBio3D software.Following, the structures were further minimized by Tripos force field implemented in Sybyl.MMFF94 charges were assigned to each atom.The minimized structures were superimposed by the atom fit method choosing the piperazinylurea nucleus as the common scaffold for alignment (Figure 4).

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 [15]) 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 attenuation factor α.

Data Set Selection and inhibitory activity
CoMSIA studies were performed on a set of 90 piperazinyl ureas derivatives reported in literature [16][17][18][19][20][21][22][23] (Table 6).The derivatives displayed potent fatty acid amide hydrolase (FAAH) inhibitors activity.The IC50 values were converted to pIC50 (-logIC50).The compounds were randomly 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 pIC50 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 was used to construct a linear correlation between the CoMFA and CoMSIA descriptors (independent variables) and the activity values (dependent variables) [24].To select the best model, the cross-validation analysis was performed by using the LOO method (and SAMPLS), which generates the square of the cross-validation coefficient (q 2 ) and the optimum number of components (N).The non-cross-validation was performed with a column filter value of 2.0 in order to speed up the analysis and reduce the noise.The q 2 , which is a measure of the internal quality of the models was obtained according to the following equation (1): Where   ,  ̅ , and   are observed, mean, and predicted activity in the training set, respectively.The models were subjected to external validation criteria according to the proposed test by Golbraikh and Tropsha [25,26], which considers a QSAR model predictive, if the following conditions are satisfied:  2 > 0.5

External validation
(2) It has been demonstrated [25] that all of the above criteria are indeed necessary to adequately assess the predictive ability of a QSAR model.Furthermore, the external predictive power of the developed 3D-QSAR models using the test set was examined by considering   2 metrics as shown below [27]: Where  2 and  0 2 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   2 should be greater than 0.5.

Applicability domain calculation
The AD was evaluated based on the simple standardization method reported by Roy et al. [13].First, each descriptor "i" for each compound "k" is standardized (Sik).Every compound must have a maximum value [Si]max(k) ≤ 3.In the case that [Si]max(k) > 3 and its minimum value [si]min(k) < 3, then the Snew(k) parameter must be calculated and has to fulfill the condition: Snew(k) =  ̅  + 1.28 *    , where  ̅  is the mean of Sik values for compound k and    is the standard deviation for such values.The software is available free of charge on the authors' website: http://dtclab.webs.com/software-toolsand http://teqip.jdvu.ac.in/QSAR_Tools/.

Conclusions
In this contribution, a 3D-QSAR CoMSIA study was carried out on an extensive database of 90 irreversible inhibitors of the enzyme FAAH of 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.Finally, the information derived from the contour maps was used to design new compounds that showed promising predicted activities (pIC50 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 promising new FAAH inhibitors.

Figure 2 .
Figure 2. A. Plots of experimental versus predicted pIC50 values for the training and test set molecules.B. Residual plots between predicted and experimental values.C. CoMSIA predictions for every molecule in the complete set.

Figure 3 .
Figure 3. CoMSIA steric (A, B), electrostatic (C, D), hydrophobic (E, F), donor (G, H) and acceptor (I, J) contour maps around compounds 66 (left) and 6 (right), the most active and less active of the series respectively.Sterically favored areas are in green and disfavored areas are in yellow.Electropositive favoured areas are in blue and electronegative favoured areas are in red.Hydrophobic favoured areas are in yellow and unfavourable areas in grey.Donor and acceptor favoured areas are in cyan and magenta respectively, and donor and acceptor unfavourable areas are in purple and red, respectively.

Figure 4 .
Figure 4.The superimposed structures of all compounds used in the CoMSIA model.

Figure 5 .
Figure 5. Histogram of frequency distribution data.

Table 1 .
Statistical parameters and Field combinations for CoMSIA.

Table 2 .
Summary of external validation parameters for CoMSIA.

2 ncv Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 8 November 2018 doi:10.20944/preprints201811.0221.v1
The values of experimental activity, predicted activity, and residual values for the best CoMSIA model is shown in Table

Table 4 .
Experimental and predicted pIC50 and residual values for analyzed compounds according CoMSIA.
t test set compound Preprints (www.

Table 5 .
The proposed structures of new molecules and their predicted pIC50 values using the CoMSIA model.

Table 6 .
Chemical structure and pIC50 values of the studied molecules.a