Binding of Chlorinated Phenylacrylonitriles to the Aryl Hydrocarbon Receptor: Computational Docking and Molecular Dynamics Simulations

Abstract

The development of ligands capable of binding to the aryl hydrocarbon receptor (AhR) and hijacking its signaling pathway is of potential use for the design of novel agents against breast cancer. To guide the synthesis of new compounds and characterize their binding to the AhR, we employed homology modeling, ligand docking, and molecular dynamics simulations. As there is currently no crystallographic information available for the structure of the AhR’s ligand-binding PAS-B domain, a homology model was developed. The location of the binding site was identified by scanning the model for concave areas and comparing them to known ligand-binding sites in proteins related to the AhR, such as the CLOCK:BMAL1 transcriptional activator complex and the hypoxia-inducible factor-2α (HIF-2α). Docking of several chlorinated phenylacrylonitriles was performed with the modeling suite MOE, identifying π-π stacking, hydrophobic, and van der Waals interactions as the driving forces for binding, an observation consistent with the hydrophobic nature of the site. Molecular dynamics simulations with one of the compounds for 100 ns verified the overall stability of a docking-predicted pose and revealed the presence of interacting water molecules in the vicinity of the ligand. Given the buried location of the ligand in the core of the receptor, this observation was somewhat unexpected, but it explained a slight shift of the ligand pose seen during the simulation.

1. Introduction

The aryl hydrocarbon receptor (AhR) is a member of the basic helix-loop-helix/Per-ARNT-SIM (bHLH/PAS) transcription factor family [1, 2, 3, 4]. In its inactive state, the AhR resides in the cytosol of the cell as a complex with a number of other proteins. This complex ensures the stability of the AhR in a high-affinity ligand-binding form and prevents the premature translocation of the receptor. Upon binding of a ligand, it dissociates from these proteins and travels to the cell nucleus, where it binds to DNA xenobiotic response elements (XREs). This in turn induces the expression of several cytochrome P450 enzymes and a sulfotransferase (typically SULT1A1) that contain XREs in their promotor sequence. These enzymes then initiate the oxidative breakdown of the offending compound.

The AhR pathway has a number of roles, including as a modulator of viral immunity and the correct functioning of the female reproductive system. Its most well-known role is a mechanism by which cells defend themselves against the toxic effects of polycyclic and polyhalogenated aromatic hydrocarbons, such as the Seveso toxin dioxin (1) (Figure 1) [5, 6].

Figure 1.

The aryl hydrocarbon receptor ligand 2,2,6,6-tetrachlorodioxin (1).

Hijacking of the pathway is based on the use of compounds capable of activating the pathway and then converting into highly reactive species such as nitrenes once being targeted by the metabolic enzymes. This process ultimately leads to DNA damage and the death of the affected cell (Figure 2) [7].

Figure 2.

The AhR pathway showing ligand binding, nuclear translocation, CYP1 activation, metabolism, and cell death. AF = Aminoflavone (2), Phort = Phortress (3), and ANI-7 = (Z)-2-(3,4-dichlorophenyl)-3-(1H-pyrrol-2-yl)acrylonitrile (4) (see Figure 3 for details).

It has been noted that the AhR detoxification process involves the active transport of a ligand, e.g., 1–4, but not the inhibition of the AhR, which would result in a buildup of toxic materials within the cell. This hijacking of the AhR signaling pathway has been proposed as a novel strategy for designing a new class of drugs against breast cancer [1, 8]. Several compound classes, such as the aromatic acrylonitriles, have shown promise in cell-based assays, displaying remarkable potency and selectivity for breast cancer cells [9, 10]. Two reported AhR ligands, Aminoflavone (2) and Phortress (3) (Figure 3), have progressed to clinical trials, demonstrating the clinical applicability of this approach [11, 12]. Based on this, we have postulated that the AhR is a promising target in the development of breast cancer-specific drugs. In particular, our early studies have demonstrated activity against triple negative breast cancer cell lines [9, 10, 13]. This makes AhR ligands, including the aromatic acrylonitriles, promising candidates for further development into novel agents against breast cancer, that act by a hitherto unexploited mechanism of action.

Figure 3.

The known AhR ligands, Aminoflavone (2) and Phortress (3), that have proceeded to clinical trials for the treatment of cancer and our recently reported lead AhR ligand, ANI-7 (4) [13].

Here, we demonstrate the use of computational tools for the elucidation of the interactions between the AhR and a targeted selection of chlorinated phenylacrylonitriles. The methods employed include homology modeling, molecular docking, and molecular dynamics (MD) simulations to model the structure of the ligand-binding domain of the AhR, identify its ligand-binding site, characterize critical ligand/receptor interactions, and study the time-dependent behavior of a ligand bound to the AhR. The results illustrate the value of computational tools for revealing the potential binding mechanism of these compounds to their target and for guiding the synthesis of novel compounds with improved properties.

2. Homology model

The sequence of the human form of the AhR was downloaded from the NCBI website (access code: NP_001612.1). Since only the ligand-binding PAS-B domain was of interest to our study, the sequence was appropriately truncated before Pro275 and after Lys397. A search in the modeling suite MOE’s structural database for suitable templates returned the structures of 4F3L [14], 3RTY [15], and 2KDK [16] as the best matches. Of these, only 4F3L, a murine transcriptional activator complex, provided complete coverage of the PAS-B domain with a sequence identity of 24.4% and a sequence similarity of 48.0% (Figure 4). Only three indels were noted in the alignment—deletions of positions 361 and 362 in the target sequence and an insertion in position 308. The absence of major gaps in the alignment is favorable for the development of homology models as it reduces the need for loop modeling and grafting, which can be challenging [17]. Model development based on the alignment in Figure 4 was performed using MOE’s default settings.

Figure 4.

Alignment of the target sequence of the human form of the AhR (NP_00161) with the sequence of a murine transcriptional activator complex, 4F3L.

The resulting homology model of the AhR (Figure 5A) was subjected to a number of quality tests, such as an analysis of the Ramachandran diagram (Figure 6) and an inspection of observed bond lengths, bond energies, and torsion angles. No abnormalities that would have questioned the quality of the model were detected.

Figure 5.

(A) Homology model for the AhR colored by secondary structure. (B) Comparison of backbone traces of homology models obtained by using the MOE modeling suite (template 4F3L) and the automated SWISS-MODEL server (template 5SY7). Coloring is according to RMSD between the two structures (green—yellow—red, in order of increasing deviation), showing very good agreement between the two models. The model obtained from the SWISS-MODEL server had a somewhat longer sequence, resulting in the gray loops at the termini that have no counterpart in the MOE model.

Figure 6.

Ramachandran diagram for the homology model for the AhR. Green ( ) symbols represent torsion angles in favored regions, whereas yellow ( ) symbols represent angles in allowed regions. No entries are present in the “forbidden” areas.

An additional check of the model’s reliability was carried out by submitting the PAS-B sequence to the automated server SWISS-MODEL [18]. The returned homology model was superimposed to the model obtained from MOE. Even though the new model was derived using a different template (5SY7, an NPAS3-ARNT complex) [19], a very good agreement between the backbones of the two structures was observed (Figure 5B), which further instilled confidence in the accuracy of the model.

3. Computational ligand docking

Before ligands could be docked into the homology model of the AhR, the exact location of the binding site had to be identified. We subjected the homology model to a binding site search, a feature implemented in MOE that screens the surface of a protein for concave areas capable of binding small molecules. Two areas large enough to accommodate a typical AhR ligand were detected: one on the surface and another one in the core of the receptor. To decide which of these two sites was more realistic, the crystal structures of the ligand/receptor complexes 3F1O [20], 3H7W [21], and 3H82 [21], all of which are proteins related to the AhR, were superimposed onto the homology model. As shown in Figure 7, all three ligands were found in an area equivalent to the binding site located at the center of the protein (Figure 5). To facilitate a convenient designation of the binding site for the subsequent docking runs, the ligand of 3F1O—N-[2-nitro-4-(trifluoromethyl)phenyl]morpholin-4-amine (5)—was copied into the file of the homology model as a point of reference.

Figure 7.

(A) Superposition of protein/ligand complexes related to the AhR onto the homology model of the AhR. Spheres delineate the putative binding site predicted by MOE that coincides with the position of the ligands seen in the crystal structures. (B) A closeup view of the ligand 5, overlaid with the spheres depicting the binding site predicted by MOE. (C) The chemical structure of the ligand of 3F1O—N-[2-nitro-4-(trifluoromethyl)phenyl]morpholin-4-amine (5).

To analyze the utility of this model for further drug development, the structures of a representative ensemble of six dichlorophenylacrylonitriles with known bioactivities (Figure 8) were modeled in MOE. Their conformational energies were minimized by molecular mechanics in conjunction with the MMFF94x force field. Docking was performed with the default settings of MOE, utilizing a flexible ligand and a mostly static receptor structure and defining the binding site by a position equivalent to that of the ligand present in 3F1O. The top-scoring pose for each ligand was considered for further analysis.

Figure 8.

Structures of six dichlorophenylacrylonitriles (4, 6–10) used for docking.

As shown in Figure 9, the docked compounds occupied a narrow and mostly hydrophobic site in the core of the AhR. Almost all ligands in the pool engaged with the AhR in a similar fashion, binding in comparable binding poses and exhibiting similar ligand/receptor interactions. Key hydrophobic contacts were observed between nonpolar regions of the ligands and the side chains (Phe21, Leu34, Phe50, Met66, Leu79, Ala93, Ile105, and Val107). Moreover, the ligand phenyl ring engaged in π-π stacking interactions with the ring of His 17. In addition, the tight fit between the ligands and the site suggested the presence of extensive favorable van der Waals interactions.

Figure 9.

(A) 3D representation of dichlorophenylacrylonitrile 7 docked into the binding site of the AhR, illustrating the central nature of the site. (B) Interaction diagram of the pyrrole ligand ANI-7 (4), showing π-π -interactions, hydrophobic contacts, and shape complementarity as the driving forces for ligand binding.

In an attempt to attain a quantitative measure of ligand-binding affinity, the docking scores of the compounds were graphed against observed potencies (Figure 10). Potencies had been obtained in cell viability assays and the underlying assumption was that ligand binding to the receptor constituted the critical step that would lead to cell toxicity and therefore would correlate with bioactivity. Figure 10 shows a reasonable correlation between the two quantities with a squared correlation coefficient of 0.79. This data is consistent with the proposed binding mode of AhR ligands, which relies predominately on hydrophobic but also on additional π-π interactions.

Figure 10.

Linear correlation between docking score and bioactivity (−log IC₅₀). Docking was performed with the standard settings of MOE.

4. Molecular dynamics simulations

We complemented our docking-based analysis by MD simulations, whose purpose was twofold. First, we wanted to ensure the stability of a docked pose by monitoring its behavior in a time-resolved system. Second, MD simulations can reveal the role of explicit solvent molecules, something that cannot be accounted for by docking. We selected compound 7 as a representative and subjected it to a simulation time of 100 ns, using the CHARMM36m force field for the protein [22] and the CHARMM general force field for the ligand [23]. The parameters for water were taken from the CHARMM-modified TIP3P water model [24, 25, 26] to match those used for the solute. The initial structure of the protein-ligand complex was obtained from docking experiments, and the simulation was performed with the software NAMD [27].

As shown in Figure 11, the differences between the poses before and after 100 ns of simulation time were minor. The overall position of the ligand did not change significantly and the only notable difference related to a slight rotation around the central axis of the molecule which placed the nitrile group in a somewhat different environment.

Figure 11.

MD simulations. (A) Ligand/AhR complex before (gray) and after (blue) 100 ns of simulation time. (B) Density of water molecules in the binding site as highlighted by orange grids. Residues are colored according to their polarity (pink: nonpolar; green: polar uncharged).

Interestingly, analysis of the MD simulation data revealed the presence of several water molecules in close proximity to the ligand. This observation was somewhat unexpected; while polar water molecules have been found in predominately hydrophobic cores of proteins, it is a rare occurrence [28, 29]. Residues exposed to water molecules included Leu315, Thr289, His 291, Gln383, Ser365, and His 337. In some cases, the solvent molecules formed bridged hydrogen bonds between the nitrile group and the ligand. The latter could explain the abovementioned slight twist of the nitrile group into a more favorable position.

5. Conclusions

Using the AhR and substituted phenylacrylonitriles as an example, we demonstrated the usefulness of a number of computational tools for the study of ligand/receptor interactions. Homology modeling gave access to the structure of a protein domain that has not yet been solved by X-ray crystallography. The most probable binding site was identified, allowing for the docking of ligands, along with a good estimate of their affinities. The identification of this docking site was consistent with subsequent compound design and biological data obtained [10]. MD simulations validated the stability of docked poses and illustrated the role of solvent molecules in the binding pocket. The value of the described techniques lies in their ability to rapidly evaluate the potential of a new ligand in silico before spending precious time and resource on its synthesis and experimental evaluation.

Acknowledgments

JR Baker is the grateful recipient of an Australian Government Research Training Program (RTP) scholarship. This work was conducted in part under the auspices of a Fulbright Fellowship awarded to S Paula.

Conflict of interest

The authors declare no conflicts of interest.