J Mol Model (2015) 21:167 DOI 10.1007/s00894-015-2708-z
ORIGINAL PAPER
Characterization of interactions and pharmacophore development for DFG-out inhibitors to RET tyrosine kinase Chunxia Gao 1 & Morten Grøtli 1 & Leif A. Eriksson 1
Received: 4 February 2015 / Accepted: 19 May 2015 # Springer-Verlag Berlin Heidelberg 2015
Abstract RET (rearranged during transfection) tyrosine kinase is a promising target for several human cancers. Abt348, Birb-796, Motesanib and Sorafenib are DFG-out multikinase inhibitors that have been reported to inhibit RET activity with good IC50 values. Although the DFG-out conformation has attracted great interest in the design of type II inhibitors, the structural requirements for binding to the RET DFG-out conformation remains unclear. Herein, the DFG-out conformation of RET was determined by homology modelling, the four inhibitors were docked, and the binding modes investigated by molecular dynamics simulation. Binding free energies were calculated using the molecular mechanics/Poisson-Bolzmann surface area (MM/PBSA) method. The trends in predicted binding free affinities correlated well with experimental data and were used to explain the activity difference of the studied inhibitors. Per-residue energy decomposition analyses provided further information on specific interaction properties. Finally, we also conducted a detailed e-pharmacophore modelling of the different RET-inhibitor complexes, explaining the common and specific pharmacophore features of the different complexes. The results reported herein will be useful in future rational design of novel DFG-out RET inhibitors. Keywords RET . DFG-out inhibitors . Molecular dynamics simulation . MM-PB(GB)SA . e-pharmacophore Electronic supplementary material The online version of this article (doi:10.1007/s00894-015-2708-z) contains supplementary material, which is available to authorized users. * Leif A. Eriksson [email protected] 1
Department of Chemistry and Molecular Biology, University of Gothenburg, 405 30 Göteborg, Sweden
Introduction RET (rearranged during transfection) is a transmembrane tyrosine kinase receptor essential for development of the sympathetic, parasympathetic and enteric nervous systems, and of the kidneys [1]. Furthermore, it is well established that deregulation and dysfunction of RET can lead to several different thyroid cancers, including multiple endocrine neoplasia 2 syndromes (MEN2A and MEN2B), familial medullary thyroid carcinoma (FMTC), and papillary thyroid carcinoma (PTC) [1]. Point mutations in the RET gene are the cause of the three cancer syndromes MEN2A, MEN2B, and FMTC, while PTC is related to chromosomal inversions or translocations. Furthermore, recent studies have demonstrated that RET is overexpressed in a subset of estrogen-receptor (ER)positive breast cancers and that crosstalk between ER and RET is important in responses to endocrine therapy [2, 3]. The tyrosine kinase domain of RET is comprised of two lobes: a smaller N-terminal lobe consisting mainly of β-strands and a single large α-helix; and a larger C-terminal lobe, which is almost exclusively α-helical. The ATP-binding site is located at the hinge region between the two lobes. The structure of the RET tyrosine kinase domain can be further divided into a number of structural regions; these include the glycine-rich loop (G-loop), the C α-helix on the Nlobe, the activation loop (A-loop), the catalytic loop (C-loop) and conserved DFG motifs on the C-lobe (Fig. 1a). In all known protein kinase structures, the DFG motif assumes a conformation with the Phe residue buried in a hydrophobic pocket in the groove between the lobes (DFG-in conformation). Recently, much attention has been devoted to the use of small molecule inhibitors to target RET in MEN2-related tumors; e.g., pyrazolopyrimidines PP1 and PP2, indolocarbazoles, 2indolinone derivatives, and the 4-anilinoquinazoline ZD6474 have been shown to display strong activity towards RET kinase
167
J Mol Model (2015) 21:167
Page 2 of 12
Fig. 1 a Ribbon representations of DFG-out RET (rearranged during transfection) tyrosine kinase structures. Key secondary structural elements are colored green (glycine-rich loop), red (Cα-helix), blue (activation loop) and pink (catalytic loop). b Alignment of DFG-in and DFGout residues. c Key residues of RET interacting with the inhibitors. d–h Illustration of RETDFG out inhibitor interaction: d Abt348, e Birb-796, f Motesanib, g Sorafenib-I, h Sorafenib-II
A
C
B
D
E
F
G
H
[4–6]. These inhibitors are all ATP competitive and bind exclusively within the ATP binding site of the active DFG-in kinase conformation; so-called type I inhibitors. Since the kinase catalytic domain is highly conserved and the ATP binding sites are similar, type I inhibitors are not specific to RET alone, but exhibit multi-kinase targeting. Instead of utilizing the allosteric pocket, which is much less conserved, a second class of inhibitors has been developed that act by stabilizing the inactive DFG-out kinase conformation [7]. This class, known as type II inhibitors, appear to display much higher selectivity. The DFG-out conformation has thus attracted significant interest in the discovery and development of selective kinase inhibitors in the past few years [8–14]. However, a recent analysis of binding modes of type II
inhibitors indicates that there may be less selectivity than initially believed [15]. It is thus important to investigate the interactions and functions of this type of inhibitor to better understand the key interactions and how selectivity may be enhanced further. Although not developed as RET inhibitors per se, several type II inhibitors have proven to be highly potent RET inhibitors, for example Sorafenib (IC50 =50 nM) [16] and Motesanib (IC50 =59 nM) [17] have both have reached clinical Phase II trials for treatment of MTC. Besides targeting RET, the latter drugs also mainly target vascular endothelial growth factor receptors (VEGFRs), platelet-derived growth factor receptor (PDGFR), and the KIT receptor protein-tyrosine kinase [18]. Two more compounds, Birb-796, a highly potent p38α inhibitor for anti-inflammatory therapy, and Abt-348, a novel
J Mol Model (2015) 21:167
Page 3 of 12 167
aurora kinase inhibitor, have also been tested and show IC50 values towards RET of 30 nM [19] and 7 nM [20], respectively (Table 1). To summarize, the four compounds cited above have been shown to have high experimental IC50 values to RET, and a study of how those compounds interact with RET is crucial for the development of more potent and selective RET inhibitors. Currently, X-ray crystal structures are available only for compounds binding to targets other than RET; i.e., Sorafenib binding to VEGFR2 (PDB code: 4ASD) [21], Motesanib binding to VEGFR2 (PDB code: 3EFL), and Birb-796 binding to p38α (PDB code: 1KV2) [13]. As there are no NMR or X-ray crystal structures of these compounds binding to RET, detailed information about the RET-type II inhibitor interaction is lacking. Such information is urgently needed in order to further improve our understanding of the structural determinants of the type II kinase inhibitor in binding and stabilizing the DFG-out conformation of RET. This will assist towards the rational design of novel, potent and selective RET kinase inhibitors. In this study, the DFG-out conformation of RET was obtained by homology modelling using VEGFR2 as template. After docking of the type-II inhibitors to the DFG-out conformation of RET, molecular dynamics (MD) simulations combined with free energy calculations were carried out to obtain
Table 1
information about binding of the inhibitors. Absolute binding free energy calculations were performed using the molecular mechanics Poisson-Bolzmann surface area (MM-PBSA) method [22]. Detailed binding free energies between the inhibitors and RET residues were calculated using a per-residue molecular mechanics/generalized Born surface area (MM-GBSA) decomposition analysis [23]. Furthermore, based on the optimized MD complex structures of RET and its inhibitors, multiple energy-based pharmacophore (e-pharmacopohore) models were constructed. We thus investigated features of the interaction between RET and type II inhibitors based on dynamic structural information from MD simulations, and outlined their specificity by analyzing the contribution of such interactions to the binding free energy of RET complexes, and their key pharmacophore descriptors.
Methods Homology model of DFG-out RET The tyrosine kinase domain (from Phe 719 to Lys 1011) of RET in its DFG-out conformation was obtained by homology
Structures and IC50 values of different inhibitors
NAME
IC50 (nM)
SORAFENIB
50 (Ref 16)
MOTESANIB
59 (Ref 17)
ABT-348
7 (Ref 20)
BIRB-796
30 (Ref 19)
STRUCTURE
167
J Mol Model (2015) 21:167
Page 4 of 12
modeling. The protein sequence of RET was obtained from the NCBI database (GeneID: 547807). The crystal structure of VEGFR2 (PDB code 4ASD, 2.03 Å resolution) [21] was chosen as the template for homology modeling calculations based on the fact that the tyrosine domain of VEGFR and RET share 47 % similarity, and that the DFG-out loop is present in its entirety in this crystal structure. A previous study dealing with construction and validation of a RET tyrosine kinase catalytic domain through homology modelling claimed that the Bopen^ and Bclosed^ states were both modelled [24]. However, both states reported in the study have the same DFG-in conformation; hence homology modelling of the RET DFT-out conformation in our study is not analogous to previous work. Intermediate homology models (25 backbone models×50 side chain models; in total 1250) were generated, and the bestscoring model as determined by the Generalized Born/ Volume Integral (GB/VI) scoring method [25], which compares the electrostatic solvation energy, was taken to be the final model. For the homology modelling, the MOE 2013.08 programme was used [26]. A Ramachandran plot of the final model (Supplementary material, Figure S1) shows that 98.4 % of the residues lie in the allowed regions, and that the 1.6 % outliers are in areas that are not in contact with the active site or affected by ligand binding. Docking The homology model of RET was superposed with the crystal structure of VEGFR2 in complex with Sorafenib (PDB code 4ASD) [21], after which the VEGFR2 structure was deleted and Sorafenib retained in the binding pocket of RET. The resulting RET–Sorafenib complex was prepared and refined using the default protein preparation wizard in Schrodinger. After ensuring chemical correctness, hydrogens were added, and side chains far from the binding cavity and not participating in salt bridges were neutralized. The hydrogen bonding network was optimized by reorienting the hydroxyl groups, water molecules, and amide groups of Asn and Gln, and selecting appropriate states and orientations of the imidazole ring in His residues. Water molecules present in the crystal structure were deleted, and termini were capped by adding ACE and NMA residues. Missing side chains and loops were added. Subsequently, the structures were minimized using the OPLS-2005 force field [27] and the Impact molecular mechanics engine, while heavy atoms were constrained. The prepared protein structure was then used for subsequent grid generation. Grids were generated by the Receptor Grid Generation panel (Schrodinger, LLC, New York, NY) which defines the receptor structure by excluding any co-crystallized ligand present, determines the position and size of the active site to be represented by receptor grids, and sets up Glide constraints. Grids were defined by centering them on the ligand using the
default box size, and with H-bonds to Glu775 and Asp892 set to be constrained. Molecular docking experiments with Abt-348, Birb-796, Motasenib and Sorafenib were carried out using Glide, as implemented in the Schrodinger package. The XP (extra precision) [28] docking function of Glide was used, granting full flexibility to the ligands. A post-docking minimization was performed on the output complexes in order to reduce the initially collected 10,000 poses ligand to five. The docked poses were then minimized using the local optimization feature in Prime. For Sorafenib, two different binding modes was observed in the different crystal structures (PDB code: 3HEG, 3GCS) [29, 30]; hence, two different docking poses were selected for the following studies. To validate the docking approach, Birb-796, Motasenib and Sorafenib (two structures) were docked back into the crystal structures of their respective receptors (Supplementary material, Figure S2). As can be seen from these benchmark dockings, the molecules overlap almost perfectly to their crystal structures, with very low RMSD values. MD simulations Atomic charges on Abt-348, Birb-796, Motesanib and Sorafenib were calculated using the restrained electrostatic potential (RESP) procedure at the HF/6-31G* level, after minimizing the molecule at the AM1 semiempirical level [31]. GAFF force field parameters [32] and the RESP partial charges were assigned using the Antechamber module in the AMBER10 package [33]. All missing hydrogen atoms were added using the LEaP module. Five conventional 100 ns MD simulations were performed on the RET-inhibitor complex structures using the Amber 99 force field [34, 35] and the GROMACS software [36]. The structures were solvated in periodic boxes with a buffer distance of 10.0 Å. A number of Na+ and Cl− ions were added to satisfy the electroneutrality condition and to give a salt concentration of 0.1 mol L−1, using the genion module in GROMACS. The obtained systems were energy minimized by steepest descent (200 steps) to remove close contacts. Position restrained simulations (2 ns duration, 1.0 fs time step, NPT ensemble, T=298 K, P=1 bar) were first performed, to enable the water molecules to reach more favorable positions. Particle-mesh Ewald (PME) [37, 38] summation was used for long-range electrostatics. A 12 Å cutoff was used for both Coulomb and Lennard-Jones interactions. The temperature and pressure was controlled through the Berendsen coupling algorithm [39], with the time constants 0.1 ps for temperature and 1.0 ps for pressure coupling. All bond lengths were constrained using the LINCS algorithm [40]. During the production simulations (200 ns duration, 1.0 fs time step, NPT ensemble, T=298 K, P=1 bar), the temperature was controlled using the Nose-Hoover thermostat [41] with a time
J Mol Model (2015) 21:167
Page 5 of 12 167
constant 0.1 ps, and the pressure was controlled using the Parrinello-Rahman barostat [42] with a time constant 1.0 ps. The remaining parameters were the same as in the position restrained simulations.
10−4 kcal mol−1 Å−2. Due to the high computational demand of this approach, only 50 snapshots that were evenly extracted from the last 2 ns MD trajectory were used to calculate the entropic contribution.
MM-PBSA calculations
Residue-inhibitor interaction decomposition
The binding free energies were calculated using the MMPBSA method implemented in AMBER10 [33]. For each complex, a total number of 500 snapshots were taken from the last 2 ns on the MD trajectory, and the free energy computed for each molecular species (complex, protein, and ligand), giving the binding free energy as the difference ΔΔGbind ¼ ΔGcomplex – ΔGprotein þ ΔGligand ð1Þ
To identify the key residues responsible for the binding of the inhibitors, free-energy decomposition was performed for each residue. Given the pairwise nature of the GB equation allowing the decomposition of polar solvation energy (ΔGele,sol) into atomic contributions in a straightforward manner, the interaction between the inhibitor and each residue was computed using the MM-GBSA decomposition process in AMBER10 [3, 23]. The binding interaction of each inhibitor–residue pair includes four terms: van der Waals contribution (ΔGvdW), electrostatic contribution (ΔGele), polar solvation contribution (ΔGele,sol), and nonpolar solvation contribution (ΔGnonpol,sol):
Each term can be estimated as follows: ΔG ¼ ΔGMM þ ΔGsol −TΔS
ð2Þ
where ΔGMM is the molecular mechanics free energy, ΔGsol is the solvation free energy, and TΔS represents the entropy term. The molecular mechanics energy was calculated from the electrostatic and van der Waals interactions, ΔGMM ¼ ΔGele þ ΔGvdW
ð3Þ
and the solvation free energy from the polar and nonpolar contributions, ΔGsol ¼ ΔGele;sol þ ΔGnonpol;sol
ð4Þ
The polar part ΔGele,sol was obtained by solving the Poisson-Boltzmann (PB) equations [43] for the MM-PBSA method, or the generalized Born (GB) model [43] for the MM-GBSA calculations. Dielectric constants of 1.0 and 80.0 were used for solute and solvent, respectively. The parameters used for the GB calculation were as developed by Onufriev et al. [44] (GBOBC, igb=2). The nonpolar solvation contribution was determined using the solvent accessible surface area (SASA, Å2) with the Amber molsurf module, according to ΔGnonpol;sol ¼ γSASA þ b
ð5Þ
where γ, representing the surface tension, and the constant b were set to 0.0072 kcal mol−1 Å−2 and 0, respectively. The conformational entropy upon ligand binding (translation, rotation, and vibration), TΔS, was calculated using normalmode analysis with the NMODE program in AMBER10 [32]. Prior to the normal-mode calculations, each snapshot was minimized for a maximum number of 100,000 cycles using the conjugate gradient method with in the presence of a distance-dependent dielectric of 4rij (rij being the distance between two atoms I and j) until the root-mean-square value of the elements of the gradient vector were less than 1.0×
ΔGinhibitor‐residue ¼ ΔGele þ ΔGvdW þ ΔGele;sol þ ΔGnonpol;sol
ð6Þ
The vdW and electrostatic interactions between inhibitor and each residue of RET were computed using the Sander program in AMBER10 [33]. The polar solvation contribution (ΔGele,sol) was calculated using the GB model [43] and the nonpolar contribution of desolvation (ΔGnonpol,sol) was computed based on SASA. All energy components were calculated using 500 snapshots extracted from the last 2 ns of each MD trajectory. e-pharmacophore hypothesis generation The structure-based pharmacophore was generated using the energy-minimized final snapshot of each simulated RETinhibitor complex. The Receptor Grid Generation tool in the Maestro software package was used to generate energy grids for all prepared protein structures and the ‘refine only’ option was used in the XP [27] docking in Glide 4.2 X. Default settings were employed to minimize and optimize the structures. On the basis of the XP descriptor information, pharmacophore features were generated using PHASE 3.4. Pharmacophore sites were generated automatically with Phase v3.0 (Schrodinger) using the default set of six chemical features: hydrogen-bond acceptor (A), hydrogen-bond donor (D), hydrophobe (H), negative ionizable (N), positive ionizable (P), and aromatic ring (R)). Hydrogen-bond acceptor sites were represented as vectors along the hydrogen-bond axis in accordance with the hybridization of the acceptor atom. Hydrogen-bond donors were represented as projected points, located at the corresponding hydrogen-bond acceptor positions in the binding site. Projected points allow the possibility
167
J Mol Model (2015) 21:167
Page 6 of 12
for structurally dissimilar active compounds to form hydrogen bonds to the same location, regardless of their point of origin and directionality. Each pharmacophore feature site was first assigned an energetic value equal to the sum of the Glide XP contributions of the atoms comprising the site, allowing sites to be quantified and ranked on the basis of the energetic terms The Glide XP descriptors include terms for hydrophobic enclosure, hydrophobically packed correlated hydrogen bonds, electrostatic rewards, π–π stacking, π–cation, and other interactions. ChemScore, hydrogen bonding, and lipophilic atom pair interaction terms were also included when the Glide XP terms for hydrogen binding and hydrophobic enclosure were zero.
Results and discussion Static structural analysis of inhibitor binding mode A large conformational change in the residues in the conserved DFG motif in the active site of the kinase is required for the binding of the DFG-out inhibitors. In the RET DFGout inhibitor complexes, the Phe side chain has moved by ~10 Å to a new position (DFG-out conformation). This movement of the Phe side chain reveals a large hydrophobic pocket in the kinase (Fig. 1b). The complexes of the four inhibitors with RET reveal that the inhibitors share the same binding mode, with RET kinase in the DFG-out conformation. The inhibitors are stabilized in a cavity formed by residues Lys728, Leu730, Val738, Ala756, Lys758, Glu775, Val778, Leu779, Ile788, Leu802, Val804, Glu805, Tyr806, Ala807, Lys808, Tyr809, Gly810, Val871, His872, Leu881, Ser891, Asp892, Phe893 (Fig. 1c). Abt-348 has the highest molecular weight among the inhibitors and forms the most interactions with RET. The pyrazol group is positioned such that the hydroxyethyl group forms two hydrogen bonds, one with the Tyr806 side chain and the other one with Lys808 backbone oxygen. The adjacent thienol group interacts with the hinge region via two hydrogen bonds, between the pyrimidine N and the backbone amino nitrogen of Ala807, and between the amino NH2-group and the backbone amino oxygen of Glu805. The urea fragment forms three H-bonds with the RET protein: the oxygen interacts with the backbone NH of the Asp892 side chain (in DFG) and both nitrogens interact with the side chain of the highly conserved residue Glu775. The 3-fluorophenyl group extends into the solvent but is still within the protein cavity (Fig. 1d). Birb-796 forms multiple hydrophobic and polar interactions with the RET protein. The morpholino oxygen accepts an H-bond from the Ala807 backbone NH in the hinge region. The naphthalene moiety is positioned in the so-called Bback pocket^ sandwiched between the side chains of Lys758, Leu795, and Val804. The adjacent urea forms three H-bonds
with the RET protein: the oxygen interacts with the backbone NH of the Asp892 side chain (DFG) and both nitrogens interact with the side chain of the highly conserved residue Glu775. The Birb-796 pyrazole-3-t-butyl moiety is positioned between the side chains of Leu865, Leu870, and His872 (Fig. 1e). Motesanib forms one hydrogen bond to the RET hinge region through the pyridine N and the backbone amino nitrogen of Ala807. The amide nitrogen and oxygen of the inhibitor forms hydrogen bonds to the δ-carboxylate oxygen of Glu775 and the amino nitrogen of Asp892, respectively. The dihydroindole group extends into the cavity formed by Leu865, Leu870 and His872 (Fig. 1f). Sorafenib was reported to bind to p38α in the DFG-out form with two different ligand conformations. We thus adopted both binding conformations to RET also in this study. In sorafenib binding mode I (Sorafenib-I), two hydrogen bonds were formed in the hinge region, one via the pyridine N to the backbone amino nitrogen of Ala807, and the other via the ethyl picolinamide nitrogen to the backbone amino oxygen of Ala807. At the opposite end of the inhibitor, the lipophilic trifluoromethyl phenyl ring inserts into a hydrophobic pocket formed between the DFG-motif and the catalytic loop by the side chains of residues Leu865, Leu870, His872 and Ile890. Importantly, the urea group of the inhibitor forms two hydrogen bonds with the protein, one via its amide nitrogen atom to the carboxylate side chain of Glu755, and the second via the carbonyl moiety to the main chain nitrogen of Asp892 of the DFG-motif (Fig. 1g). In sorafenib binding mode II (Sorafenib-II), sorafenib loses the hydrogen bonding to the hinge region due to a rotation of the pyrimidine ring, whereas the remaining part of the inhibitor displays the same interaction pattern as in sorafenib binding mode I (Fig. 1h).
MD simulations Overall structural flexibility and stability MD simulations of the five complexes in explicit aqueous solution were run for a duration of 100 ns. To explore the stability of the complexes and to ensure the correctness of the sampling method, root-mean-square deviations (RMSD) from the starting structures were analyzed. The RMSDs show that the studied systems reach equilibrium within 20–30 ns, and indicate that the systems are stable (Fig. 2a). Analyses of root-mean-square fluctuations (RMSF) versus the residue number for the complex systems show that the protein Fig. 2 a Root-mean-square deviations (RMSDs) of RET Cα atoms for all five complexes. b Root-mean-square fluctuations (RMSFs) of each RET residue for all five complexes. c Distance fluctuations of protein and inhibitor atoms forming hydrogen bonds during the 100 ns molecular dynamics (MD) simulations
J Mol Model (2015) 21:167
A
C
Page 7 of 12 167
B
167
J Mol Model (2015) 21:167
Page 8 of 12
structures share similar RMSF distributions and similar trends in their dynamic features. The active site region (Lys728, Leu730, Val738, Ala756, Lys758, Val778, Leu779, Ile788, Leu802, Val804, Glu805, Tyr806, Ala807, Lys808, Gly810, Val871, His872, Leu881, Ser891, Asp892, Phe893) is relatively rigid, and this suggests all inhibitors have similar interaction mechanism with RET (Fig. 2b). The main motions occur, as expected, in the loop regions. The important hydrogen bond interactions of the RET inhibitors, as outlined above, were analyzed by monitoring fluctuations in donor–acceptor distances. A hydrogen bond was defined by having a heavy atom distance
ORIGINAL PAPER
Characterization of interactions and pharmacophore development for DFG-out inhibitors to RET tyrosine kinase Chunxia Gao 1 & Morten Grøtli 1 & Leif A. Eriksson 1
Received: 4 February 2015 / Accepted: 19 May 2015 # Springer-Verlag Berlin Heidelberg 2015
Abstract RET (rearranged during transfection) tyrosine kinase is a promising target for several human cancers. Abt348, Birb-796, Motesanib and Sorafenib are DFG-out multikinase inhibitors that have been reported to inhibit RET activity with good IC50 values. Although the DFG-out conformation has attracted great interest in the design of type II inhibitors, the structural requirements for binding to the RET DFG-out conformation remains unclear. Herein, the DFG-out conformation of RET was determined by homology modelling, the four inhibitors were docked, and the binding modes investigated by molecular dynamics simulation. Binding free energies were calculated using the molecular mechanics/Poisson-Bolzmann surface area (MM/PBSA) method. The trends in predicted binding free affinities correlated well with experimental data and were used to explain the activity difference of the studied inhibitors. Per-residue energy decomposition analyses provided further information on specific interaction properties. Finally, we also conducted a detailed e-pharmacophore modelling of the different RET-inhibitor complexes, explaining the common and specific pharmacophore features of the different complexes. The results reported herein will be useful in future rational design of novel DFG-out RET inhibitors. Keywords RET . DFG-out inhibitors . Molecular dynamics simulation . MM-PB(GB)SA . e-pharmacophore Electronic supplementary material The online version of this article (doi:10.1007/s00894-015-2708-z) contains supplementary material, which is available to authorized users. * Leif A. Eriksson [email protected] 1
Department of Chemistry and Molecular Biology, University of Gothenburg, 405 30 Göteborg, Sweden
Introduction RET (rearranged during transfection) is a transmembrane tyrosine kinase receptor essential for development of the sympathetic, parasympathetic and enteric nervous systems, and of the kidneys [1]. Furthermore, it is well established that deregulation and dysfunction of RET can lead to several different thyroid cancers, including multiple endocrine neoplasia 2 syndromes (MEN2A and MEN2B), familial medullary thyroid carcinoma (FMTC), and papillary thyroid carcinoma (PTC) [1]. Point mutations in the RET gene are the cause of the three cancer syndromes MEN2A, MEN2B, and FMTC, while PTC is related to chromosomal inversions or translocations. Furthermore, recent studies have demonstrated that RET is overexpressed in a subset of estrogen-receptor (ER)positive breast cancers and that crosstalk between ER and RET is important in responses to endocrine therapy [2, 3]. The tyrosine kinase domain of RET is comprised of two lobes: a smaller N-terminal lobe consisting mainly of β-strands and a single large α-helix; and a larger C-terminal lobe, which is almost exclusively α-helical. The ATP-binding site is located at the hinge region between the two lobes. The structure of the RET tyrosine kinase domain can be further divided into a number of structural regions; these include the glycine-rich loop (G-loop), the C α-helix on the Nlobe, the activation loop (A-loop), the catalytic loop (C-loop) and conserved DFG motifs on the C-lobe (Fig. 1a). In all known protein kinase structures, the DFG motif assumes a conformation with the Phe residue buried in a hydrophobic pocket in the groove between the lobes (DFG-in conformation). Recently, much attention has been devoted to the use of small molecule inhibitors to target RET in MEN2-related tumors; e.g., pyrazolopyrimidines PP1 and PP2, indolocarbazoles, 2indolinone derivatives, and the 4-anilinoquinazoline ZD6474 have been shown to display strong activity towards RET kinase
167
J Mol Model (2015) 21:167
Page 2 of 12
Fig. 1 a Ribbon representations of DFG-out RET (rearranged during transfection) tyrosine kinase structures. Key secondary structural elements are colored green (glycine-rich loop), red (Cα-helix), blue (activation loop) and pink (catalytic loop). b Alignment of DFG-in and DFGout residues. c Key residues of RET interacting with the inhibitors. d–h Illustration of RETDFG out inhibitor interaction: d Abt348, e Birb-796, f Motesanib, g Sorafenib-I, h Sorafenib-II
A
C
B
D
E
F
G
H
[4–6]. These inhibitors are all ATP competitive and bind exclusively within the ATP binding site of the active DFG-in kinase conformation; so-called type I inhibitors. Since the kinase catalytic domain is highly conserved and the ATP binding sites are similar, type I inhibitors are not specific to RET alone, but exhibit multi-kinase targeting. Instead of utilizing the allosteric pocket, which is much less conserved, a second class of inhibitors has been developed that act by stabilizing the inactive DFG-out kinase conformation [7]. This class, known as type II inhibitors, appear to display much higher selectivity. The DFG-out conformation has thus attracted significant interest in the discovery and development of selective kinase inhibitors in the past few years [8–14]. However, a recent analysis of binding modes of type II
inhibitors indicates that there may be less selectivity than initially believed [15]. It is thus important to investigate the interactions and functions of this type of inhibitor to better understand the key interactions and how selectivity may be enhanced further. Although not developed as RET inhibitors per se, several type II inhibitors have proven to be highly potent RET inhibitors, for example Sorafenib (IC50 =50 nM) [16] and Motesanib (IC50 =59 nM) [17] have both have reached clinical Phase II trials for treatment of MTC. Besides targeting RET, the latter drugs also mainly target vascular endothelial growth factor receptors (VEGFRs), platelet-derived growth factor receptor (PDGFR), and the KIT receptor protein-tyrosine kinase [18]. Two more compounds, Birb-796, a highly potent p38α inhibitor for anti-inflammatory therapy, and Abt-348, a novel
J Mol Model (2015) 21:167
Page 3 of 12 167
aurora kinase inhibitor, have also been tested and show IC50 values towards RET of 30 nM [19] and 7 nM [20], respectively (Table 1). To summarize, the four compounds cited above have been shown to have high experimental IC50 values to RET, and a study of how those compounds interact with RET is crucial for the development of more potent and selective RET inhibitors. Currently, X-ray crystal structures are available only for compounds binding to targets other than RET; i.e., Sorafenib binding to VEGFR2 (PDB code: 4ASD) [21], Motesanib binding to VEGFR2 (PDB code: 3EFL), and Birb-796 binding to p38α (PDB code: 1KV2) [13]. As there are no NMR or X-ray crystal structures of these compounds binding to RET, detailed information about the RET-type II inhibitor interaction is lacking. Such information is urgently needed in order to further improve our understanding of the structural determinants of the type II kinase inhibitor in binding and stabilizing the DFG-out conformation of RET. This will assist towards the rational design of novel, potent and selective RET kinase inhibitors. In this study, the DFG-out conformation of RET was obtained by homology modelling using VEGFR2 as template. After docking of the type-II inhibitors to the DFG-out conformation of RET, molecular dynamics (MD) simulations combined with free energy calculations were carried out to obtain
Table 1
information about binding of the inhibitors. Absolute binding free energy calculations were performed using the molecular mechanics Poisson-Bolzmann surface area (MM-PBSA) method [22]. Detailed binding free energies between the inhibitors and RET residues were calculated using a per-residue molecular mechanics/generalized Born surface area (MM-GBSA) decomposition analysis [23]. Furthermore, based on the optimized MD complex structures of RET and its inhibitors, multiple energy-based pharmacophore (e-pharmacopohore) models were constructed. We thus investigated features of the interaction between RET and type II inhibitors based on dynamic structural information from MD simulations, and outlined their specificity by analyzing the contribution of such interactions to the binding free energy of RET complexes, and their key pharmacophore descriptors.
Methods Homology model of DFG-out RET The tyrosine kinase domain (from Phe 719 to Lys 1011) of RET in its DFG-out conformation was obtained by homology
Structures and IC50 values of different inhibitors
NAME
IC50 (nM)
SORAFENIB
50 (Ref 16)
MOTESANIB
59 (Ref 17)
ABT-348
7 (Ref 20)
BIRB-796
30 (Ref 19)
STRUCTURE
167
J Mol Model (2015) 21:167
Page 4 of 12
modeling. The protein sequence of RET was obtained from the NCBI database (GeneID: 547807). The crystal structure of VEGFR2 (PDB code 4ASD, 2.03 Å resolution) [21] was chosen as the template for homology modeling calculations based on the fact that the tyrosine domain of VEGFR and RET share 47 % similarity, and that the DFG-out loop is present in its entirety in this crystal structure. A previous study dealing with construction and validation of a RET tyrosine kinase catalytic domain through homology modelling claimed that the Bopen^ and Bclosed^ states were both modelled [24]. However, both states reported in the study have the same DFG-in conformation; hence homology modelling of the RET DFT-out conformation in our study is not analogous to previous work. Intermediate homology models (25 backbone models×50 side chain models; in total 1250) were generated, and the bestscoring model as determined by the Generalized Born/ Volume Integral (GB/VI) scoring method [25], which compares the electrostatic solvation energy, was taken to be the final model. For the homology modelling, the MOE 2013.08 programme was used [26]. A Ramachandran plot of the final model (Supplementary material, Figure S1) shows that 98.4 % of the residues lie in the allowed regions, and that the 1.6 % outliers are in areas that are not in contact with the active site or affected by ligand binding. Docking The homology model of RET was superposed with the crystal structure of VEGFR2 in complex with Sorafenib (PDB code 4ASD) [21], after which the VEGFR2 structure was deleted and Sorafenib retained in the binding pocket of RET. The resulting RET–Sorafenib complex was prepared and refined using the default protein preparation wizard in Schrodinger. After ensuring chemical correctness, hydrogens were added, and side chains far from the binding cavity and not participating in salt bridges were neutralized. The hydrogen bonding network was optimized by reorienting the hydroxyl groups, water molecules, and amide groups of Asn and Gln, and selecting appropriate states and orientations of the imidazole ring in His residues. Water molecules present in the crystal structure were deleted, and termini were capped by adding ACE and NMA residues. Missing side chains and loops were added. Subsequently, the structures were minimized using the OPLS-2005 force field [27] and the Impact molecular mechanics engine, while heavy atoms were constrained. The prepared protein structure was then used for subsequent grid generation. Grids were generated by the Receptor Grid Generation panel (Schrodinger, LLC, New York, NY) which defines the receptor structure by excluding any co-crystallized ligand present, determines the position and size of the active site to be represented by receptor grids, and sets up Glide constraints. Grids were defined by centering them on the ligand using the
default box size, and with H-bonds to Glu775 and Asp892 set to be constrained. Molecular docking experiments with Abt-348, Birb-796, Motasenib and Sorafenib were carried out using Glide, as implemented in the Schrodinger package. The XP (extra precision) [28] docking function of Glide was used, granting full flexibility to the ligands. A post-docking minimization was performed on the output complexes in order to reduce the initially collected 10,000 poses ligand to five. The docked poses were then minimized using the local optimization feature in Prime. For Sorafenib, two different binding modes was observed in the different crystal structures (PDB code: 3HEG, 3GCS) [29, 30]; hence, two different docking poses were selected for the following studies. To validate the docking approach, Birb-796, Motasenib and Sorafenib (two structures) were docked back into the crystal structures of their respective receptors (Supplementary material, Figure S2). As can be seen from these benchmark dockings, the molecules overlap almost perfectly to their crystal structures, with very low RMSD values. MD simulations Atomic charges on Abt-348, Birb-796, Motesanib and Sorafenib were calculated using the restrained electrostatic potential (RESP) procedure at the HF/6-31G* level, after minimizing the molecule at the AM1 semiempirical level [31]. GAFF force field parameters [32] and the RESP partial charges were assigned using the Antechamber module in the AMBER10 package [33]. All missing hydrogen atoms were added using the LEaP module. Five conventional 100 ns MD simulations were performed on the RET-inhibitor complex structures using the Amber 99 force field [34, 35] and the GROMACS software [36]. The structures were solvated in periodic boxes with a buffer distance of 10.0 Å. A number of Na+ and Cl− ions were added to satisfy the electroneutrality condition and to give a salt concentration of 0.1 mol L−1, using the genion module in GROMACS. The obtained systems were energy minimized by steepest descent (200 steps) to remove close contacts. Position restrained simulations (2 ns duration, 1.0 fs time step, NPT ensemble, T=298 K, P=1 bar) were first performed, to enable the water molecules to reach more favorable positions. Particle-mesh Ewald (PME) [37, 38] summation was used for long-range electrostatics. A 12 Å cutoff was used for both Coulomb and Lennard-Jones interactions. The temperature and pressure was controlled through the Berendsen coupling algorithm [39], with the time constants 0.1 ps for temperature and 1.0 ps for pressure coupling. All bond lengths were constrained using the LINCS algorithm [40]. During the production simulations (200 ns duration, 1.0 fs time step, NPT ensemble, T=298 K, P=1 bar), the temperature was controlled using the Nose-Hoover thermostat [41] with a time
J Mol Model (2015) 21:167
Page 5 of 12 167
constant 0.1 ps, and the pressure was controlled using the Parrinello-Rahman barostat [42] with a time constant 1.0 ps. The remaining parameters were the same as in the position restrained simulations.
10−4 kcal mol−1 Å−2. Due to the high computational demand of this approach, only 50 snapshots that were evenly extracted from the last 2 ns MD trajectory were used to calculate the entropic contribution.
MM-PBSA calculations
Residue-inhibitor interaction decomposition
The binding free energies were calculated using the MMPBSA method implemented in AMBER10 [33]. For each complex, a total number of 500 snapshots were taken from the last 2 ns on the MD trajectory, and the free energy computed for each molecular species (complex, protein, and ligand), giving the binding free energy as the difference ΔΔGbind ¼ ΔGcomplex – ΔGprotein þ ΔGligand ð1Þ
To identify the key residues responsible for the binding of the inhibitors, free-energy decomposition was performed for each residue. Given the pairwise nature of the GB equation allowing the decomposition of polar solvation energy (ΔGele,sol) into atomic contributions in a straightforward manner, the interaction between the inhibitor and each residue was computed using the MM-GBSA decomposition process in AMBER10 [3, 23]. The binding interaction of each inhibitor–residue pair includes four terms: van der Waals contribution (ΔGvdW), electrostatic contribution (ΔGele), polar solvation contribution (ΔGele,sol), and nonpolar solvation contribution (ΔGnonpol,sol):
Each term can be estimated as follows: ΔG ¼ ΔGMM þ ΔGsol −TΔS
ð2Þ
where ΔGMM is the molecular mechanics free energy, ΔGsol is the solvation free energy, and TΔS represents the entropy term. The molecular mechanics energy was calculated from the electrostatic and van der Waals interactions, ΔGMM ¼ ΔGele þ ΔGvdW
ð3Þ
and the solvation free energy from the polar and nonpolar contributions, ΔGsol ¼ ΔGele;sol þ ΔGnonpol;sol
ð4Þ
The polar part ΔGele,sol was obtained by solving the Poisson-Boltzmann (PB) equations [43] for the MM-PBSA method, or the generalized Born (GB) model [43] for the MM-GBSA calculations. Dielectric constants of 1.0 and 80.0 were used for solute and solvent, respectively. The parameters used for the GB calculation were as developed by Onufriev et al. [44] (GBOBC, igb=2). The nonpolar solvation contribution was determined using the solvent accessible surface area (SASA, Å2) with the Amber molsurf module, according to ΔGnonpol;sol ¼ γSASA þ b
ð5Þ
where γ, representing the surface tension, and the constant b were set to 0.0072 kcal mol−1 Å−2 and 0, respectively. The conformational entropy upon ligand binding (translation, rotation, and vibration), TΔS, was calculated using normalmode analysis with the NMODE program in AMBER10 [32]. Prior to the normal-mode calculations, each snapshot was minimized for a maximum number of 100,000 cycles using the conjugate gradient method with in the presence of a distance-dependent dielectric of 4rij (rij being the distance between two atoms I and j) until the root-mean-square value of the elements of the gradient vector were less than 1.0×
ΔGinhibitor‐residue ¼ ΔGele þ ΔGvdW þ ΔGele;sol þ ΔGnonpol;sol
ð6Þ
The vdW and electrostatic interactions between inhibitor and each residue of RET were computed using the Sander program in AMBER10 [33]. The polar solvation contribution (ΔGele,sol) was calculated using the GB model [43] and the nonpolar contribution of desolvation (ΔGnonpol,sol) was computed based on SASA. All energy components were calculated using 500 snapshots extracted from the last 2 ns of each MD trajectory. e-pharmacophore hypothesis generation The structure-based pharmacophore was generated using the energy-minimized final snapshot of each simulated RETinhibitor complex. The Receptor Grid Generation tool in the Maestro software package was used to generate energy grids for all prepared protein structures and the ‘refine only’ option was used in the XP [27] docking in Glide 4.2 X. Default settings were employed to minimize and optimize the structures. On the basis of the XP descriptor information, pharmacophore features were generated using PHASE 3.4. Pharmacophore sites were generated automatically with Phase v3.0 (Schrodinger) using the default set of six chemical features: hydrogen-bond acceptor (A), hydrogen-bond donor (D), hydrophobe (H), negative ionizable (N), positive ionizable (P), and aromatic ring (R)). Hydrogen-bond acceptor sites were represented as vectors along the hydrogen-bond axis in accordance with the hybridization of the acceptor atom. Hydrogen-bond donors were represented as projected points, located at the corresponding hydrogen-bond acceptor positions in the binding site. Projected points allow the possibility
167
J Mol Model (2015) 21:167
Page 6 of 12
for structurally dissimilar active compounds to form hydrogen bonds to the same location, regardless of their point of origin and directionality. Each pharmacophore feature site was first assigned an energetic value equal to the sum of the Glide XP contributions of the atoms comprising the site, allowing sites to be quantified and ranked on the basis of the energetic terms The Glide XP descriptors include terms for hydrophobic enclosure, hydrophobically packed correlated hydrogen bonds, electrostatic rewards, π–π stacking, π–cation, and other interactions. ChemScore, hydrogen bonding, and lipophilic atom pair interaction terms were also included when the Glide XP terms for hydrogen binding and hydrophobic enclosure were zero.
Results and discussion Static structural analysis of inhibitor binding mode A large conformational change in the residues in the conserved DFG motif in the active site of the kinase is required for the binding of the DFG-out inhibitors. In the RET DFGout inhibitor complexes, the Phe side chain has moved by ~10 Å to a new position (DFG-out conformation). This movement of the Phe side chain reveals a large hydrophobic pocket in the kinase (Fig. 1b). The complexes of the four inhibitors with RET reveal that the inhibitors share the same binding mode, with RET kinase in the DFG-out conformation. The inhibitors are stabilized in a cavity formed by residues Lys728, Leu730, Val738, Ala756, Lys758, Glu775, Val778, Leu779, Ile788, Leu802, Val804, Glu805, Tyr806, Ala807, Lys808, Tyr809, Gly810, Val871, His872, Leu881, Ser891, Asp892, Phe893 (Fig. 1c). Abt-348 has the highest molecular weight among the inhibitors and forms the most interactions with RET. The pyrazol group is positioned such that the hydroxyethyl group forms two hydrogen bonds, one with the Tyr806 side chain and the other one with Lys808 backbone oxygen. The adjacent thienol group interacts with the hinge region via two hydrogen bonds, between the pyrimidine N and the backbone amino nitrogen of Ala807, and between the amino NH2-group and the backbone amino oxygen of Glu805. The urea fragment forms three H-bonds with the RET protein: the oxygen interacts with the backbone NH of the Asp892 side chain (in DFG) and both nitrogens interact with the side chain of the highly conserved residue Glu775. The 3-fluorophenyl group extends into the solvent but is still within the protein cavity (Fig. 1d). Birb-796 forms multiple hydrophobic and polar interactions with the RET protein. The morpholino oxygen accepts an H-bond from the Ala807 backbone NH in the hinge region. The naphthalene moiety is positioned in the so-called Bback pocket^ sandwiched between the side chains of Lys758, Leu795, and Val804. The adjacent urea forms three H-bonds
with the RET protein: the oxygen interacts with the backbone NH of the Asp892 side chain (DFG) and both nitrogens interact with the side chain of the highly conserved residue Glu775. The Birb-796 pyrazole-3-t-butyl moiety is positioned between the side chains of Leu865, Leu870, and His872 (Fig. 1e). Motesanib forms one hydrogen bond to the RET hinge region through the pyridine N and the backbone amino nitrogen of Ala807. The amide nitrogen and oxygen of the inhibitor forms hydrogen bonds to the δ-carboxylate oxygen of Glu775 and the amino nitrogen of Asp892, respectively. The dihydroindole group extends into the cavity formed by Leu865, Leu870 and His872 (Fig. 1f). Sorafenib was reported to bind to p38α in the DFG-out form with two different ligand conformations. We thus adopted both binding conformations to RET also in this study. In sorafenib binding mode I (Sorafenib-I), two hydrogen bonds were formed in the hinge region, one via the pyridine N to the backbone amino nitrogen of Ala807, and the other via the ethyl picolinamide nitrogen to the backbone amino oxygen of Ala807. At the opposite end of the inhibitor, the lipophilic trifluoromethyl phenyl ring inserts into a hydrophobic pocket formed between the DFG-motif and the catalytic loop by the side chains of residues Leu865, Leu870, His872 and Ile890. Importantly, the urea group of the inhibitor forms two hydrogen bonds with the protein, one via its amide nitrogen atom to the carboxylate side chain of Glu755, and the second via the carbonyl moiety to the main chain nitrogen of Asp892 of the DFG-motif (Fig. 1g). In sorafenib binding mode II (Sorafenib-II), sorafenib loses the hydrogen bonding to the hinge region due to a rotation of the pyrimidine ring, whereas the remaining part of the inhibitor displays the same interaction pattern as in sorafenib binding mode I (Fig. 1h).
MD simulations Overall structural flexibility and stability MD simulations of the five complexes in explicit aqueous solution were run for a duration of 100 ns. To explore the stability of the complexes and to ensure the correctness of the sampling method, root-mean-square deviations (RMSD) from the starting structures were analyzed. The RMSDs show that the studied systems reach equilibrium within 20–30 ns, and indicate that the systems are stable (Fig. 2a). Analyses of root-mean-square fluctuations (RMSF) versus the residue number for the complex systems show that the protein Fig. 2 a Root-mean-square deviations (RMSDs) of RET Cα atoms for all five complexes. b Root-mean-square fluctuations (RMSFs) of each RET residue for all five complexes. c Distance fluctuations of protein and inhibitor atoms forming hydrogen bonds during the 100 ns molecular dynamics (MD) simulations
J Mol Model (2015) 21:167
A
C
Page 7 of 12 167
B
167
J Mol Model (2015) 21:167
Page 8 of 12
structures share similar RMSF distributions and similar trends in their dynamic features. The active site region (Lys728, Leu730, Val738, Ala756, Lys758, Val778, Leu779, Ile788, Leu802, Val804, Glu805, Tyr806, Ala807, Lys808, Gly810, Val871, His872, Leu881, Ser891, Asp892, Phe893) is relatively rigid, and this suggests all inhibitors have similar interaction mechanism with RET (Fig. 2b). The main motions occur, as expected, in the loop regions. The important hydrogen bond interactions of the RET inhibitors, as outlined above, were analyzed by monitoring fluctuations in donor–acceptor distances. A hydrogen bond was defined by having a heavy atom distance
