Combined Pharmacophore Modeling, 3D-QSAR, Homology Modeling and Docking Studies on CYP11B1 Inhibitors

The mitochondrial cytochrome P450 enzymes inhibitor steroid 11β-hydroxylase (CYP11B1) can decrease the production of cortisol. Therefore, these inhibitors have an effect in the treatment of Cushing’s syndrome. A pharmacophore model generated by Genetic Algorithm with Linear Assignment for Hypermolecular Alignment of Datasets (GALAHAD) was used to align the compounds and perform comparative molecular field analysis (CoMFA) with Q2 = 0.658, R2 = 0.959. The pharmacophore model contained six hydrophobic regions and one acceptor atom, and electropositive and bulky substituents would be tolerated at the A and B sites, respectively. A three-dimensional quantitative structure-activity relationship (3D-QSAR) study based on the alignment with the atom root mean square (RMS) was applied using comparative molecular field analysis (CoMFA) with Q2 = 0.666, R2 = 0.978, and comparative molecular similarity indices analysis (CoMSIA) with Q2 = 0.721, R2 = 0.972. These results proved that all the models have good predictability of the bioactivities of inhibitors. Furthermore, the QSAR models indicated that a hydrogen bond acceptor substituent would be disfavored at the A and B groups, while hydrophobic groups would be favored at the B site. The three-dimensional (3D) model of the CYP11B1 was generated based on the crystal structure of the CYP11B2 (PDB code 4DVQ). In order to probe the ligand-binding modes, Surflex-dock was employed to dock CYP11B1 inhibitory compounds into the active site of the receptor. The docking result showed that the imidazolidine ring of CYP11B1 inhibitors form H bonds with the amino group of residue Arg155 and Arg519, which suggested that an electronegative substituent at these positions could enhance the activities of compounds. All the models generated by GALAHAD QSAR and Docking methods provide guidance about how to design novel and potential drugs for Cushing’s syndrome treatment.


Introduction
Cortisol is a principal glucocorticoid [1] that is not only used in the treatment of inflammation, allergy, collagen diseases, asthma, adrenocortical deficiency, shock, and some neoplastic conditions, but also exhibits many physiological functions in the regulation of metabolism of life substances, blood pressure and cardiovascular function [2]. Biosynthesis of cortisol take place in the adrenal cortex whose final step (conversion from 11-deoxycortisol) is catalyzed by the mitochondrial cytochrome P450 enzyme steroid 11β-hydroxylase (CYP11B1) [1,3]. (Please reorder references numbers, you jumped 2 and 3 before 4) The secretion of cortisol is precisely controlled by adrenocorticotropic hormone (ACTH) [2,4,5] within the negative feedback cycle of hypothalamic-pituitary-adrenal axis [1]. However, pathological changes in adrenals and the upstream regulating switches can cause an overproduction of cortisol, which is known as Cushing's syndrome. Cushing's syndrome patients mainly show a "moon face", sanguine temperament appearance, obesity, acne, purple lines, high blood pressure, secondary diabetes and osteoporosis, etc. It's a hormonal disorder caused by prolonged exposure to high levels of circulating glucocorticoids such as cortisol [6,7].
Normally, the surgical removal of adrenal or pituitary tumors is used for the treatment of hypercortisolism [8]. However, as mentioned above, CYP11B1 can promote the synthesis of cortisol. Therefore, inhibition of CYP11B1 as the pharmacological approach to block cortisol biosynthesis represents a treatment for Cushing's syndrome [9]. Inhibitors of cortisol biosynthesis, such as ketoconazole, etomidate, and metyrapone have been used in the clinic [4], however, all of them show severe side effects due to the fact that they are unselective. Metyrapone is the only drug reported to be a relatively selective CYP11B1 inhibitor.
In recent years, a series of mitochondrial cytochrome P450 (CYP) superfamily receptors [10], such as CYP1A2, CYP17, CYP19, CYP11B2 [7,11,12], were used to analysis the combination of ligand compounds in different molecular docking studies [13]. Few studies published so far have used the pharmacophore modeling or 3D-QSAR approaches for modeling ligand interactions with the CYP11B1 receptor [14]. In addition, a three-dimensional model of CYP11B1 has not been published yet. Therefore, in the present study we aimed to build and validate homology models of CYP11B1 [15] and then run a docking procedure which indicates the active groups and atoms of inhibitory compounds. In order to analyse the molecular shape of CYP11B1 inhibitors [14], a series of compounds with good inhibitory activities synthesized by Hartmann [15] were collected to establish 3D-QSAR models using CoMFA and CoMSIA. The combination of the pharmacophore model [16] GALAHAD and CoMFA methods also helped establish the structure-activity relationship (SAR) of CYP11B1 inhibitors [17].

GALAHAD Modeling Results
Once GALAHAD modeling based on the training set compounds was completed (Figure 1a), a total of 20 pharmacophore models were generated using common feature hypothesis generation approach. The generated pharmacophore model contained six hydrophobic regions and one acceptor atom They are featured as the cyan balls and the green ball, respectively [18,19]. According to the principals, if all energy parameters had the same level, the Pareto rank would be taken into consideration. Thus, model 1 with a Pareto rank of 0 was selected as the best template to do the CoMFA analysis. The result of this model is Q 2 = 0.658, R 2 = 0.959, F = 82.102, SEE = 0.154 and had two components. The model obtained from combined method had an acceptable balance between energy, pharmacophoric coherence and pharmacosteric overlap statistically. Therefore model 1 was used as the template in aligning the full dataset and did partial least square (PLS) analysis [20,21]. The contour plots between observed and predicted activities of all compounds were shown in Table 1. Nearly, all of compounds were located on the trend line (Figure 1c), indicating that the proposed model was able to successfully predict compounds in test set. In the CoMFA study, the contour maps (Figure 1d) of the pharmacophore model indicated that the blue and yellow contours located around site A would be electropositive groups. The yellow contours located at the B site around the hydrophobe group indicated that a bulky substituent would not be tolerated. The green contours around acceptor atoms and hydrophobe groups indicated that a bulky substituent would be tolerated.

3D-QSAR Modeling Result
The superimposition of all molecules aligned with a common substructure is shown in Figure 2. The aligned molecules were used to generate CoMFA and CoMSIA models which were developed using a combination of different fields, and the statistically significant models were reported with statistical parameters shown in Table 2.  The CoMFA model using both steric and electrostatic fields gave Q 2 of 0.666, R 2 of 0.978, F of 270.441 and SEE of 0.159 values with six components. In the CoMSIA study, the first five models using a single field indicated that the steric field is the most important one. The combination of the steric, electrostatic, hydrophobic, and electrostatic fields led to the S + E + H + A model (Q 2 = 0.721, N = 6), providing the best overall model. The best model led to the highest R 2 , F and the lowest SEE (R 2 = 0.972, F = 209.908 and SEE = 0.180). For these reasons, we considered S+E+H+A to be the best possible combination.

Predictive Power of 3D-QSAR Analyses
CoMFA and CoMSIA analysis developed QSAR models using the training set of CYP11B1 inhibitors. The predicted and observed activities of these compounds were obtained by using the best model that was given in Table 1. The contour plots between observed and predicted bioactivities of training and test set were shown in Figure 3. Most of compounds were located near the trend line, implying the proposed model is able to successfully predict the activities of compounds, which indicated that these 3D-QSAR models are reliable and powerful in predicting pIC50 values. The 3D contour maps were generated to represent the 3D-QSAR results produced by the CoMFA and CoMSIA methods. The different field contributions of COMFA and COMSIA models were illustrated with compound 26 (Figure 4). The results indicated that most of the contours were located at the right side of the compound structures (A and B groups). In the CoMFA model, the contributions of the steric and electrostatic fields to activity were 40.9% and 59.1%. The yellow and blue contours located at the A site indicated that a bulky substituent would not be tolerated and electropositive groups would be favorable [22]. However, bulky substituents and electropositive groups would not be favorable to the B group (Figure 4a) [23]. In the CoMSIA model, the contributions of the steric, electrostatic, hydrophobic and acceptor fields were 13.2%, 40.7%, 22.0% and 24.2%, respectively. The percentages of different fields indicated that the steric and electrostatic fields almost make the same contribution. The combination of these four fields provided the most predictive model for CYP11B1 inhibitors, and the electrostatic field was the most significant for bioactivity prediction. The contour maps of the steric and electrostatic fields of CoMSIA (Figure 4b) were generally in accordance with the field distribution of the CoMFA maps (Figure 4a). The hydrophobic and hydrogen bond acceptor field contour maps of CoMSIA implied that hydrogen bond acceptor substituents were disfavored at the A and B sites, while, hydrophobic groups were favored at the B site (Figure 4c) [24]. The small and electropositive substituent such as an aliphatic amine group would be tolerated at the A site. Small, electronegative and hydrophobic substituent such as trifluoromethanesulfonamide group would be tolerated at B site.

Homology Modeling Result
The length of the CYP11B1 sequence was 574 aa and the most suitable template was the A chain of the CYP11B2 protein (PDB code: 4DVQ). The BLASTP alignment between the CYP11B2 template and CYP11B1 sequences is shown in Figure 5a, revealing 81% identity and 95% consensus similarity. Then, after the structural-based alignment of CYP11B1 from 4DVQ-A, the initial 3D structure of CYP11B1 was obtained from a homology modeling procedure shown in Figure 5d. Evaluation of the homology models consists of Profiles_3D scores and Ramachandran plot analysis [25] (Figure 5c). Profiles_3D scores shown in the Ramachandran Plot of the CYP11B1 model showed a good distribution of 574 amino acid residues of CYP11B1. About 97.0% of the residues featured with green plots were most favored and an additionally allowed region, and only 3.0% residues featured with red plots were an erroneously allowed region.

Docking Analysis
Sixty two (62) CYP11B1 inhibitors were evaluated by docking scores after Surflex-dock. According to the rule, compounds with scores above 7 were considered as active in the docking study. In this paper compounds 33 and 37 scored 7.8544 and 7.1932 (Figure 6a,b) and showed a similar docking mode in the active site of the receptor. The binding cavity of CYP11B1 is formed by residues Arg155, Phe175, Arg519, Cys521, Leu522. Hydrogen atoms of the amino group of residue Arg155 and Arg519 are found near N3 of the imidazolidine ring, thus suggesting that electronegative substituents at these positions could enhance the activities of compounds binding to CYP11B1 [26,27]. These results matched well with the electrostatic contour maps of the CoMFA and CoMSIA models suggesting that an electronegative substituent such as a trifluoromethanesulfonamide would be positive at the B site as mentioned above (Figures 1d and 4a,b).

GALAHAD
Pharmacophore modeling by GALAHAD [11,29] alignment and CoMFA analysis serve as useful tools to produce pharmacophore models of CYP11B1 inhibitors and predict the inhibitory properties of compounds [30,31]. Ten structurally representative compounds with high activities marked with "#" were selected as training set (Table 3) to generate pharmacophore models using the GALAHAD program. The default parameters were used for the GALAHAD runs with the population size turned to 75 and max generations turned to 30. The generated models were evaluated by the test set of CoMFA modeling based on the pharmacophore alignment [30,32]. Therefore the most reasonable pharmacophore model was selected to predict the bioactivities of components and analyse the structures of CYP11B1 inhibitors [33].

3D-QSAR Modeling
The structures of CYP11B1 inhibitors used for the 3D-QSAR study were randomly divided into a training set (42 molecules) and a test set (20 molecules) [34] (Table 1). All structures were energy minimized using the Powell gradient algorithm under the Tripos force field with Gasteiger-Marsili atomic partial charges. The most potent CYP11B1 inhibitor (compound 26) was selected as the alignment template (Figure 7), which was also used as template to build all the other inhibitors with the atom root mean square (RMS) approach. 3D-QSAR models were constructed by using CoMFA and CoMSIA methods based on the molecular alignment. The default values of the parameters of the CoMFA and CoMSIA methods were used. The CoMFA method was performed using steric and electrostatic fields with standard 30 kcal/mol cutoffs. In the CoMSIA study, besides steric and electrostatic fields, three other different fields were calculated: hydrophobic, hydrogen bond donor, and hydrogen bond acceptor [34][35][36]. A series of models were constructed with an increasing number of partial least squares (PLS) analysis factors. The numbers of components in the PLS models were optimized by using the cross-validated correlation coefficient (Q 2 ), non-cross-validated correlation coefficient (R 2 ), standard error estimate (SEE) and F-statistic values (F), etc., which were obtained from the leave-one-out (LOO) cross-validation procedures [37,38]. According these parameters, the best model was chosen to predict bioactivities of compounds. In this work the best 3D-QSAR model was graphically represented by field contour maps, and the coefficients were generated using the StDev⁄Coeff field type.

Homology Modeling and Docking Analysis
The amino acid sequence of CYP11B1 (GenBank: AAH96285.1) was obtained from the National Center for Biotechnology Information Database (NCBI) [39]. Identification of candidate templates was performed by sequence similarity search using BLAST search (NCBI Sever) protocol with default values and each target was searched and downloaded from the NCBI database. The CYP11B1 sequence was aligned to the templates and homology models of CYP11B1 built with the default parameters. Subsequently, the models were analyzed based on Profiles_3D scores and Ramachandran plots which indicate the percentage of amino acids located in the disallowed regions. The interaction of small molecule ligands with a protein, which implied atomic-detail accuracy including position and conformation, was obtained from the docking procedure. Docking and scoring were performed using the Sybyl-X2.1 software (Tripos Inc., St. Louis, MO, USA, 2014) Surflex-Dock method. The docking receptor CYP11B1 was constructed as described in Homology modeling, while the ligands set consisted of all 62 inhibitors. Polar hydrogen atoms were added to both protein and ligand structures. The other parameters of Surflex-Dock methods were used as default values. The best docking pose was selected according to the total score.

Conclusions
CYP11B1 plays a crucial role in the biosynthesis of cortisol which can cause a series of diseases known as Cushing's syndrome. Therefore, CYP11B1 inhibitors that can be regarded as a pharmacological approach to block cortisol biosynthesis have become another treatment for Cushing's syndrome. GALAHAD is a novel pharmacophore screening module, which generates models to analyse the pharmacophore features of ligands. A CoMFA study based on the pharmacophore (GALAHAD) alignment has been developed to derive the structure-activity relationships, which also indicate the interaction between CoMFA fields and pharmacophore features of ligands.
In this study, CoMFA and CoMSIA models were built using the alignment based on the atom root mean square (RMS), which explored the structure-activity relationship of the CYP11B1 inhibitors. These models with excellent consistency manifested good predictive ability for the test compounds. The contour maps also identified the important features contributing to interactions between the different fiedls and the active site of CYP11B1 inhibitors.
Homology modeling method was used to construct a 3D model of the CYP11B1 protein. Then, the Surflex-Dock analysis was used to evaluate the binding activities between the protein and ligand compounds. The binding site indicated the interaction and combination between ligand groups and amino acid residues, which acted as a guide in the prediction and design of CYP11B1 inhibitors.
The combined pharmacophore modeling and 3D-QSAR modeling methods indicated that small and electropositive substituents would be tolerated at the A site. Meanwhile, small, electronegative and hydrophobic substituents would be better at the B site. The docking results indicated that electronegative substituents at the B position could enhance the activities of compounds binding to CYP11B1.