FEP+ calculations predict a stereochemical SAR switch for first-in-class indoline NIK inhibitors for multiple myeloma

We recently reported the discovery and biological characterization of JNJ64290694 (Figure 1) as a potent NIK kinase inhibitor with strong anti-tumor activity in aggressive lymphoid malignancies, including multiple myeloma (MM), that are driven by the noncanonical NF-κB pathway [1]. Additional reports previously documented the optimization of NIK inhibitors for immune and auto-immune diseases [2–4].

JNJ64290694 potently inhibits phospho-IKKα, prevents nuclear accumulation of p52/RelB and represses the associated NF-κB gene program selectively in MM cell lines with genetic activation of the noncanonical NF-κB pathway. Proliferation of NIK translocated, TRAF3 or BIRC3 mutant MM cell lines is preferentially inhibited by JNJ64290694 over MM cell lines which lack genetic activation of the noncanonical NF-κB pathway. A single, oral dose of JNJ64290694 to mice bearing an NIK-translocated MM tumor inhibits phospho-IKKα and represses P52-mediated transcription of NF-κB-regulated genes. Daily oral dosing of JNJ64290694 completely inhibits expansion of NIK or TRAF3 mutant multiple myeloma tumors, with no signs of toxicities in these mouse models.

The type I kinase inhibitor cyanoindoline JNJ64290694 was generated by medicinal chemistry optimization of an azaindoline hit obtained from a biochemical NIK autophosphorylation high-throughput screening assay [5,6].The binding potency SAR developed in the hit-to-lead and lead optimization phases explored different regions of the molecular system (see Figure 1 for chemical structures of compounds and naming codes included in this study), including: the R/S stereoswitch (i.e., the R enantiomers have in general better potency compared with the corresponding S enantiomers), the azaindoline/cyanoindoline switch (i.e., cyanoindolines have in general better potency compared with the corresponding azaindolines), the adenine pocket filling pyrimidine versus 5-fluoro/chloropyrimidine versus triazine cores, and the front pocket filling aminophenyl/aminopyridyl substitutions (see Supplementary Figure 1 for illustration of the SAR regions). An SAR switch was herein defined as a structural modification which was repeatedly observed to drive the potency in one same direction. Compound design was supported during hit-to-lead and lead optimization phases in a prospective and predictive manner by in silico docking, FEP+ and WaterMap calculations [7–9]. The calculations were used to explore SAR questions or to prioritize chemistry plans at the moment when they were established. More specifically the R/S stereoswitch calculations were predictive. Initially the team was working with racemic compounds and requested computational support to address the preference. Herein we report on the computational support provided by FEP+ and WaterMap calculations and show how this type of computational chemistry methods in combination with structural analysis were useful during hit-to-lead and lead optimization phases to gain better understanding of the different SAR regions and especially the R/S stereoswitch.

Computational methods

All calculations were run using the Schrödinger (NY, USA) software suite 2017-2 [10].

Glide SP docking

Glide docking approximates a complete systematic search of the conformational, orientational and positional space of the docked ligand in a rigid protein receptor binding pocket [7,11]. An initial rough positioning and scoring phase that dramatically narrows the search space is followed by torsionally flexible energy optimization on an OPLS3-based nonbonded potential grid for a few hundred surviving candidate poses [12]. The very best candidates are further refined via a Monte Carlo sampling. Selection of the best docked pose uses the model energy function GlideScore that combines empirical and force-field-based terms.

The higher resolution protein structure of the S-isomer cyanoindoline compound 4S/3694 (2.31 Å resolution PDB: 6Z1T) was used for docking [13]. The protein and ligands were respectively prepared with the Protein Preparation Wizard and the LigPrep software using default settings. The ligands were modeled in the tautomeric and neutral form as depicted in Figure 1. Water molecules within the 5 Å contact sphere of the ligand were kept within the structure. Glide SP (Standard Precision) docking was used with positional core-constraints for the pyrimidine-azaindoline or pyrimidine-cyanoindoline scaffolds (0.5 Å tolerance) as observed in the crystal structures.

FEP+ simulations

FEP+ provides a computational workflow for protein–ligand binding-affinity calculations using Schrödinger's molecular dynamics/replica exchange capabilities in the Desmond MD software with the FEP/REST (free energy perturbation/replica exchange with solute tempering) algorithm which is a single molecular topology method (Supplementary Figure 2) [8,14,15].

The herein carried out simulations used the OPLS3 force field for the protein and the ligands along with the SPC water model [12]. The torsional parameters for ligand torsions that were not included in the default OPLS3 force field were generated using the Force Field builder tool. The Force Field builder tool fitted the parameters to quantum chemical data based on the construction of model fragments containing the relevant torsional structures. The perturbation pathways were either manually edited or automatically generated using the Mapper tool. In the automated set-up, the maximum common substructure between any pair of compounds were generated and their similarity was measured. Then ligand pairs with high similarity scores were connected by edges, where each edge represents one FEP calculation that will be performed between the two ligands. FEP+ calculations were run using starting poses generated as described above with default settings Glide SP docking based on the crystal structure of the S-isomer cyanoindoline compound 4S/3694 [13]. Crystallographic water molecules within 5 Å of the ligand were kept in the system. FEP+ calculations were done using the ‘-solvate’ keyword to enable Grand Canonical Monte Carlo positioning of additional water molecules into new cavities formed upon alchemical transitions. The systems were solvated in a water box with buffer width of 5 Å for the complex simulations and 10 Å for the solvent simulations. The systems were relaxed and equilibrated using the default relaxation protocol. A total of 16 λ windows were used for each FEP/REST calculation. The production λ stages lasted each 5 or 25 ns for both the complex and the solvent simulations using NPT ensemble conditions. Replica exchanges between neighboring λ windows were attempted every 1.2 ps. The Bennett acceptance ratio method was used to calculate the free energy [16]. Errors were estimated for each free energy calculation using both bootstrapping and the Bennett acceptance ratio analytical error prediction [17]. After the simulations were completed, the hysteresis along closed thermodynamic cycles was calculated and best estimates of the free energies and the errors for the predictions were calculated using the cycle closure algorithm.

WaterMap

Through a workflow of restrained molecular dynamics simulations, solvent clustering and statistical thermodynamic analysis, WaterMap computes various properties of water molecules, including location, occupancy, enthalpy, entropy and free energy [18,19].

Two ns explicit-solvent MD simulation of the apo-protein or a protein–ligand complex were performed using the program Desmond. The OPLS3 force field was used for protein and ligand atoms together with the TIP4P water model [12]. In the truncated protein set-up mode, only residues within 10 Å of the ligand were simulated as restrained flexible residues. The coordinates of the protein were restrained with a 5.0 kcal/mol/Ų harmonic potential applied to the initial positions of the heavy atoms, which ensured convergence of the water sampling around the protein conformation of interest.

The WaterMap truncation is a real truncation of the system. All residues with at least one heavy atom within the truncation distance from any heavy atom in the ligand are retained. Broken bonds are reformed and capped with hydrogen atoms. This has no effect on the dynamics since the restraints keep the shape of the protein despite missing most of it.

Water molecules from approximately 2000 equally spaced frames from the MD simulation were then spatially clustered to form localized hydration sites, and the thermodynamic properties of those sites with the bulk solvent as reference state were computed.

WaterMap simulations were done with both crystal structures Des-F 3R/4076 and 4S/3694 as input, both in the apo and the liganded forms. The protein structures were prepared with the Protein Preparation Wizard. Water molecules resolved in the crystal structures within 10 Å contact sphere of the ligands were kept included in the simulations and treated as solvent.

The excess free energy of individual hydration sites on a protein surface can be either positive (unfavorable) or negative (favorable) compared with a site in bulk solvent. A positive excess ΔG means a water molecule at this position that is unstable and displacement by a ligand will contribute favourably to ligand binding. A negative excess ΔG means a water that is stable and that it would cost energy to displace it from this position. The fact that the excess ΔG can be positive can lead to the question why unstable water do not simply leave the binding site and produce in the binding site a vacuum cavity [9]. This vacuum cavity formation is associated with a penalty that has been estimated to be around 6–8 kcal/mol. Thus, unstable water molecules will not leave a cavity unless their excess ΔG is more positive than the formation of the vacuum cavity [20,21].

Results & discussion

FEP+ calculations

At Janssen, we have been investigating the feasibility of FEP+ calculations in ongoing drug discovery programs, regularly finding accuracies to within 1 kcal/mol of experiment, in agreement with initial reports from Schrödinger [22–25]. In the NIK project FEP+ calculations were initially performed to address local SAR questions, including four to ten compounds in each FEP+ map. While such small FEP+ maps may suffer from the viewpoint of statistical significance, they can play an important role in decision making during hit-to-lead or lead optimization chemistry.

As shown in Supplementary Figure 3, initial results of the FEP+ calculations showed a correlation for the azaindoline to cyanoindoline switch and the R/S stereochemical switch with DiscoverRx Kinome Scan pK_D data [26] and pIC₅₀ data from the AlphaScreen NIK autophosphorylation assay (SText 1). The cyanoindoline complexes were found approximately 1.6 kcal/mol more stable than the corresponding azaindoline complexes. The R stereoisomer complexes were found approximately 1 kcal/mol more stable than the corresponding S stereoisomer complexes.

This latter correlation was also consistently observed in retrospective FEP+ calculations (Figure 2, Tables 1 & Supplementary Table 1) including seven enantiomeric pairs of compounds shown in Figure 1. The comparison to the Glide docking score, which provides a simplified description of the binding energy, showed that it was not reliable to predict K_D values and could not predict the preference of R versus S stereoisomers; the reversed order was actually predicted based on the Glide docking score. The FEP+ calculations were robust to changes in the simulation time for individual λ windows – Supplementary Table 2. The default 5 ns per λ step simulations were strongly correlated with longer 25 ns per λ step simulations (R² = 0.96). Furthermore, the simulations were robust for changes of the FEP mapper network which is the map of performed alchemical transitions – Supplememtary Figure 4 & Table 2. A manually edited map was compared with the automatically generated map using default settings and correlated very well (R² = 0.98). The manually edited map included the direct transitions between the 7 R/S enantiomer pairs which were not all included in the automated map. As a result of this, we observed that intrinsic prediction accuracy (as computed by the Bennett error) and the cycle closure error were slightly improved using the manual maps. The same observations applied comparing the 25 ns per λ step to the 5 ns per λ step simulations using the automated map.

No clear correlation was observed for the other SAR regions mentioned in the introduction, and also not for a comprehensive retrospective FEP+ map including 23 compounds addressing all SAR regions (R² = 0.08 correlation with DiscoverRx Kinome Scan pK_{D s} – data not shown). Noteworthy was also that cell-based assays (e.g., inhibition of proliferation of JJN3 cells R² = 0.34 – see Supplementary Text 2) did not correlate to the predicted free energy of binding; pointing to the fact that FEP calculations should be compared with direct experimental binding data and not to convoluted readouts.

Structural analysis of crystal structures

The analysis of crystal structure complexes helped to understand the origin of the differences in binding energies – Figure 3 [13]. The structures are typical type I DFG-in and C-α-Helix-in structures [27]. The structure of the R-isomer azaindoline compound Des-F(fluoro) 3R/4076 in complex with the NIK protein catalytic kinase domain, which was resolved up to 2.42 Å, showed that the 2-(phenylamino)-pyrimidinyl-azaindoline scaffold made polar contacts to protein backbone groups and conserved polar side chain residues. As expected, the 2-aminopyrimidyl group interacted with hinge residues Glu 470 and Leu 472. The azaindoline NH group interacted with the sidechain of the DFG-in motif Asp534 residue and the azaindoline acceptor N atom made a water (W2) bridged interaction to the ammonium group of the conserved catalytic Lys429. Most interestingly, the switch hydroxyl group was hydrogen bonded to the main chain carbonyl group of Asp519. In addition, the hydroxyl group made a hydrogen bond to one solvent accessible water molecule (W42) which was bridged to the side-chain carboxylate group of Asp519 and the side-chain OH group of Ser476. Three additional water molecules were resolved in the proximity of the ligand; W7 was hydrogen bonded to the main chain NH of Arg408. An additional interaction existed between the N-tetrahydropyranylamide NH group and the main chain carbonyl group of Leu406. In comparison to this, the corresponding structure of the S-isomer cyanoindoline compound 4S/3694 (2.31 Å) showed identical interactions with the kinase hinge motif. In both R and S isomer structures, the core scaffolds occupied identical positions and structure the 4S/3694 was used as input for the docking and FEP+ computations. However, the interactions of the cyanoindoline group were different. The cyano group was directly hydrogen bonded to the ammonium group of conserved catalytic Lys429. There was no space for the above-mentioned W2 molecule; it would clash with the cyano group. Similar to the 3R/4076 structure, the cyanoindoline NH group interacted with the sidechain of the DFG-in motif Asp534 residue. The switch hydroxyl group formed a different hydrogen bond, in this case to the main chain carbonyl group of Arg408. In addition, the hydroxyl group made a hydrogen bond to a solvent accessible water molecule (W17) which was bridged to the main chain NH of Arg408. The interactions of the switch hydroxyl group in the analyzed R and S stereoisomer structures were thus comparable in nature, for example – hydrogen bonds with protein atoms and water molecules. For the R isomer, the degree of structural organization appeared to be higher, with more structural water molecules being resolved in the proximity of the ligand. It was unclear whether this higher organization resulted in an enthalpic advantage or an entropic disadvantage. However, the net effect appeared to be that the R complexes were stabilized compared with corresponding S complexes.

WaterMap simulations

To explore the importance of the water network further, we carried out WaterMap simulations with both crystal structures as input, both in the apo and liganded forms – Supplementary Figure 5. The superimposition of the compounds Des-F 3R/4076 and 4S/3694 with the WaterMaps generated for their respective apo structures showed that upon ligand binding a number of high energy waters are excluded, which likely contributed to the binding affinity of the compounds. Given the high structural similarity of the two crystal structures we observed only minor differences for the two apo maps. The main differences were located at the cyanosubstitution point and could potentially be explained by the differences between the starting structures which did not converge in the 2 ns short restrained simulations. Significant differences due to the difference in substitution were observed for the liganded WaterMap simulations, especially for the front pocket region of the Des-F 3R/4076 compound. The water molecule bridging the azaindoline Des-F 3R/4076 and the catalytic Lys429 is a high energy water molecule (W2 +8.33 kcal/mol). In the case of the cyanoindoline 4S/3694, there was a high energy water molecule (W9 +6.68 kcal/mol) at a slightly displaced position and not bridging to the catalytic lysine. In agreement with this, we observed that in the two matched molecular pairs included in the dataset that involved the azaindoline to cyanoindoline switch, the cyanoindoline compounds had better potencies than the corresponding azaindoline compounds, both in prediction and in experiment (see Supplementary Figure 3). For low energy water molecules (ΔG <= -3.0 kcal/mol) within 5 Å of the ligand binding spheres, we observed that there are no major differences for the two liganded forms while there were differences for the apo forms with most notably the splitting in the Des-F 3R/4076 structure of one structural water cluster bond in both liganded forms between the sidechains of Asn520 and Asp534 (see W05 and W96 in Figure 3).

Analysis of the dynamic ligand–protein residue interactions

To investigate the potential importance of differences in specific ligand–protein contacts, we analyzed the dynamic ligand-protein residue interaction diagrams for end point λ replicates of the seven direct R/S transitions included in the FEP+ simulations using the manually defined Map – Supplementary Figure 6. As shown in Table 2, we observed that the interaction of the hydroxyl group in the R enantiomer with the Asp519 was in six out of seven direct R/S transitions more persistent over time in the simulation than the interaction of the hydroxyl group in S enantiomer with the Arg408 residue. This difference in persistence was a plausible reason for the binding energy differences of R versus S enantiomers. This information should ideally be augmented by the main average energy of the hydrogen bonds. It was interesting to note that in four out of the seven transitions this persistence difference of the switching hydroxyl group interaction was corresponding to a persistence difference of the interaction of the hinge Leu472 NH group with the central pyrimidine ring; indicating the presence of microsite-microsite correlation and cooperativity in the multi-microsite ligand–protein interaction. Such microsite-microsite correlations and cooperativity have their origins in the covariance of atomic fluctuations and are at the basis of non-additive SAR [28–31].

Conclusion

In conclusion, we reported the capability of FEP+ calculations to accurately predict a stereochemical SAR switch in the optimization of first-in-class indoline NIK inhibitors for multiple myeloma. The analysis of the dynamic interaction diagrams showed that for the R isomers, the H-bond interaction of the switching hydroxyl group with the main chain carbonyl group of Asp519 was more persistent than the H-bond interaction of the hydroxyl group with the main chain carbonyl group of Arg408 observed in the S isomers. The differences in structure and dynamics of the systems thus explained the observed differences of free binding energies. As a practical guideline, we learned from our study that the FEP+ maps should be improved or redesigned to include all transitions between matched molecular pairs corresponding to the SAR switches one aims to investigate.

In specific cases, FEP+ calculations can be helpful to address local potency related SAR questions which are of interest in hit-to-lead and lead optimization. Given that FEP+ is a holistic method, all contributing enthalpic and entropic factors are integrated. Further progress is expected to result from the investigation of methods to decompose the differences in binding free energy into ligand fragment, protein residue and solvent components [32–34]. This should enable the detection of potentially dominant interaction contributions in a more quantitative manner than the currently available qualitative analysis of interaction diagrams and be potentially useful to further direct the potency optimization of compounds.