Effect of Kinase Inhibiting RNase Attenuator (KIRA) Compounds on the Formation of Face-to-Face Dimers of Inositol-Requiring Enzyme 1: Insights from Computational Modeling

Abstract

Inositol-requiring enzyme 1α (IRE1α) is a transmembrane dual kinase/ribonuclease protein involved in propagation of the unfolded protein response (UPR). Inositol-requiring enzyme 1α is currently being explored as a potential drug target due to the growing evidence of its role in variety of disease conditions. Upon activation, IRE1 cleaves X-box binding protein 1 (XBP1) mRNA through its RNase domain. Small molecules targeting the kinase site are known to either increase or decrease RNase activity, but the allosteric relationship between the kinase and RNase domains of IRE1α is poorly understood. Subsets of IRE1 kinase inhibitors (known as “KIRA” compounds) bind to the ATP-binding site and allosterically impede the RNase activity. The KIRA compounds are able to regulate the RNase activity by stabilizing the monomeric form of IRE1α. In the present work, computational analysis, protein–protein and protein–ligand docking studies, and molecular dynamics simulations were applied to different IRE1 dimer systems to provide structural insights into the perturbation of IRE1 dimers by small molecules kinase inhibitors that regulate the RNase activity. By analyzing structural deviations, energetic components, and the number of hydrogen bonds in the interface region, we propose that the KIRA inhibitors act at an early stage of IRE1 activation by interfering with IRE1 face-to-face dimer formation thus disabling the activation of the RNase domain. This work sheds light on the mechanism of action of KIRA compounds and may assist in development of further compounds in, for example, cancer therapeutics. The work also provides information on the sequence of events and protein–protein interactions initiating the unfolded protein response.

1. Introduction

The accumulation of misfolded proteins in the endoplasmic reticulum (ER) triggers an evolutionarily conserved intracellular signaling pathway called the unfolded protein response (UPR) [1]. The UPR is mainly an adaptive response to re-establish ER proteostasis through the signaling of three transmembrane proteins, inositol-requiring protein 1α (IRE1α, hereafter IRE1), protein kinase R (PKR)-like ER kinase (PERK), and activating transcription factor 6 (ATF6) [2]. Inositol-requiring protein 1α represents the most evolutionarily conserved branch of the UPR; it is a multidomain protein with an N-terminal luminal domain and cytosolic kinase and RNase domains connected by a transmembrane linker [2]. Upon activation, IRE1 dimerizes, trans-autophosphorylates, and oligomerizes thereby activating the IRE1 RNase domain which is able to remove a 26 nucleotide intron from X-Box binding protein 1 (XBP1) mRNA. Upon re-ligation by RNA-splicing ligase RtcB, spliced XBP1 (XBP1s) is formed [2]. The XBP1s protein is a potent transcription factor, and its target genes in turn permit the ER to adapt to stress conditions [2]. There are several models that explain how ER stress is sensed by IRE1 [3]. One of the established models links its trans-autophosphorylation with face-to-face dimer orientation as the first step of IRE1 activation (Figure 1A) and subsequent transition to back-to-back dimer (Figure 1B) and larger oligomeric structure formation to stimulate the RNase domain activity [3].

IRE1 signaling is implicated in the etiology of various diseases [2,4], including cancer, where tumor cells activate the UPR to avoid apoptosis and survive [5]. Thus, IRE1 has been the recent focus of several drug discovery projects in cancer research [6]. Current IRE1 modulators can be categorized as (i) direct IRE1 RNase inhibitors [6,7]; (ii) ATP-competitive kinase inhibitors that activate RNase [8]; and (iii) ATP-competitive kinase inhibitors that decrease RNase activity [8,9], the latter also known as “kinase inhibiting RNase attenuators” (KIRAs). Optimized KIRA compounds include KIRA6 [10], KIRA7 [4], and KIRA8 [11]. Both KIRA6 and KIRA7 possess an imidazopyrazine scaffold [10], whereas KIRA8 is a sulfonamide compound with high selectivity [11].

The observation that KIRAs allosterically inhibit IRE1 RNase domain was confirmed by a competitive in vitro assay where ATP-competitive RNase activators were shown to completely restore the RNase activity in the presence of KIRA [8]. The model proposed by Feldman et al. [8] speculates that KIRAs stabilize the DFG-out kinase domain conformation and helix-αC displacement that makes it incompatible with back-to-back dimer formation thereby leading to kinase and RNase inhibition.

Herein, we propose an alternative and, we think, a more likely model to rationalize the link between kinase binding and RNase domain inhibition by the KIRA compounds based on molecular level analyses of structures and dynamics of IRE1 in complex with KIRA compound 16 (Figure 1C); numbering from the original paper [9] (hereafter referred to as KIRA).

Using protein–protein docking [12], protein–ligand docking [13], and molecular dynamics (MD) simulations [13], we investigated the structure and dynamics of the two relevant IRE1 dimer structures (Figure 1) face-to-face and back-to-back which are critical for trans-autophosphorylation and RNase activity, respectively [3]. Using models in the presence and absence of KIRA bound to the kinase pocket, we propose a mechanistically based model of how KIRAs inhibit the RNase activity of IRE1 by allosteric interaction with the kinase domain.

2. Methods

2.1. Selection and Preparation of IRE1 Crystal Structure

The IRE1 crystal structures with PDB codes 4U6R (KIRA-bound monomer), 3P23 (ADP-bound face-to-face dimer), and 4YZC (staurosporine-bound back-to-back dimer) were prepared using the Schrödinger protein preparation wizard [14]. Using the 3P23 PDB structure, the face-to-face dimer was obtained by deleting Chain C and D and generating missing loops using Prime [15]. Hydrogen atoms were added and the protonation and tautomeric states of Asp, Glu, Arg, Lys, and His were adjusted to match a pH of 7.4. Possible orientations of Asn and Gln residues were generated. Finally, the IRE1 dimer and monomer structures were subjected to geometry refinements using the OPLS3 force field [16] in restrained minimizations.

2.2. Protein-Protein Docking

A protein–protein docking analysis was performed to understand the structural basis of IRE1 recognition and regulation mediated by KIRA. Initially, five different protein–protein docking programs were chosen: SwarmDock [17], ZDOCK [18], HsymDock [19], PatchDock [20], and ClusPro [21]. The ability of the programs to reproduce the native IRE1 dimers (PDB code 4YZC for back-to-back and 3P23 for face-to-face) was checked. Starting from the known crystallographic dimer structures of IRE1 in back-to-back (PDB code: 4YZC) and face-to-face (3P23) conformers, we split the crystal structure into monomers and explored the capabilities of the software to reproduce the two dimer complexes. The best performing program, SwarmDock [17], was then used for predicting the structures of KIRA-bound dimers using the 4U6R PDB monomer structure. A schematic representation of the method used is shown in Figure 2.

2.3. KIRA Preparation for Docking Studies

The co-crystallized ligand from PDB structure 4U6R (KIRA) shown in Figure 1C, was extracted and used for docking studies. The KIRA compound was prepared using LigPrep [22] in the Schrödinger suite [23]. The OPLS3 force field [16] was used for KIRA preparation steps and possible protonation and ionization states were assigned using Ionizer at pH 7.4.

2.4. Molecular Docking of KIRA

KIRA was docked in the kinase pocket of the IRE1 dimer back-to-back (PDB code: 4YZC) and face-to-face (PDB code: 3P23) structures using the Glide program [24] in Schrödinger [23] with the receptor grids prepared using the OPLS3 force field [16]. The molecule was docked in the kinase domain of both monomers of each dimer. The IRE1 dimers in apo form were obtained by deleting the small organic molecules and ions present in the crystal structures, from the kinase active site (i.e., ADP and Mg²⁺ in the 3P23 PDB structure and staurosporine in the 4YZC PDB structure, respectively). Two grid centers per dimer (i.e., one per each kinase active site) were prepared, each cubic grid with a side length of 20 Å. The grid center was set at the centroid of Lys599, a residue crucial for the kinase activity. The XP (extra precision) docking mode and flexible ligand sampling were employed in the docking procedure. All other parameters were set to default values.

2.5. Molecular Dynamics Simulations

The stability of the native crystallographic dimer structures of the IRE1 in back-to-back (PDB 4YZC) and face-to-face (PDB 3P23) conformer were compared with the predicted KIRA-bound dimer forms, based on 300 ns MD simulations. For KIRA bound dimers, we used both complexes with KIRA docked in the kinase active site of the dimer structures 4YZC and 3P23 and dimer structures generated from KIRA-co-crystallized IRE1 monomers.

For the MD simulations, the following steps were followed:

• Systems preparation: systems included the experimental IRE1 dimer structures (PDB 4YZC, 3P23), predicted dimers (from KIRA-bound monomer) (Section 2.2.), and KIRA-docked dimer forms (Section 2.4.). The systems were prepared separately as discussed in Section 2.1.
• Ligand parameterization: ligands (ADP, staurosporine, and KIRA) were parametrized with GAFF as implemented in Ambertools2018 using the Antechamber interface tool [25]. The AM1-BCC atomic point charges [26] were calculated using Antechamber [27].
• Molecular dynamics simulation protocol: MD simulations were performed using the GROMACS 5.1 package [28] with the AMBER14SB force field for the protein [29]. The systems were explicitly solvated using cubic water boxes with cell borders placed at least 10 Å away from the protein or ligand atoms using TIP3P water [30] under periodic boundary conditions. The rational for the choice of the 10 Å cutoff distances was to place the protein or ligand atoms at a distance longer than the non-bonded interactions cut-off (i.e., 8 Å). The systems were first neutralized and Na⁺/Cl^– counter ions were added to give a physiological salt concentration of 0.154 M. All simulation runs consisted of energy minimization until the force was less than 1000 kJ mol⁻¹ nm⁻¹, 200 ps under NVT conditions subjected to position-restrained equilibration on the heavy atoms of IRE1, snf 200 ps equilibration and 300 ns of classical molecular dynamics simulation under NPT conditions. The simulations were run in triplicate (referred to as Replica 1, 2, and 3). In all simulations, the temperature was kept at 300 K by the velocity rescaling thermostat [31] with a coupling constant of 0.1 ps and pressure at 1.01325 bar using the Parrinello–Rahman barostat [32] with a coupling time of 5.0 ps, excluding NVT pre-simulation steps. Constraints were applied on all bonds using the LINCS algorithm [33]. The leap-frog algorithm [34] was employed in the simulations with integration timesteps of 2 fs.

The structural deviations during the MD simulation were analyzed using RMSD, number of distinct hydrogen bonds, and energy terms such as electrostatic (Ele) and van der Waals (vdW) interactions using built-in tools in the GROMACS 5.1 package [28]. For analyzing the dimer interface RMSD, an index file was created with specific residues. The dimer interface was defined as any pair of Cα atoms from one monomer within 10 Å of the other in the face-to-face and back-to-back dimers [35,36].

2.6. Data Availability

All simulation protocols are provided as tarballs (.tar.gz) freely accessible at zenodo.org as DOI: 10.5281/zenodo.3368654.

There are eighteen tarball (.tar.gz) files, three replicas for each of the systems investigated.

The contents of each tarball is as follows:

• a source PDB (.pdb) file
• leap.log—commands used to create the. prmtop and. inpcrd files
• two AMBER parameter/topology (.prmtop) and an AMBER coordinate (.inpcrd) file
• .mdp file used for performing all the minimisation, relaxation, equilibration, and production run steps
• Executable script (i.e., job009) that was used to perform the production run
• trajectory (.xtc) files for each independent MD simulation

3. Results and Discussion

3.1. Protein–Ligand Docking Analysis

KIRA-bound IRE1 dimer crystal structures are not available. To predict the KIRA binding in dimers, the compound was docked in the IRE1 back-to-back (PDB code: 4YZC) and face-to-face (3P23) structures. Visual inspection of the docked poses revealed a different binding mode of KIRA compared to the 4U6R PDB structure (Figure S1). In particular, in the face-to-face dimer, KIRA binding was mainly stabilized by electrostatic interaction with Asp711, Asp688, and Lys690 (Figure S2A,B) and in the back-to-back dimer with Glu651 and Cys645 (Figure S2C,D), while the co-crystallized KIRA interacted mainly with Lys599, Glu651, Cys645, Phe712, and Ile642 as reported in previous studies [9,37].

3.2. Protein–Protein Docking Analysis

To address questions regarding the effect of KIRA binding on the IRE1 dimer formation, and if KIRA is able to structurally interfere with either the face-to-face or back-to-back dimer form, or both, we assessed if the KIRA-bound monomer structure (PDB code: 4U6R) was capable of forming dimer structures.

In order to identify an appropriate protein–protein docking program, a series of docking experiments were performed using five freely available blind protein–protein docking programs. In the absence of experimental data (to the best of our knowledge) which could suggest interface residues critical for IRE1–IRE1 interaction, we considered the blind protein–protein docking approach best suitable for our scope. Starting from the known crystallographic dimer structures of IRE1 in back-to-back (PDB code: 4YZC) and face-to-face (3P23) conformers, we split the crystal structure into monomers and tried to reproduce the dimer complexes with the programs. The docking results are shown in

Table 1. The RMSD^a results (in angstrom) of five different protein–protein docking approaches to reproduce the known IRE1 dimer complexes.

	Face-to-Face Dimer (PDB Code: 3P23)	Back-to-Back Dimer (PDB Code: 4YZC)
SwarmDock	1.39	3.56
ZDOCK	12.48	3.32
HsymDock	3.12	13.25
PatchDock	24.33	29.49
ClusPro	3.58	31.01

^a Root mean square deviation (RMSD) was calculated for Cα atoms by superimposing the top-scored docked pose generated by the program with the crystallographic structures.

and Table S1.

After evaluating the RMSD over the five top-scored docked poses for all programs studied here (Table S1), we note that the top-scored docking pose in almost all cases was also the one with the lowest RMSD. As observed in a recent benchmark study on the accuracy of free protein–protein docking methods [12], it is not surprising that the five docking programs displayed large variations in the results. Of the five docking programs tested in the current study, SwarmDock was able to reproduce the native crystallographic back-to-back (PDB code: 4YZC) and face-to-face (PDB code: 3P23) dimers structures of the IRE1 to an RMSD of 3.56 and 1.39 Å, respectively, as shown in Figure 3. Based on acceptable quality levels in reproducing the near-bound dimer states (Figure 3), SwarmDock was thus also used for predicting dimers of the KIRA-bound monomer structures (PDB 4U6R) in the back-to-back and face-to-face orientations.

We first evaluated steric clashes at the interchain region by superposing the KIRA-bound monomer (PDB ID 4U6R) onto each monomer of the native crystallographic structures of the IRE1 back-to-back and face-to-face dimers, respectively (Figure S3). The KIRA-bound dimer forms thus generated produced several steric clashes at the interchain level, especially at the helix-αC and the activation segment of the KIRA-bound face-to-face dimer (Figure S3). In particular, the backbone of Glu604 and the carboxylate group of Glu735 created steric clashes with the guanidino group of Arg600 and the carboxylate of Glu735, located in the helix-αC and the activation segment, respectively (Figure S3). In the back-to-back dimer, the guanidine groups of Arg627 and Arg905 created potential steric clashes with the guanidine groups of Arg627 and Arg905 located in the other monomer of the symmetric complex (Figure S3).

In order to generate more appropriate KIRA-bound dimer forms (either IRE1 back-to-back or face-to-face), we performed protein–protein docking using the 4U6R PDB structure. The monomer of the 4U6R PDB structure was first superposed on each monomer of the native crystallographic structures of the IRE1 in back-to-back (PDB code: 4YZC) and face-to-face (PDB code: 3P23) forms in order to generate 3D coordinates for the docking input. Results from the protein–protein docking calculations of the KIRA-bound dimer forms performed using SwarmDock are shown in Figure 4. The best docked results were analyzed further through MD simulations.

The best-scoring docked pose (in terms of Cα atom RMSD) generated for the face-to-face dimer shown in Figure 4A was very far from the experimental one with an RMSD of 24.24 Å. No structures with a better RMSD were found within the 5 top-scoring docking poses (RMSD 24.24, 25.29, 34.11, 27.49, and 25.40 Å, respectively). Since the program was shown to successfully predict the crystallographic dimer forms, the high RMSD values were not likely a result of bad sampling. Rather, the KIRA-induced conformational changes rendered the system incapable of appropriate dimer formation as supported by our analysis of the steric clashes at the interchain level (Figure S3A).

A comparison of the best back-to-back docking poses obtained using the 4U6R PDB structure and the crystal structures of the back-to-back dimer was also performed (Figure 4B). The RMSD analysis for the best docking pose gave a value near 4 Å, in line with the RMSD value identified when using the native monomer structure (Table 1).

3.3. MD Simulations Analysis: Influence of KIRA on the Face-to-Face Dimer

Three different IRE1 face-to-face dimers were explored further, namely, the native crystal dimer structure (PDB code: 3P23), the native structure (PDB 3P23) with KIRA docked in the active sites, and the protein–protein docked pose of PDB 4U6R in the face-to-face dimer. The stability of each of the three systems was analyzed during three MD replicas, each replica being of a 300 ns duration.

To assess the structural stabilities, RMSD values of each dimer and of the dimer interfaces were calculated (Figure 5). For each replica, the calculations were done by considering the structures present in the minimized, equilibrated systems as the reference points. The RMSD values were analyzed as the functions of simulation time. The three replicas for the native face-to-face crystal dimer (PDB code: 3P23) revealed that the IRE1 dimer was stable as evidenced by low and relatively constant RMSD values of the three independent trajectories (Figure 5A,D). For the KIRA-docked face-to-face dimer, slightly higher values of RMSD (Figure 5B,E) were seen, supporting our first speculation of the effects exerted by KIRA on the formation of IRE1 face-to-face dimers. In contrast, the predicted PDB 4U6R dimer (Figure 4A) displayed higher RMSD values during the MD simulations (Figure 5C,F) even though it already started from a very large RMSD value (Figure 4A), indicating that this dimer explores even more distorted structures compared to the native IRE1 face-to-face dimer. Conformations with large interface RMSD were visually examined and, by comparing the distance of the center of mass (COM) of each RNase domain in each dimer with the native face-to-face dimer, we analyzed the impact of KIRA on the stability of the system (Figure S4). The COM distances of the RNase domains in the PDB 3P23 dimer with the KIRA docked, and in the protein–protein docked pose of PDB 4U6R in the face-to-face dimer form, were significantly higher compared to that of the native IRE1 dimer (Figure S4). This data further corroborates the impact of KIRA on the stability of the system.

To further analyze the system deviation, the interaction energy and the number of H-bonds between the monomers were computed (Figure 6). Interaction energy analysis among the monomers revealed a smaller energetic stabilization of the KIRA-docked face-to-face dimer and the protein–protein docked pose of PDB 4U6R in the face-to-face dimer form compared to the native face-to-face crystal dimer (PDB code: 3P23) (Figure 6). The same trend was confirmed by the H-bond analysis with a higher number of H-bonds occurring in the native IRE1 face-to-face dimer (average number of H-bonds over the simulation was 14) compared to the protein–protein docked pose and KIRA-docked face-to-face dimer, respectively, with both having an average number of 11 H-bonds (Figure 6). The overall analysis confirms our initial interpretation of the different stabilities of the three IRE1 face-to-face dimers investigated, with the native one being the more stable, indirectly reflecting the impact of KIRA on the stabilization of the IRE1 face-to-face dimer.

Moreover, energetic analysis of KIRA and ADP in each kinase-active site of the face-to-face dimer revealed a lower energetic stabilization of ADP compared to KIRA, with each ligand being able to interact favorably with the IRE1 active site pocket for the entire simulation time during all three replicas (Figure S5).

3.4. MD Simulations Analysis: Influence of KIRA on the Back-to-Back Dimer

To investigate the impact of KIRA on the IRE1 back-to-back dimer, three different systems were also considered here, namely, the native back-to-back dimer crystal structure (PDB code: 4YZC), the native dimer structure (PDB code: 4YZC) with KIRA docked, and the protein–protein docked pose of PDB 4U6R in the back-to-back dimer form, respectively. The stabilities of the three systems were studied during three MD replicas, each replica being 300 ns in length.

As reported in Section 3.2, RMSD values of each dimer and the RMSD for the IRE1 dimer interface (Figure 7 and Figure S6, respectively) were calculated to assess the structural stability of the dimer structures. For each replica, the structure present in the minimized, equilibrated system was used as the reference point and RMSD values were analyzed as functions of simulation time. The three replicas for the native back-to-back crystal dimer (PDB code: 4YZC) revealed that the IRE1 dimer was stable as evidenced by the low and relatively constant RMSD values of the three independent trajectories (Figure 7A and Figure S6A). In contrast to the face-to-face dimer systems investigated, the KIRA-docked back-to-back dimer (PDB code: 4YZC) and the protein–protein docked pose of PDB 4U6R in the back-to-back dimer form showed similar structural stabilities as the native IRE1 back-to-back dimer (Figure 7B,C and Figure S6B,C).

Also, for the three back-to-back systems, the interaction energy and the number of H-bonds among the monomers were analyzed (Figure 8). Energetic analysis of the occurring interchain interactions shows similar energetic stabilization of all three back-to-back dimers structures and in all three replicas (Figure 8). The same trend was confirmed by the H-bond analysis with an equal number of H-bonds occurring between the native IRE1 back-to-back dimer compared to the protein–protein docked pose and KIRA-docked back-to-back dimer (Figure 8). The overall analysis of these three IRE back-to-back dimers confirms our initial interpretation of the similar stabilities of the three dimer models. To further validate this hypothesis, the last frame of each MD simulation was superposed with the native back-to-back dimer revealing an overall similar IRE1 back-to-back active conformation (Figure S7).

Moreover, energetic analysis of KIRA and staurosporine in each kinase active site of the back-to-back dimer revealed similar energetic stabilization of staurosporine compared to KIRA with each ligand being able to interact favorably with the IRE1 pocket active site during the entire simulation in the three replicas (Figure S8).

Finally, inspection of the number of H-bonds between the face-to-face and back-to-back dimers (Figure 6 and Figure 8) suggests that the back-to-back dimer was stabilized by a much larger number of H-bonds (at least two-fold more). This correlates with the lower structural fluctuation (i.e., RMSD analysis) of the back-to-back compared to the face-to-face dimer.

4. Conclusions and Perspective

We investigated the impact of KIRA binding on IRE1 dimer structures. Unexpectedly, the docking and MD simulations studies revealed that KIRA can bind to the kinase pocket of IRE1 in both the native face-to-face and back-to-back forms. A detailed analysis of the IRE1 monomer–monomer interactions process for the face-to-face dimer in the presence of KIRA revealed energetic destabilization, suggesting that the binding of KIRA already affects the system at the stage of face-to-face dimer formation. Given that IRE1 activation appears to be dependent on the close communication between the kinase and RNase domains, the data lead us to believe that KIRA has a prominent role at the early stage of IRE1 activation by destabilizing face-to-face dimer formation. This will impair the trans-autophosphorylation process thus preventing IRE1 from reaching the RNase active back-to-back structure.

The proposed mechanism of blocking the trans-autophosphorylation provides a molecular level validation of available experimental data where KIRA compounds inhibit IRE1 phosphorylation [8]. This is further supported by the experimental observation that, upon the inclusion of KIRAs, Western blotting reveals the formation of IRE1 monomers only [8].

Moreover, it is worthwhile to underline that KIRA-bound IRE1 dimer crystal structures are not available. Currently, one of our major efforts is to consider protein structural variations during KIRA molecular docking and explore the importance of the activation segment during the dimerization process to complement the aforementioned results with a detailed characterization and rationalization of the IRE1 dimerization impairment with KIRA compounds.

The data reported herein provide another small piece of information towards the understanding of IRE1 activity and the structural evidence of KIRA’s role in the IRE1 inhibition process, representing a stimulus to explore and better understand IRE1 signaling and accelerate the IRE1 drug design campaign.