The Discovery of Potentially Selective Human Neuronal Nitric Oxide Synthase (nNOS) Inhibitors: A Combination of Pharmacophore Modelling, CoMFA, Virtual Screening and Molecular Docking Studies

Neuronal nitric oxide synthase (nNOS) plays an important role in neurotransmission and smooth muscle relaxation. Selective inhibition of nNOS over its other isozymes is highly desirable for the treatment of neurodegenerative diseases to avoid undesirable effects. In this study, we present a workflow for the identification and prioritization of compounds as potentially selective human nNOS inhibitors. Three-dimensional pharmacophore models were constructed based on a set of known nNOS inhibitors. The pharmacophore models were evaluated by Pareto surface and CoMFA (Comparative Molecular Field Analysis) analyses. The best pharmacophore model, which included 7 pharmacophore features, was used as a search query in the SPECS database (SPECS®, Delft, The Netherlands). The hit compounds were further filtered by scoring and docking. Ten hits were identified as potential selective nNOS inhibitors.


Introduction
Nitric oxide (NO) is one of the most studied biological signaling molecules and is produced by catalysis from nitric oxide synthase (NOS), which converts L-arginine to L-citrulline, and produces this tiny, short-lived molecule. To date, there are three distinct isoforms of NOS: neuronal NOS (nNOS), endothelial NOS (eNOS) and inducible NOS (iNOS). nNOS and eNOS are constitutively expressed and depend on increases in external calcium and binding of a calcium/calmodulin complex for activation. nNOS and eNOS play an important role in neurotransmission and smooth muscle relaxation, respectively, and iNOS is expressed during bacterial infection, tumor cell cytolysis and inflammation [1][2][3].
As an inorganic reactive free radical gas, NO is believed to be involved in a number of physiological processes such as inflammation, neurotransmission, blood pressure regulation, platelet aggregation and pain [4][5][6]. However, overproduction of NO has been implicated in numerous disease states [7]. In particular, excess NO in the central nervous system from nNOS activity can lead to many neurological disease states including neurodegeneration during Alzheimer's and Parkinson's diseases [8], altered spinal transmission of neuropathic pain [9,10], and progression of migraine and chronic tension-type headaches [11]. Consequently, an inhibitor of nNOS has the potential to be therapeutic in these diseases; however, the functions of eNOS in blood pressure regulation and iNOS in immune responses must be preserved, and the selective inhibition of nNOS has been the challenge of many researchers in the past decade [12][13][14][15].
Recently, several categories of selective inhibitors of nNOS have been designed and developed for the treatment of central nervous system (CNS) disorders [3,[15][16][17][18][19][20][21][22]. Some showed significant efficacy in the rat Chung model of neuropathic pain [22] and in a rodent model of dural inflammation relevant to migraine pain [22]. X-ray structures of nNOS co-crystallized with various ligands [23][24][25] provided insights into the essential structural elements and motifs central to its catalytic mechanism and mode of binding. These findings provide useful information about the interaction between the ligands and the residues near the binding site and can be utilized to design even more selective and potent drug-like NOS inhibitors.
Virtual screening based on a pharmacophore model as a 3D search query has been successfully employed as an efficient alternative to high throughput screening approaches for the development of new compounds with the desired biological properties [26]. Pharmacophore modeling can be used to analyze the common functional groups responsible for specific drug receptor interactions or as a prelude to three dimensional quantitative structure activity relationship (3D-QSAR) analyses that are aligned accordingly with a set of known active compounds in 3D space. 3D-QSAR has been successfully applied in drug discovery and design. As a popular QSAR method, Comparative Molecular Field Analysis (CoMFA) [27] studies incorporate 3D information of the ligands by searching for the sites on molecules that are capable of being modified into more specific ligands. As a useful methodology for studying interaction mechanisms, receptor based molecular docking analysis can be used as a complementary tool to prioritize the hits from the pharmacophore-based virtual screening [28].
In the present study, a 3D pharmacophore model for nNOS inhibitors was assembled and the generated model was used as a search query in the SPECS database containing 197,000 compounds. The virtual screening approach, in combination with pharmacophore modeling and molecular docking can be used to identify and design novel nNOS inhibitors with high selectivity. These molecules may be potential lead compounds for future drug development.

Pharmacophore Results
Twenty pharmacophore models were generated using SYBYL X 1.3 (Tripos Associates Inc., St. Louis, MO, USA). Table 1 lists the parameters of each model. Specificity is a logarithmic indicator of the expected discrimination for each query and is based on the number of features it contains, their allotment across partial match constraints, and the degree to which the features are separated in space. Strong models should have a high Specificity value. Generally, the Specificity value should be at least 5 in a pharmacophore model used as the query for a UNITY flex search [29]. For this study, MODEL 012, MODEL 019, and MODEL 003 had the high Specificity values of 5.138, 5.128 and 4.8580, respectively. These models yielded reasonable pharmacophore models. The N_HITS column shows the actual number of ligands hit by the model query, with the majority of the models matching at least 5 ligands. The value in the FEATS column indicates the total number of features possessed by each model. All of the models had six or more features except for MODEL 011. The retained models had a PARETO rank value of 0, indicating that a single model is not superior to any other. The HBOND term is a measure of the overall pharmacophoric similarity among the ligand conformers. The STERIC term is a measure of the overall steric similarity among the ligand conformers; this term is basically the same as the HBOND term. The ENERGY term indicates the total energy (using the Tripos force field) of all molecules in the training set. The most significant pharmacophore hypothesis was characterized by the conflicting demands of maximizing pharmacophore consensus, maximizing steric consensus, and minimizing conformer potential energy [30]. We constructed a 3D plot to visualize the Pareto surface and select the best pharmacophore model ( Figure 1). Considering only the ENERGY and STERICS criteria, the best model is shown in the upper left-hand corner of the graph in Figure 1b, where the ENERGY is low and the STERICS score is high. In terms of ENERGY and HBOND criteria, the best model is shown in the lower part of the graph in Figure 1c, where the ENERGY is low and the HBOND score is high. Finally, in terms of ENERGY and HBOND, the best model is shown in the upper part of the graph, where both scores are high ( Figure 1d). Among the considered models, MODEL_012 (represented with a red cross in Figure 1) has the optimal position because it fulfills the three criteria and it has the highest Specificity value [31].
The best GALAHAD MODEL 012 is displayed in Figure 2. All of the aligned conformers represent low-energy conformations of the molecules, and the final alignment shows a satisfactory superimposition of the pharmacophoric points. Cyan, magenta, green and red spheres indicate hydrophobes, donor atoms, acceptor atoms and positive nitrogens, respectively. Model 012 includes 7 pharmacophore features: three hydrophobes (HY_1, HY_2 and HY_3), one donor atom (DA_4), one acceptor atom (AA_5) and two positive nitrogens (NP_6 and NP_7). The magenta sphere is covered by a green sphere because the donor atom and the acceptor atom are in the same position in this molecule.

CoMFA (Comparative Molecular Field Analysis) Statistical Results
We used MODEL 012 as a template to align all molecules. The generated steric and electrostatic fields were scaled by the CoMFA-Standard scaling method in SYBYL with the default energy cutoff value. The CoMFA model yielded a good cross-validated correlation coefficient (q 2 ) of 0.513 with an optimized component value of 4, which suggests that the model should be a useful tool for predicting the IC 50 values. A high non-cross-validated correlation coefficient (r 2 ncv ) of 0.933 with a low standard error estimate (SEE) of 0.134 and an f value of 149.950 were obtained. The steric and electrostatic contributions were 45.1% and 54.9%, respectively. The predicted activities for the inhibitors are listed in Table 2 and the correlation between the predicted activities and the experimental activities is depicted in Figure 3. The predictive correlation coefficient (r 2 pred ) was 0.742 for the test set. The statistical results indicate that the CoMFA model is a reliable predictor.    The CoMFA steric and electrostatic contour maps are shown in Figure 4 using compound 41 as a reference structure. In Figure 4a, the blue contour indicates regions in which an increase of positive charge enhances the activity, and the red contour indicates regions in which more negative charges are favorable for activity. The two large blue contours around the red sphere indicate that the substituent in this region should be electron deficient for increased binding affinity with a protein.
Another small blue contour is found around the guanidine isosteric group indicating that a negatively charged substituent in this area is unfavorable. The CoMFA model showed the same result as the pharmacophore hypothesis. In Figure 4b, the steric field is represented by green and yellow contours, in which the green contours indicate regions where a bulky group is favorable and the yellow regions represent regions where a bulky group will decrease activity. In this case, the green contours around the substituent R demonstrated that bulky groups enhance the binding affinity of the nNOS. Most compounds with high activities in this dataset have the same such properties. The CoMFA contour maps and the predicted result further indicated that MODEL 012 can be used as a theoretical screening tool that is able to discriminate between active and inactive molecules [31].

Virtual Screening
The pharmacophore based virtual screening was conducted to find potential nNOS inhibitors. A stepwise virtual screening procedure was applied, wherein the pharmacophore based virtual screening was followed by drug-likeness evaluation, screening of the pharmacophore query, QFIT (The QFIT score is a value between 0 and 100, where 100 is best and represents how close the ligand atoms match the query target coordinates within the range of a spatial constraint tolerance) scoring filtration, and a molecular docking study. The sequential virtual screening flowchart we employed is depicted in Figure 5, in which the reduction in the number of hits for each screening step is shown.

Database Searching
Flexible 3D screening was performed using the UNITY tool to screen the SPECS database [32], which contains approximately 197,000 compounds. The database query was generated based on the pharmacophore MODEL 012. The database was restricted with Lipinski's rule. In general, this rule describes molecules that have drug-like properties. Drug-likeness is a property that is most often used to characterize compound libraries such as combinatorial or screening libraries that are screened to find novel lead chemical compounds [33]. According to this rule, we used simple molecular descriptors, such as molecular weight (≤500), hydrophobicity (MLogP ≤ 4.15) and the number of H-bond donor (≤5) and acceptor atoms (≤10), as the first filter to select the molecules with good absorption or permeation [34]. The remaining 223 compounds were further screened on the basis of QFIT to reduce the dataset, where QFIT is the pharmacophore match between the query and hit [35].

Molecular Docking
To predict the appropriate binding conformation for nNOS inhibitors and the reported hit compounds from virtual screening, Surflex Docking (Tripos Associates Inc., St. Louis, MO, USA) was used to generate an ensemble of docking conformations. The top 62 hit compounds with the highest QFIT score from screening after UNITY filtering were further screened using molecular docking into the binding site of nNOS to select the compounds with the ability to form favorable interactions with the active site. The docked compounds were filtered based on scoring function and interaction with the crucial residues [24] in the binding site. Finally, ten compounds were selected on the basis of the dock score and favorable interactions with the key residues. The results of the hit compounds with their dock score and QFIT values are shown in Table 3. Among the active compounds reported [24], AG_205/36953325 has a similar linker length and two aromatic ring centers on both ends. There is also a hydrogen bond donor in the aromatic ring center and at least one hydrogen bond donor on the linker. The phenolic hydroxyl makes a hydrogen bond with the NOS active site GLU592, which is conserved in all mammalian NOS isoforms [24], and the hydrophobic phenyl ring π-stacks with the heme (HEM801) next to the GLU. The long, flexible linker extending from the phenyl ring allows the 2-phenyl-2,3-dihydro-1H-pyrazole to reach and to π-stack with TYR706 ( Figure 6). The binding mode of this hit compound is similar to that of the reported co-crystallized compound [24] and indicates that the identified hit compounds may have the same mechanism of action as known nNOS inhibitors.

Compounds and Biological Data
Fifty-nine novel nNOS inhibitors were taken from the literature [15][16][17]21,22] with their biological activities in terms of IC 50 values; 49 compounds were used as a training set and the remaining 10 compounds were used as a test set, based on random selection. The compounds in the test set have a range of biological activity values similar to that of the training set. The IC 50 values of the inhibitors were converted into pIC 50 (log (1/IC 50 )) and used as dependent variables in the Pharmacophore generation and CoMFA calculations. The structures of the compounds and their pIC 50 values are given in Table 2. All molecular modeling calculations were conducted using SYBYL X 1.3 (Tripos Associates Inc.). Molecular building was performed with a molecule sketch program in the same software. The molecular geometry of each compound was first minimized using the standard Tripos molecular mechanics force field with 0.01 kcal/(mol A) energy gradient convergence criterion. Partial atomic charges were calculated by the Gasteiger-Hückel method and energy minimizations were performed using the Conjugate Gradient method with 1000 iterations [36,37].

Pharmacophore Generation
Pharmacophore models were generated and analyzed using the GALAHAD module. In this study, ten compounds (13, 16, 19, 21, 24, 30, 40, 41, 43, and 52) were selected to carry out the pharmacophore hypothesis, and the genetic algorithm was used to create conformers for all molecules. The compounds that were selected to generate the pharmacophore hypothesis are highly active. All of the ligands were aligned with a population size value of 60, a maximum generation value of 60 and a value of molecular required hitting of 5. Twenty models were generated with default parameters.

CoMFA Field Calculation Partial Least Square Analysis
The standard CoMFA procedure as implemented by SYBYL X 1.3 (Tripos Associates Inc.) was performed. Each set of aligned molecules was positioned inside a 3D cubic lattice with a grid spacing of 2.0 A (default distance) in all Cartesian directions and was generated to enclose the molecule aggregate. The fields generated were automatically scaled by the CoMFA standard in SYBYL. The partial least squares (PLS) methodology was used to derive a linear relationship for the CoMFA, and cross-validation was performed using the leave-one-out (LOO) method to choose the optimum number of components (ONC) and assess the statistical significance of each model. In PLS, the independent variables were the CoMFA descriptors, and the pIC 50 values were used as dependent variables. The ONC was the number of components that led to the highest cross-validated correlated correlation coefficient q 2 (or r 2 cv ). Non-cross-validation was performed to calculate conventional r 2 ncv using the same number of components [38][39][40].

Virtual Screening
The selected pharmacophore model was validated and converted into a UNITY query for pharmacophore guided virtual screening studies. The query was screened against the SPECS database.
The flex search option was implemented to perform virtual screening. Primary filters such as Lipinski's rule of five were applied to reduce the dataset. Further screening of the hits was carried out using the Surflex Dock in SYBYL.
The docking study was performed to validate the hits obtained from the virtual screening. The crystal structure of nNOS was retrieved from the RCSB Protein Data Bank (PDB code: 4EUX) [24]. The nNOS structure was utilized in subsequent docking experiments without energy minimization. Protein structures were prepared using the biopolymer module of SYBYL. Hydrogen atoms were added to the structure, atom types and charges were assigned using AMBER7 FF99 force field, and side chain amides were modified. The ligand method was used as the mode of construction for the protomol, threshold and bloat using default values to determine the extent of the protomol.

Conclusions
nNOS is a therapeutic target for central nervous system diseases that has attracted interest from pharmaceutical companies and researchers. Selective inhibition of nNOS activity represents an exciting drug approach for the development of new therapeutic agents to treat neurodegenerative diseases. In this study, we described a rational strategy for identifying novel nNOS inhibitors using a pharmacophore-based virtual screening protocol. The best pharmacophore model (MODEL 012) was established and showed good statistical parameters in the validation process. MODEL 012 was further employed as a 3D search query to screen the SPECS compound database. Molecular docking studies were also performed to improve the reliability and accuracy of the virtual screening. Ten hit compounds were identified as potential selective nNOS inhibitors and exhibited good search scoring, high docking scores, similar binding mode to experimentally proven compounds and favorable drug-like properties. The pharmacophore models developed in this work, and the information gained about the interactions between nNOS and the potential selective inhibitors, indicated that the combination of pharmacophore, molecular docking, and virtual screening efforts is a successful approach for identifying effective inhibitory compounds that may have an impact on future experimental studies in selective nNOS inhibition. The identified hit compounds were structurally different from available inhibitors and may serve as potential leads or starting points for structural optimization to identify novel nNOS inhibitors.