Homology Model of Human Retinoic Acid Metabolising Enzyme Cytochrome P450 26A1 (CYP26A1): Active Site Architecture and Ligand Binding

Abstract

Homology models of cytochrome P450 RA1 (CYP26A1) were constructed using three human P450 structures, CYP2C8, CYP2C9 and CYP3A4 as templates for the model building. Using MOE software the lowest energy CYP26A1 model was then assessed for stereochemical quality and side chain environment. Further active site optimisation of the CYP26A1 model built using the CYP3A4 template was performed by molecular dynamics to generate a final CYP26A1 model. The natural substrate, all-trans-retinoic acid (atRA), and inhibitor R115866, were docked into the model allowing further validation of the active site architecture. Using the docking studies structurally and functionally important residues were identified with subsequent characterisation of secondary structure. Multiple hydrophobic interactions, including the side chains of TRP112, PHE299, PHE222, PHE84, PHE374 and PRO371, are important for binding of atRA and R115866. Additional hydrogen bonding interactions were noted as follows: atRA – C = O of the atRA carboxylate group and ARG86; R115866 – benzothiazole nitrogen and the backbone NH of SER115.

Keywords::

• CYP26A1
• homology model
• MOE
• all-trans-retinoic acid (atRA)
• docking studies
• inhibitor interactions
• inhibition
• R115866

Introduction

Retinoic acid (RA) regulates genes involved in cell proliferation, differentiation and apoptosis [1], and exerts activity by binding to transcription-regulatory factors in the cell nucleus known as RAR (retinoic acid receptor) and RXR (retinoid X receptor), each having subtypes α, β and δ [2]. RAR receptors are activated by all-trans- (atRA) and 9-cis- (9-cis-RA) isomers of RA, whereas RXR receptors are only activated by 9-cis-RA. Upon RA binding the activated receptor transcriptionally regulates its target genes by binding to specific response elements in promoter regions of RA target genes (Figure 1) [3,4].

Figure 1 Alcohol dehydrogenases (ADH) and short-chain dehydrogenase/reductase catalyse the oxidation of retinol to retinaldehyde, which is subsequently oxidised by aldehyde dehydrogenases (ALDH) to retinoic acid. Retinoic acid, all-trans (atRA) and 9-cis (9cRA) isomers are further metabolised by cytochrome P450 enzymes (including CYP26A1, B1 and C1) to inactive polar metabolites resulting in retinoid excretion.

Retinoic acid is synthesised from vitamin A (retinol), which is oxidised through retinal by dehydrogenases in the cytoplasm of target cells, in low yields, to all-trans-retinoic acid (Figure 1). atRA is subsequently metabolised by human liver and intestine cytochrome P450s to the inactive 4-hydroxy-RA and thence by dehydrogenases to the partially active 4-keto-RA and inactive polar metabolites (Figure 1) [5].

The P450s responsible for 4-hydroxylation of atRA in the human liver are CYP2C8 as a major contributor as well as 3A7, 3A5, 3A4, 2C9 and 1A1 [6]. However, in living tissues, atRA administration induces another RA-metabolising enzyme, CYP26 [7], which recognises only atRA as its substrate, and the expression of this isozyme can be induced by atRA both in vitro and in vivo [8]. Three members of the CYP26 family have now been identified: CYP26A1 [7] and CYP26B1 [9], which metabolise atRA in the embryo and adult and, more recently [10], CYP26C1 that may have a role in the metabolism/isomerisation of both all-trans and 9-cis isomers of RA.

Retinoids in general have been used for some time in the treatment of psoriasis, cystic acne, cutaneous malignancies due to hyperkeratinisation as well as in the treatment of photo-damaged skin [11,12]. Retinoic acid has been used in a number of clinical situations, especially oncology and dermatology, atRA has also been shown to improve the efficacy of other treatments such as radiation, cisplatin and interferon therapies [13,14].

Although the amino acid sequence of CYP26A1 has been characterised no crystal structure is available. A homology model for human CYP26A1 would provide valuable information regarding active site architecture and allow docking studies to analyse ligand binding interactions, the information obtained from these studies would assist with rationale drug design of retinoid mimetics and CYP26A1 inhibitors. Importantly, crystal structures for three human cytochrome P450 enzymes, CYP2C8 [15], CYP2C9 [16] and CYP3A4 [17] have recently become available allowing the use of a human crystal template for the homology building.

Materials and methods

Computational approaches

All molecular modelling studies were performed on a RM Innovator Pentium IV 2.4 GHz running either Linux Fedora Core 3 or Windows XP using Molecular Operating Environment (MOE) 2004.03 [18] and FlexX module in SYBYL 7.0 [19] molecular modelling software.

Molecular dynamics simulations were performed with GROMACS 3.2 [20,21] and the Gromacs force field in a NVT (canonical) environment. Individual ligand/protein complexes obtained from the docking results were soaked in a triclinic water box and minimised using a steepest descent algorithm to remove unfavorable van der Waals contacts. The system was then equilibrated via a 20 ps MD simulation at 300°K with restrained ligand/protein complex atoms. Finally, a 800 ps simulation was performed at 300°K with a time step of 2 fs and hydrogen atoms constrained with a LINCS algorithm. Visualisation of the dynamics trajectories was performed with the VMD software package, version 1.8.3 [22].

All the minimisations were performed with MOE until RMSD gradient of 0.05 Kcal mol^{− 1}Å^{− 1} with the forcefield specified and the partial charges were automatically calculated.

Homology searching

The protein sequence of human CYP26A1 was obtained from the ExPASy server (O43174) [23]. Homologous proteins with known crystal structures were found by performing a PSI-BLAST search (comparison matrix, BLOSUM62; E-threshold, 10), using the ExPASy server, aligning the query sequence (CYP26A1) against sequences in the Protein Data Bank (PDB) [24]. The alignment parameters and thresholds used for screening candidate homologues were used with their default values and BLOSUM62 comparison matrix.

Multiple sequence and structure alignment

MOE-Align, which implements a modified version of the alignment methodology originally introduced by Needleman and Wunsch [25], was employed. The sequence alignment was performed with an alignment constraint between the haem cysteine residue of the query sequence (CYP26A1) and the corresponding haem cysteine residue of the templates. All the default settings in the MOE-Align panel were used for the sequence alignment.

D Model Building

Homology models were built using MOE-Homology using a single template approach with CHARMM22 forcefield. The three human P450 crystal structures CYP2C8, CYP2C9 and CY3A4 were selected as the templates for constructing CYP26A1 models. Ten intermediate models were generated and the final model was taken as the Cartesian average of all the intermediate models. The haem was positioned using the same coordinates as in the template and the complex model was energy minimised.

Model Validation

Stereochemical quality of the polypeptide backbone and side chains was evaluated using Ramachandran plots obtained from the RAMPAGE server [26]. Amino acid environment was evaluated using Verify3D [27] and Errat plots[28]. Verify 3D assesses environment of the sidechain based on the solvent accessibility of the sidechain and the fraction of the sidechain covered by polar atoms. Errat assesses the distribution of different types of atoms with respect to one another in the protein model. Validation data for the templates CYP2C8, CYP2C9 and CYP3A4 was used as the baseline to assess the respective models.

Docking

Ligands were docked within the active site of the homology model using the FlexX docking programme of SYBYL, performed with the default values. The active site was defined by all the amino acid residues within a 6.5 Å distance from TRP112, VAL116, THR304, VAL370 and GLY373, including the haem in a heteroatom file.

The output of FlexX docking was visualised in MOE and the scoring.svl script [29] was used to identify interaction types between ligand and protein.

Results and discussion

Comparative modelling methods use structural templates that have the highest sequence homology with the target protein. Homologous proteins were identified by scanning the protein sequence of CYP26A1 [30], obtained from the ExPASy server [23], against 3D structures deposited in the Protein Data Bank (PDB) [24] using PSI-BLAST [31]. The search returned amino acid sequences of different P450s isolated from different species, including the three human P450s recently deposited in the PDB. The percent sequence identity, chain length and E-value for the homologous P450 sequences are shown in Table I. CYP3A4 (1TQN), CYP2C8 (1PQ2) and CYP2C9 (1R90) were promising templates with percent sequence identity ranging from 24-26% and good E-values.

Table I.  P450 sequences identified by PSI-BLAST search using CYP26A1 as the query sequence.

Protein	PDB code	PSI-BLAST score^a	Sequence identity^b	% Sequence identity	Chain length	E-value
CYP51	1EA1-A	275	115/447	25	455	1.3e-22
CYP102A1	1JPZ-A	253	113/432	25	473	5.7e-20
CYP3A4	1TQN-A	249	111/437	26	486	1.8e-19
CYP2C9	1R90-A	171	52/203	25	477	3.8e-13
CYP2B4	1PO5-A	164	71/285	24	476	6.9e-16
CYP2C8	1PQ2-A	156	67/277	24	476	1.1e-11
CYP2C5	1DT6-A	153	53/219	24	473	4.5e-12
CYP108	1CPT	153	42/133	31	412	4.0e-10

^a The PSI-BLAST score for an alignment is calculated by summing the scores for each aligned position and the scores for gaps.

^b (Number of identical residues)/(length of sequence fragment identified by PSI-BLAST).

All three sequences have a similar chain length to the query sequence CYP26A1 and are high resolution structures. Structural alignment of the CYP26A1 with the three human P450 templates was performed using ClustalW [32] on the EBI server [33].

Alignment allowed the assignment of substrate recognition sites and P450 binding motifs (ETLR, PERF and haem binding domains) (Figure 2).

Figure 2 CYP26A1 sequence and ClustalW (1.82) alignment with CYP2C9, CYP2C8 and CYP3A4. The colour coding of amino acid type: red, small + hydrophobic (AVFPMILWY); blue, acid (DE); magenta, basic (RHK); green, hydroxyl + amine + basic (STYHCNGQ). Substrate recognition sites (SRS) and conserved P450 domains (ETLR, PERF and Haem binding domain) are indicated. “*” means that the residues are identical, “:” means that conserved substitutions have been observed, “.” means that semi-conserved substitutions are observed.

Models, based on each template using a single template alignment approach, were generated using Molecular Operating Environment (MOE) [18] software as described in the Experimental section. To evaluate the quality of the modelled structures, the lowest energy model generated from each template was subjected to a number of checks. Stereochemical quality was assessed using Ramachandran plots [34], using the Cambridge RAMPAGE server [26], and amino acid environment was assessed using Verify3D [27] and Errat [28], using the UCLA-DOE server [35] (Table II).

Table II.  Validation results for the lowest energy CYP26A1 models (produced from the three different templates) and the 3D structural templates.

	Ramachandran plot^a (%)	Errat (%)	Verify 3D^b (total score)
Model
CYP26A1-CYP3A4 model	81.4	84.5	121
CYP26A1-CYP2C8 model	82.3	84.5	137
CYP26A1-CYP2C9 model	82.4	82.3	126
3D Template (Resolution)
CYP3A4 (2.05 Å)	91.2	93.7	178
CYP2C8 (2.70 Å)	92.3	89.6	183
CYP2C9 (2.00 Å)	93.8	90.8	161

Validation results for the final CYP26A1 models built from the 1jpzB template.

^a Percentage of residues with ϕ, ψ conformation in the ’most favoured’ regions of the Ramachandran plot.

^b The total Verify 3D score summed over all of the residues.

Validation results would suggest that all three models performed equally well in terms of mainchain stereochemistry and amino acid environment, therefore substrate and ligand docking was required for further validation of the active site architecture.

Ligand Docking

atRA was found to dock satisfactorily in the active site of all three models, with atRA orientated in a position that favours 4-hydroxylation at C-4. However, only the CYP26A1 model built using the CYP3A4 template was able to accommodate the inhibitor ligand (S)-R115866 [8] in an orientation that would allow coordination between the nitrogen of the triazole heterocycle and the haem iron transition metal. However, docking studies had indicated that TRP112 was partially obstructing coordination between the triazole nitrogen of the ligand R115866 and the transition metal resulting in a distance of 5.50 Å (Figure 3a). Docking studies do not take into account protein flexibility therefore molecular dynamics were performed on the active site containing (S)-R115866 resulting in an optimised active site architecture. This resulted in the TRP112 residue positioned in a more favourable conformation with respect to ligand binding (Figure 3b) with the triazole nitrogen perpendicular to the haem iron at a distance of 2.6 Å.

Figure 3 CYP26A1 model active site (a) before and (b) after active site optimisation with the (S)-R115866 bound inhibitor.

Furthermore, (S)-R115688 establishes a hydrogen bond between the benzothiazole nitrogen and the NH of SER115, as well as several hydrophobic interactions with the side chains of TRP112, PHE374, PHE84, PHE299, VAL116, PRO371 and other residues (Figure 4).

Figure 4 Interactions between (S)-R115866 and the CYP26A1 active site.

The natural substrate atRA was docked with the optimised active site. Figure 5 shows the putative active site of CYP26A1 containing the bound substrate orientated for 4-hydroxylation with the C4 atom positioned above the haem iron at a distance of 5.3 Å, this distance would accommodate a water molecule between the 4-position of atRA and the haem iron. atRA interacts with amino acid residues at the active site by multiple hydrophobic interactions, including the side chains of TRP112, PHE299, PHE222, PHE84, PHE374 and PRO371. Hydrogen bonding interactions between the carbonyl group of atRA and ARG86 hold the molecule within the hydrophobic tunnel.

Figure 5 atRA docked in CYP26A1 active site.

The secondary structure of the CYP26A1 model was determined using Swiss-PdbViewer 3.7 [36] with identification of α–helices, β-sheets, coils and loops as shown in Table III. The secondary structure is derived from the template structure CYP3A4 therefore further comparison with the template allows identification of specific helices and β-sheets. Studies have shown that the active site region (Figure 2) of the P450s is conserved across the entire family [37] and the middle section of the I-helix (Figure 2) is also well conserved despite low overall sequence homology [38]. However, although similarities do exist across the family, there are flexible regions in each P450 which have significantly different amino acid composition, such as the B'-helix, which may or may not be present, and F-G loop regions.

Table III.  3. Comparison of CYP26A1 model and CYP3A4 template secondary structure.

CYP26A1 Model		CYP3A4 Template
Residues	Secondary structure	Residues	Secondary structure
LEU58-MET60	Helix	LEU32-LEU36	Helix
ARG64-TYR75	A-Helix	PHE57-TYR68	A-Helix
TYR79-LEU83	β-Sheet	LYS70-GLY77	β-Sheet
ARG86-VAL91	β-Sheet	GLN79-THR85	β-Sheet
ALA94-LEU101	B-Helix	PRO87-LEU94	B-Helix
ASP105-LEU107		MET114-LYS115	B'-Helix
SER132-ARG142	C-Helix	GLU124-ARG128	C-Helix
ARG146-SER169	D-Helix	LEU132-PRO135
ARG173-VAL177	D'-Helix	SER139-GLU165	D-Helix
PHE186-GLY195	E-Helix	LYS173-SER186	E-Helix
ASP227-VAL228	F'-Helix	PRO202-LYS208	F-Helix
GLY232-ALA239	G-Helix	PRO218-THR224	F'-Helix
GLU248-GLU262		PHE228-LEU236	G-Helix
GLN288-GLY317	I-Helix	ARG243-LYS257
HIS321-LYS331	J-Helix	PHE271-GLN279	H-Helix
ASP346-ILE349	J'-Helix	ASP292-LEU321	I-Helix
VAL359-LEU366	K-Helix	PRO325-VAL338	J-Helix
GLY373-ALA377	β-Sheet	TYR347-LEU351	J'-Helix
TRP392-SER397	β-Sheet	TYR355-LEU366	K-Helix
CYS399-HIS402	K'-Helix	MET371-CYS377	β-Sheet
PRO422-GLU425	Helix	VAL381-ILE383	β-Sheet
GLU446-ARG461	L-Helix	MET386-ILE388	β-Sheet
ASP464-GLN466	β-Sheet	VAL392-SER398	β-Sheet
THR473-LYS475	β-Sheet	TYR399-LEU401	K'-Helix
TYR481-VAL483	β-Sheet	GLU417-ARG418	Helix
ARG489-THR491	β-Sheet	MET445-LEU460	L-Helix
		ASN462-LYS466	β-Sheet
		LYS476-SER478	β-Sheet
		GLN484-GLU486	β-Sheet
		VAL489-ARG496	β-Sheet

A good overlap between the model and template structure was found in particular for the A-, E- G-, I-, J-, J'- and K-helices. The model differs from the template in that the B'-, F- and H-helices are not identified in the model. As seen in Figure 6, the large I-helix lies above the haem and is connected to the J-helix through a loop. The active site cysteine 442 lies below the haem within a coil and it is also worth noting that all the residues involved in the binding of atRA and (S)-R115688, are located in the substrate recognition sites (SRSs)(Figure 2).

Figure 6 3D Structure of the CYP26A1 model. α-Helices are in red, β-sheets are in yellow, loops are cyan and coils are grey.

Conclusions

The CYP26A1 model has been validated for stereochemical and amino acid environment quality using appropriate programmes, with further validation of active site architecture achieved by docking studies with the natural substrate and an inhibitor to determine the key residues involved in the ligand binding. A drug design program of CYP26A1 inhibitors using this CYP26A1 model is now ongoing based on the results obtained from the work presented here.

Acknowledgements

For MSG we would like to acknowledge the Embassy of the Arab Republic of Egypt for the award of a PhD scholarship, and for SWY we would like to acknowledge the ORS Awards Scheme for a United Kingdom Scholarship.