Discovery of Novel Imidazopyridine GSK-3β Inhibitors Supported by Computational Approaches

The interest of research groups and pharmaceutical companies to discover novel GSK-3β inhibitors has increased over the years considering the involvement of this enzyme in many pathophysiological processes and diseases. Along this line, we recently reported on 1H-indazole-3-carboxamide (INDZ) derivatives 1–6, showing good GSK-3β inhibition activity. However, they suffered from generally poor central nervous system (CNS) permeability. Here, we describe the design, synthesis, and in vitro characterization of novel imidazo[1,5-a]pyridine-1-carboxamide (IMID 1) and imidazo[1,5-a]pyridine-3-carboxamide (IMID 2) compounds (7–18) to overcome such liability. In detail, structure-based approaches and fine-tuning of physicochemical properties guided the design of derivatives 7–18 resulting in ameliorated absorption, distribution, metabolism, and excretion (ADME) properties. A crystal structure of 16 in complex with GSK-3β enzyme (PDB entry 6Y9S) confirmed the in silico models. Despite the nanomolar inhibition activity, the new core compounds showed a reduction in potency with respect to INDZ derivatives 1–6. In this context, Molecular Dynamics (MD) and Quantum Mechanics (QM) based approaches along with NMR investigation helped to rationalize the observed structure activity relationship (SAR). With these findings, the key role of the acidic hydrogen of the central core for a tight interaction within the ATP pocket of the enzyme reflecting in good GSK-3β affinity was demonstrated.


Introduction
Glycogen synthase kinase-3 (GSK-3) is a ubiquitous, multifunctional and constitutively active serine/threonine kinase, involved in numerous pathways and over one hundred different putative substrates [1].
GSK-3β is a significantly networked enzyme, which regulates a plethora of physiological functions in peripheral tissues and the central nervous system (CNS), spanning from metabolism and cell cycle regulation to brain development. Accordingly, balanced GSK-3β regulation and signaling is key to human physiology but often compromised in diverse disease conditions such as diabetes, cancers, and brain disorders [6,7]. As for CNS-related pathologies, aberrant GSK-3β activity has been reported Despite few derivatives showing promising in vitro profiles, enhanced blood brain barrier (BBB) penetration was still desirable to hinder CNS disorders. Indeed, the BBB is a dynamic and functional structure which poses a vast challenge in the development of drugs acting on the CNS [20]. The major route of drug permeation through the BBB is passive diffusion, which can be boosted by fine-tuning several physicochemical properties. A variety of approaches has been put forward to guide the design of BBB penetrant compounds based upon these features [21][22][23].
On these bases, starting from the INDZ core, compound physicochemical property evaluation and computer-aided methods were applied to support a scaffold hopping approach leading to the identification of novel imidazo[1,5-a]pyridine-1-carboxamide (IMID 1) and imidazo [1,5-a]pyridine-3-carboxamide (IMID 2) derivatives 7-18 ( Figure 1 and Table 1) as putative GSK-3β inhibitors with improved in vitro CNS permeability. Specifically, structure-based approaches were applied to predict ligand-protein interactions, which were confirmed by the high-resolution co-crystal structure of compound 16 with human GSK-3β enzyme (PDB entry 6Y9S). Additionally, computational investigations based on Molecular Dynamics (MD) simulations and Quantum Mechanics (QM) helped to retrospectively rationalize the structure activity relationships (SARs) of the in vitro tested compounds. NMR-based studies were also carried out to further substantiate the observed experimental data. Despite few derivatives showing promising in vitro profiles, enhanced blood brain barrier (BBB) penetration was still desirable to hinder CNS disorders. Indeed, the BBB is a dynamic and functional structure which poses a vast challenge in the development of drugs acting on the CNS [20]. The major route of drug permeation through the BBB is passive diffusion, which can be boosted by fine-tuning several physicochemical properties. A variety of approaches has been put forward to guide the design of BBB penetrant compounds based upon these features [21][22][23].
On these bases, starting from the INDZ core, compound physicochemical property evaluation and computer-aided methods were applied to support a scaffold hopping approach leading to the identification of novel imidazo[1,5-a]pyridine-1-carboxamide (IMID 1) and imidazo [1,5-a]pyridine-3-carboxamide (IMID 2) derivatives 7-18 ( Figure 1 and Table 1) as putative GSK-3β inhibitors with improved in vitro CNS permeability. Specifically, structure-based approaches were applied to predict ligand-protein interactions, which were confirmed by the high-resolution co-crystal structure of compound 16 with human GSK-3β enzyme (PDB entry 6Y9S). Additionally, computational investigations based on Molecular Dynamics (MD) simulations and Quantum Mechanics (QM) helped to retrospectively rationalize the structure activity relationships (SARs) of the in vitro tested compounds. NMR-based studies were also carried out to further substantiate the observed experimental data.

Derivatives 1-6 and Rationale for New Design
The INDZ core was previously identified as an attractive scaffold for the development of novel GSK-3β inhibitors [16][17][18][19]. Through in silico and crystallographic studies [16,17,19], it was found to bind in the ATP binding domain of GSK-3β enzyme between the N-and C-terminal lobes (Figure 2a). The INDZ moiety is located at the adenine binding site and contacts the hinge residues via three hydrogen bonds (H-bonds) referred to as: 1) "deep" engaging the hydrogen (Ha) at position N-1 of the core with Asp133 carbonyl group, 2) "central" involving the N-2 of the core and Val135 N-H group, and 3) "outer" between the hydrogen (Hb) of the carboxamide group and Val135 carbonyl (the three interactions are represented in Figure 2a as red, green, and orange dashed lines, respectively). In addition, R1 and R2 substituents appended to the central core are oriented towards Arg141 in the external solvent accessible part of the binding pocket and the inner cavity in proximity of Lys85, respectively.

Derivatives 1-6 and Rationale for New Design
The INDZ core was previously identified as an attractive scaffold for the development of novel GSK-3β inhibitors [16][17][18][19]. Through in silico and crystallographic studies [16,17,19], it was found to bind in the ATP binding domain of GSK-3β enzyme between the N-and C-terminal lobes (Figure 2a). The INDZ moiety is located at the adenine binding site and contacts the hinge residues via three hydrogen bonds (H-bonds) referred to as: 1) "deep" engaging the hydrogen (Ha) at position N-1 of the core with Asp133 carbonyl group, 2) "central" involving the N-2 of the core and Val135 N-H group, and 3) "outer" between the hydrogen (Hb) of the carboxamide group and Val135 carbonyl (the three interactions are represented in Figure 2a as red, green, and orange dashed lines, respectively). In addition, R 1 and R 2 substituents appended to the central core are oriented towards Arg141 in the external solvent accessible part of the binding pocket and the inner cavity in proximity of Lys85, respectively.
Compounds 1-6, published ahead as potent GSK-3β inhibitors [17][18][19], are reported herein as representative derivatives of the INDZ series and served as starting point for an optimization exercise towards the development of novel compounds with improved CNS penetration. Compounds 1-6 showed outstanding GSK-3β inhibitory activity in the nanomolar range, with derivatives 2 and 4 being the most potent of the set (IC 50 = 4 and 9 nM, respectively). Interestingly, both compounds 2 and 4 carry an isopropoxypyridinyl group as R 2 substituent, resulting from 3 to 5-fold more potent than the corresponding pyridinyl analogues 1 and 3, respectively. The high-resolution X-ray crystal structure of 2 in complex with GSK-3β (2.08 Å, PDB entry 6Y9R) supported the expected binding mode (Figure 2b). Here, the piperidine chain was solvent exposed towards Arg141, the INDZ was anchored to the hinge domain via the tridentate H-bond interaction, the pyridine formed a H-bond with the catalytic Lys85, and the isopropoxyl group faced the ribose binding site. In proximity of the ligand, a network of hydration sites was found involving residues like Ile62, Leu132, Thr138, Gln185, Asn186, and Asp200. Interestingly, MD simulations (see below) highlighted that the isopropoxyl group of 2 might be also involved in additional water-mediated interactions, justifying the boost in potency with respect to the pyridinyl analogue 1. Replacement of the R 2 pyridinyl group of 3 with a (hydroxymethyl)pyridinyl (5) or a di-F-phenyl (6) also retained double-digit nanomolar inhibition activity. Finally, when comparing the enzymatic potency of the pyridinyl compound 1 against 3 or the isopropoxypyridinyl derivative 2 against 4, it appeared that the two matching pairs resulted almost equipotent, suggesting that the R 1 substituents do not significantly impact on the GSK-3β inhibitory activity. This is in line with the observed INDZ binding mode, where the R 1 group is solvent oriented.  [16,17,19]. The 2D structure of the INDZ core with explicit Ha and Hb atoms is also showed; (b) X-ray co-crystal structure of inhibitor 2 in complex with GSK-3β enzyme (PDB entry 6Y9R). The deep, central and outer H-bond interactions are represented as red, green, and orange dashed lines, respectively. The most relevant residues of the ATP binding site are showed in light blue (for clarity, light transparency is used for some amino acids). Relevant water molecules around inhibitor 2 are also displayed. Both figures were prepared with VMD 1.9.4 [24].
Compounds 1-6, published ahead as potent GSK-3β inhibitors [17][18][19], are reported herein as representative derivatives of the INDZ series and served as starting point for an optimization exercise towards the development of novel compounds with improved CNS penetration. Compounds 1-6 showed outstanding GSK-3β inhibitory activity in the nanomolar range, with derivatives 2 and 4 being the most potent of the set (IC50 = 4 and 9 nM, respectively). Interestingly, both compounds 2 and 4 carry an isopropoxypyridinyl group as R2 substituent, resulting from 3 to 5-fold more potent than the corresponding pyridinyl analogues 1 and 3, respectively. The high-resolution X-ray crystal structure of 2 in complex with GSK-3β (2.08 Å, PDB entry 6Y9R) supported the expected binding mode (Figure 2b). Here, the piperidine chain was solvent exposed towards Arg141, the INDZ was anchored to the hinge domain via the tridentate H-bond interaction, the pyridine formed a H-bond with the catalytic Lys85, and the isopropoxyl group faced the ribose binding site. In proximity of the ligand, a network of hydration sites was found involving residues like Ile62, Leu132, Thr138, Gln185, Asn186, and Asp200. Interestingly, MD simulations (see below) highlighted that the isopropoxyl group of 2 might be also involved in additional water-mediated interactions, justifying the boost in potency with respect to the pyridinyl analogue 1. Replacement of the R2 pyridinyl group of 3 with a (hydroxymethyl)pyridinyl (5) or a di-F-phenyl (6) also retained double-digit nanomolar inhibition activity. Finally, when comparing the enzymatic potency of the pyridinyl compound 1 against 3 or the isopropoxypyridinyl derivative 2 against 4, it appeared that the two matching pairs resulted almost equipotent, suggesting that the R1 substituents do not significantly impact on the GSK-3β inhibitory activity. This is in line with the observed INDZ binding mode, where the R1 group is solvent oriented.  [16,17,19]. The 2D structure of the INDZ core with explicit Ha and Hb atoms is also showed; (b) X-ray co-crystal structure of inhibitor 2 in complex with GSK-3β enzyme (PDB entry 6Y9R). The deep, central and outer H-bond interactions are represented as red, green, and orange dashed lines, respectively. The most relevant residues of the ATP binding site are showed in light blue (for clarity, light transparency is used for some amino acids). Relevant water molecules around inhibitor 2 are also displayed. Both figures were prepared with VMD 1.9.4 [24].
Alongside enzymatic potency evaluation, preliminary brain permeability needs also to be assessed for compounds to be successfully advanced in the CNS therapeutic area. Indeed, BBB penetration is still the bottleneck in brain drug development and is the single most important factor limiting the future growth of neurotherapeutics. The MDCK-MDR1 cell line is considered a good surrogate of the in vivo BBB model. Therefore, we report herein on compounds 2, 3, 4, and 6 that, being the most potent derivatives of the dataset, were selected as INDZ representatives for experimental testing in an in vitro MDCK-MDR1 permeability assay (Table 1) [25]. The permeability coefficients (P app ) were calculated by measuring the transport of the compounds from the apical (A) to the basolateral (B) direction (Table 1 and Equation (1)). Efflux Ratio (ER) was also measured (see Equation (2)) to quantify efflux transport across a cell monolayer expressing P-glycoproteins (P-gp), by assessing transport in both directions (A-B and B-A) [26].
Interestingly, compounds 3 and 6 showed moderate permeability (P app = 15.8 × 10 −6 and 13.3 × 10 −6 cm/s, respectively), whereas compounds 2 and 4 resulted only poorly permeable (P app < 1.5 × 10 −6 cm/s). Moreover, high ER values were obtained for derivatives 2 and 4 and lower ratios for compounds 3 and 6. The observed experimental data could be rationalized by analysis of the physicochemical profile. Indeed, several properties have consistently been found to be relevant for optimal brain exposure, such as lipophilicity, expressed by the calculated logarithm of the octanol/water partition coefficient (cLogP) or calculated logarithm of the octanol/water partition coefficient at physiological pH 7.4 (cLogD), topological polar surface area (TPSA), ionization state (pK a ), number of H-bond donors (HBD) and acceptors (HBA), rotatable bond (RB) count, and molecular weight (MW). Up to now, many published works endeavor to provide guidance in terms of physicochemical property range for the design of CNS penetrant molecules [27]. Many of these studies have focused on the profile of marketed CNS drugs. As an example, the median values derived from an analysis of marketed CNS drugs conducted at Ely Lilly define a CNS drug profile as: cLogP = 3.1, cLogD = 1.9, TPSA = 40.5 Å, most basic pK a = 7.9, HBD = 1, HBA = 4, RB count = 4, and MW = 300. Obviously, these numbers are median values, and there are many marketed drugs with properties at either end of the range [27]. Overall, CNS drug property studies, despite providing different values, came to the common conclusion that CNS drugs are generally smaller, more lipophilic, with fewer HBDs, and lower TPSA compared with oral non-CNS drugs [21]. At this regard, cLogD (pH = 7.4), TPSA, the most basic pK a , HBD, and MW of compounds 1-6 were calculated here with ACD/Percepta 2017.2, for permeability estimation. Overall, the higher P app values of derivative 3 and 6 with respect to 2 and 4 could be ascribed to a more desirable CNS profile, characterized by higher lipophilicity (in particular 6), reduced H-bond potential (TPSA), and lower molecular size (Table 1). Moreover, derivatives 2 and 4 showed higher ER compared to 3 and 6, a result in agreement with the observation that increasing the numbers of HBA and MW would confer greater likelihood of P-gp efflux [28,29].
In the context of CNS drug design, the HBD count, an important H-bond capacity descriptor, is a critical CNS differentiating parameter [21][22][23]. Indeed, it has been reported that increased HBD count is associated with lower passive solubility and higher risk of interaction with efflux transporters. Particularly, compounds with HBD > 2 are highly unlikely to be brain penetrant. Therefore, we conceived a new scaffold hopping design strategy to reduce the HBD count of the INDZ core (HBD = 2) and boost the BBB penetration of the current series. In particular, we replaced the indazole N-Ha with a C-Ha group and shuffled the ring nitrogen atoms resulting in the IMID 1 and IMID 2 cores (Table 1). These modifications were intended to reduce the acidity of the Ha atom so to enhance permeability, while retaining the tridentate H-bond interaction with the hinge residues and enzymatic potency.
We expected SAR to be reasonably transferrable from the INDZ to the IMID 1 and IMID 2 series; thus, new derivatives 7-18 were obtained by decorating the IMID 1 and IMID 2 cores with the same R 1 and R 2 substituents discussed above for INDZs 1-6.

In Silico Investigation of Novel 7-18 Derivatives
Taking advantage of the crystal structure of 2 in complex with GSK-3β, we ran Glide Standard Precision (SP) ligand docking [30,31] and Molecular Dynamics (MD) simulations [32] to evaluate the binding mode and the stability of the new series compared to the related INDZ compounds (see Materials and Methods for details). It followed that the docked ligands accommodated the central core in the ATP binding site similarly to the INDZ group, while major variability in the orientations of R 1 and R 2 groups was observed. Such conformational variability was mainly ascribed to: (i) the inherent flexibility of the R 1 appended moieties, (ii) the open solvent-exposed binding site occupied by the R 1 groups, and (iii) the wide binding site filled by the R 2 residues. Figure 3 highlights the good overlap of INDZ-IMID 1-IMID 2 triads bearing the same R 1 and R 2 groups.
Deep visual inspection and energetics considerations (from −8.95 to −7.72 kcal/mol docking score) led to the selection of one pose per ligand for further studies (Table 2). For a comprehensive understanding of the IMID 1 and IMID 2 binding mode, the RMSD of the triplets was first calculated based on the common substructure of the core scaffold and carboxamide portion (excluding hydrogen atoms), while the more flexible R 1 and R 2 groups were disregarded (Table 2). In the match pair analysis, the indazole was considered as reference structure. core in the ATP binding site similarly to the INDZ group, while major variability in the orientations of R1 and R2 groups was observed. Such conformational variability was mainly ascribed to: i) the inherent flexibility of the R1 appended moieties, ii) the open solvent-exposed binding site occupied by the R1 groups, and iii) the wide binding site filled by the R2 residues. Figure 3 highlights the good overlap of INDZ-IMID 1-IMID 2 triads bearing the same R1 and R2 groups.  (Table 2). For a comprehensive understanding of the IMID 1 and IMID 2 binding mode, the RMSD of the triplets was first calculated based on the common substructure of the core scaffold and carboxamide portion (excluding hydrogen atoms), while the more flexible R1 and R2 groups were disregarded (Table 2). In the match pair analysis, the indazole was considered as reference structure.  The new cores showed little displacement, in the range 0.15-0.67 Å RMSD, indicating very high similarity to each other in the binding site. For sake of completeness, RMSD of the common substructure of the whole molecules is also reported in Table 2, showing more variability within each triplet when R 1 and R 2 groups are considered.
In this work, MD was intended as a postprocessing tool to validate and refine Glide SP docking solutions. In addition, we aimed to evaluate the complex stability since co-crystal structures were unavailable at that time. The basic assumption was that a bad docking pose would have generated an unstable MD trajectory, during which the ligand could even leave the binding site. Conversely, a meaningful docking pose would have displayed stable and specific interactions with the target, showing a low RMSD over time, with respect to the starting conformation [33].
Based on the docking results, compounds seemed to reasonably occupy the ATP binding site of GSK-3β kinase. To confirm this hypothesis, we ran MD simulations of all derivatives using Desmond (see Materials and Methods). To define the optimum MD time scale, 500 ns simulations of 2, 8, and 14 in complex with GSK-3β enzyme were initially run. After a deep analysis, we found that the first 100 ns were representative of the overall 500 ns trajectories. In fact, the main structural events were generally preserved and the interaction fraction with the most relevant residues of the binding site were comparable ( Figure S1 in Supplementary Materials). Therefore, a 100 ns time scale was picked for all compounds as a good compromise between structural sampling thoroughness and computational cost. A single 100 ns MD simulation per compound was assumed to provide exhaustive conformational sampling given that the same combination of the R1 and R2 groups (considered the most fluctuating portions of the compounds) was simulated three times within each triplet.
Both protein and ligand RMSDs from the initial structures were calculated for each MD trajectory showing good structural stability of the whole systems ( Figure 4). In detail, the former spanned 0 to 2.66 Å and the latter 0 to 4.25 Å. Larger ligand RMSD values were observed in compounds bearing the R 1 piperidine group, which fluctuated in the solvent-exposed pocket. For each complex, RMSD of protein and ligand resulted to be quite comparable, meaning that the ligand thoroughly occupied the initial binding site without diffusing away. Ligand RMSF was also evaluated to investigate the most flexible moieties of the inhibitors ( Figure S2). In all systems, R 1 and R 2 groups showed the highest RMSF values compared to the stable core scaffolds.
Molecules 2020, 25, x 10 of 29 evaluated to investigate the most flexible moieties of the inhibitors ( Figure S2). In all systems, R1 and R2 groups showed the highest RMSF values compared to the stable core scaffolds. The pattern of H-bond interaction involving the hinge residues was essentially maintained along all the simulations. In detail, Val135 carbonyl interacted with carboxamide N-Hb via stable H-bond with comparable persistence between INDZ and corresponding IMID 1 and IMID 2 compounds along the MD trajectories. Conversely, the interaction occurrence of the central H-bond within each triad was variable, with INDZ derivatives contacting Val135 more persistently than IMID 1 and IMID 2 analogues (Table S1). The inner H-bond involving Asp133 was found in all systems although interaction distances and angles changed from INDZ to IMID 1 and IMID 2 core scaffolds (see below). The pattern of H-bond interaction involving the hinge residues was essentially maintained along all the simulations. In detail, Val135 carbonyl interacted with carboxamide N-Hb via stable H-bond with comparable persistence between INDZ and corresponding IMID 1 and IMID 2 compounds along the MD trajectories. Conversely, the interaction occurrence of the central H-bond within each triad was variable, with INDZ derivatives contacting Val135 more persistently than IMID 1 and IMID 2 analogues (Table S1). The inner H-bond involving Asp133 was found in all systems although interaction distances and angles changed from INDZ to IMID 1 and IMID 2 core scaffolds (see below). These results suggest that INDZ derivatives are more robustly anchored to the hinge residues than IMID 1 and IMID 2 analogues. Nevertheless, none of the simulations run on the three cores showed ligands moving away from the binding pocket. The amide carbonyl was exposed to water molecules through which it interacted with Ile62 and Gln185 to a different extent in all trajectories spanning from 10 up to 50% persistence ( Figure S3). Albeit not preserved during the entire simulations, this along with the outer H-bond seemed to stabilize the amide group preventing its rotation. Different results were obtained for the analyzed R 1 and R 2 groups. In 1-7-13 and 2-8-14 systems, the highly flexible piperidine group fluctuated in the open solvent-exposed pocket and contacted Gln72 and Pro136 via intermittent direct or indirect interactions ( Figure S3). A water-bridged H-bond involving the R 1 oxanyl ring was mostly observed with Pro136 in the remaining systems ( Figure S3). Other solvent-mediated contacts were found along all MD simulations between the R 1 group and Tyr134, Glu137, and Arg141. Overall, such indirect bonds randomly observed along the simulations reflected the high flexibility of the R 1 groups and the open binding pocket in which they lied. On the other side, aromatic and heteroaromatic rings variably interacted with the surrounding residues. In the 1-7-13 triad, the pyridine system stably occupied the region around the positively charged Lys85 residue, but optimal geometry for a stable contact was found to a variable extent in the three systems. Pi-pi interaction between the heteroaromatic ring and Phe67 was detected along the trajectories moving the nitrogen atom away from Lys85. The meta-pyridinyl substituent of 2, 8, and 14 derivatives was found to engage Asp200 from the highly conserved motif DFG via a water-mediated H-bond. Consequently, the pyridinyl ring fairly rotated to ease the rearrangement of the 5-substituent preventing a stable contact of its nitrogen atom with Lys85. De facto, this charged-assisted H-bond occurred no more than 30% of the simulation time in 2 and 8 trajectories and was not found in the 14 system. Moving to derivatives 3 and 15, the unsubstituted pyridinyl ring preserved the contact with the positively charged Lys85 for over 70% MD simulations. Instead, rotation of pyridine ring in 9 moved away the nitrogen atom resulting in 50% occurrence of such an H-bond. The R 2 group was also placed towards Phe67 in a pi-pi stacking interaction. Next, compounds 10 and 16 showed none or weak interaction with Lys85, while the isopropoxyl group confirmed the anticipated interaction with Asp200 via a water molecule. The corresponding INDZ derivative 4 contacted Lys85 more persistently and the oxygen atom of the meta substituent underwent also an indirect interaction with Gln185 ( Figure S3). The indirect contact with Asp200 occurred around 50% of the simulation time in all systems with the isopropoxyl group. Compared to the interaction pattern observed in unsubstituted pyridinyl derivatives, this additional bond could contribute to ameliorated affinity for GSK-3β. In the 5-11-17 triad, the R 2 hydroxymethyl group contacted Gln185 through a direct or water-mediated H-bond ( Figure S3). An indirect interaction of the ortho substituent with Asn186 was also observed. In addition, Lys85 occasionally contacted the pyridinyl nitrogen of 5. Finally, the 6-12-18 triad mostly orientated the halogen atoms of the di-F-phenyl group towards Lys85 so promoting an electrostatic-like interaction ( Figure S3).
In 1-18 systems, hydrophobic contacts were also observed where ligands engaged Ile62, Val70, Ala83, Leu132, Leu188, and Tyr134 at different extents. A timeline representation of the contacts (H-bonds, hydrophobic, ionic, water bridges) with protein is reported in Figure S4.
On the whole, ligand docking and MD simulations provided valuable structural information about the new IMID 1 and IMID 2 molecules. The initial binding poses appeared to be stable with preserved contacts across the MD trajectories, with subtle diversity of central and deep H-bonds. Thus, the newly designed compounds were basically suggested as comparable to the previous INDZ analogues.

In Vitro Data of Novel IMID 1 and IMID 2 Derivatives
Driven by the promising results, compounds 7-18 were synthesized as reported in Chemistry Section 2.1. Their biological activity was determined using an in vitro GSK-3β assay described elsewhere [17,19]. It followed that 7-18 exhibited good inhibition activity mostly in the nanomolar scale with an approximate 8-to 50-fold increase in IC 50 with respect to 1-6 ( Table 1). In vitro MDCK-MDR1 permeability of IMID 1 and IMID 2 derivatives was also experimentally assessed (Table 1). Derivatives 12, 15, and 18 resulted in improved permeability compared to their INDZ analogues. In this respect, a decrease of TPSA was observed moving from INDZ 6 to the corresponding analogues 12 (IMID 1) and 18 (IMID 2), and from 3 to 15. The SAR analysis highlighted the detrimental impact of the positively charged piperidine ring of 2-8-14 triplet and derivative 13 on CNS permeability. In fact, these compounds showed low P app values. When the piperidine was replaced with an oxanyl ring, while keeping the isopropoxypyridinyl group on the other side (2 vs . 4 and 8 vs. 10), a low CNS penetration was still observed, indicating a coexisting unfavorable contribution (high TPSA) of the R 2 moiety to permeability. The comparison of the 4-10 pyridinyl derivatives with the corresponding 6-12 di-F-phenyl molecules confirmed an ameliorated BBB permeability when lipophilicity raises.
In the new series of compounds, a general reduction of the ER values was also observed. This was consistent with published guidelines according to which removal of H-bond groups or reduction of H-bond strength can reduce P-gp efflux [34]. Moreover, it is known that compounds with the most basic pK a of less than 8.0 have a much higher likelihood of being identified as a non-P-gp substrate [34]. This unambiguously emerged from the measured ER of 8, 13, and 14 and by comparing derivatives 8 and 10 and 13 and 15.
Hence, the optimization of physicochemical properties and the introduction of focused chemical modifications to contain the critical HBD count of IMID 1 and IMID 2 series enabled to successfully overcome the indazole issues for CNS penetration.

Retrospective In Silico Analyses and NMR Studies to Rationalize in Vitro Data
This far, ligand docking and MD simulations were described driving the INDZ optimization campaign in a predictive way. To strengthen the results of the calculations, we generated a high-resolution X-ray structure of 16 in complex with the protein (PDB entry 6Y9S). In line with our prediction, 16 confirmed to inhibit GSK-3β enzyme by binding within the ATP pocket. The tridentate interaction with the hinge residues, the oxanyl ring towards Arg141, R 2 group contacting Lys85, and a pool of water molecules towards the open cavity and the back pocket were found. Compared to the docked structure, a very good consensus was observed with core RMSD = 0.35 Å (Figure 5), thus validating our models. The experimental binding mode also confirmed the stable water-bridged contact of the isopropoxyl group with Asp200 detected in the MD simulation ( Figure S3).
Thus, we carried out other in silico analyses starting from MD simulations to elucidate retrospectively the factors responsible for the unanticipated potency decrease of the new series. Due to the highly similar inhibition activity, the following is not meant to provide a rank of IMID 1 and IMID 2 compounds. Our results are intended to compare the new derivatives with the INDZ analogues instead, regardless of the tiny IC 50 differences within each core series.
First, the average binding free energy (∆G_bind) based on the molecular mechanics-generalized born surface area (MM-GBSA) was calculated from each production trajectory [35,36]. Interestingly, the INDZ molecules always resulted in better MM-GBSA scores compared to the corresponding IMID 1 and IMID 2 derivatives, in agreement with the in vitro potency (Table 3). In fact, energetics terms turned out to be generally more favorable for the INDZ compounds within each triplet.
Molecules 2020, 25, x 12 of 29 a pool of water molecules towards the open cavity and the back pocket were found. Compared to the docked structure, a very good consensus was observed with core RMSD = 0.35 Å (Figure 5), thus validating our models. The experimental binding mode also confirmed the stable water-bridged contact of the isopropoxyl group with Asp200 detected in the MD simulation ( Figure S3). Thus, we carried out other in silico analyses starting from MD simulations to elucidate retrospectively the factors responsible for the unanticipated potency decrease of the new series. Due to the highly similar inhibition activity, the following is not meant to provide a rank of IMID 1 and IMID 2 compounds. Our results are intended to compare the new derivatives with the INDZ analogues instead, regardless of the tiny IC50 differences within each core series.
First, the average binding free energy (ΔG_bind) based on the molecular mechanics-generalized born surface area (MM-GBSA) was calculated from each production trajectory [35,36]. Interestingly, the INDZ molecules always resulted in better MM-GBSA scores compared to the corresponding IMID 1 and IMID 2 derivatives, in agreement with the in vitro potency (Table 3). In fact, energetics terms turned out to be generally more favorable for the INDZ compounds within each triplet.   In both the IMID series, we supposed the drop in potency to be attributable to the less acidic C-Ha group of the central core. Together with other moieties, it tightly anchors the compounds in the binding pocket, but the nature of the H-bond is weaker compared to the one involving conventional INDZ N-Ha. The general characteristics of this type of interaction to take into account are: (i) the distance between the hydrogen and the acceptor (A), (ii) the angle formed by the donor (D), the hydrogen itself and A, and (iii) the strength [37,38]. To different levels, these components were herein investigated across the molecules.
We examined how the binding mode changes as the core scaffolds are varied, in terms of distance between each hinge binding atom of the ligands and the partner heavy atom of Asp133 and Val135 along the MD simulations. The analysis revealed that moving from INDZ to IMID 1 and IMID 2 analogues, an elongation of the deep H-bond occurred, while the outer and central interactions remained approximately constant in length (red versus green and orange scatter plots in Figure 6). The stretching of the deep H-bond was found in IMID 1 and IMID 2 compounds of all triplets regardless of the R 1 and R 2 moieties. This suggests a possible rotation of the ligand about the hinge binding residues which might affect the SAR and elucidate the potency reduction. As anticipated in the MD discussion, persistence percentage of the central H-bond is higher in INDZ compounds compared to IMID 1 and IMID 2 analogues. The optimally spaced deep H-bond of the INDZ system might reflect on such high persistence. Conversely, the central H-bond is more scattered along the simulations of IMID 1 and IMID 2 derivatives as a consequence of the elongated interaction of the C-Ha.
investigated across the molecules.
We examined how the binding mode changes as the core scaffolds are varied, in terms of distance between each hinge binding atom of the ligands and the partner heavy atom of Asp133 and Val135 along the MD simulations. The analysis revealed that moving from INDZ to IMID 1 and IMID 2 analogues, an elongation of the deep H-bond occurred, while the outer and central interactions remained approximately constant in length (red versus green and orange scatter plots in Figure 6). The stretching of the deep H-bond was found in IMID 1 and IMID 2 compounds of all triplets regardless of the R1 and R2 moieties. This suggests a possible rotation of the ligand about the hinge binding residues which might affect the SAR and elucidate the potency reduction. As anticipated in the MD discussion, persistence percentage of the central H-bond is higher in INDZ compounds compared to IMID 1 and IMID 2 analogues. The optimally spaced deep H-bond of the INDZ system might reflect on such high persistence. Conversely, the central H-bond is more scattered along the simulations of IMID 1 and IMID 2 derivatives as a consequence of the elongated interaction of the C-Ha. In the MD trajectories, optimal D-H-A angles promoting strong interaction of the INDZ N-Ha with the hinge portion were also observed, while the values lowered moving to IMID 1 and IMID 2 cores (Figure 7). Based on the observed arrangements, we may ascribe the H-bond between IMID C-Ha and Asp133 carbonyl group to a weakened interaction compared to INDZ N-Ha, so substantiating the affinity decrease of IMID 1 and IMID 2 derivatives versus the GSK-3β enzyme.
Lastly, the H-bond acidity prediction was investigated via QM methods using Jaguar software [39]. The Vα(r) descriptor of the core Ha and the carboxamide Hb was calculated. This is an effective and useful predictor of H-bond strength defined in Peter Kenny's work and recently implemented in Schrödinger suite [40]. Larger values correlate with stronger HBDs, while HBAs are characterized by Based on the observed arrangements, we may ascribe the H-bond between IMID C-Ha and Asp133 carbonyl group to a weakened interaction compared to INDZ N-Ha, so substantiating the affinity decrease of IMID 1 and IMID 2 derivatives versus the GSK-3β enzyme.
Lastly, the H-bond acidity prediction was investigated via QM methods using Jaguar software [39]. The Vα(r) descriptor of the core Ha and the carboxamide Hb was calculated. This is an effective and useful predictor of H-bond strength defined in Peter Kenny's work and recently implemented in Schrödinger suite [40]. Larger values correlate with stronger HBDs, while HBAs are characterized by low numbers. Although Kenny's analysis never really considered C-H donors, the principle of the calculation is still well defined. Here, all ligands from Glide SP docking were used as input for Vα(r) calculation. (Table 4). The calculation indicated that the INDZ core has the strongest H-bond potential in all compounds, with a very strong N-Ha donor. This illustrates the importance of this moiety at the hinge for good affinity, but it equally explains the reason why it is felt to be an issue for CNS penetration. The calculations also showed that the C-Ha of the IMID 1 and IMID 2 cores are weaker H-bond donors, with some variability between the two new series. In each triplet, there seemed to be a large Vα(r) difference when Ha was compared between INDZ derivatives and the corresponding IMID 1 and IMID 2 analogues (∆Vα(r) > 40). On the other hand, the Vα(r) values between related IMID 1 and IMID 2 compounds within the same triplet resulted in a minor gap. These in silico results appeared to be in line with the experimental IC 50 values, where potency reduction of one or more orders of magnitude was observed moving from INDZ compounds to the new series, while IMID 1 and IMID 2 derivatives were quite comparably potent as GSK-3β inhibitors. About Hb, no big changes were observed from INDZ to IMID 1 and IMID 2 scaffolds. This suggested that the acidic N-Ha/C-Ha group was a significant discriminant of the cores for inhibition activity, while the carboxamide N-Hb did not make a real difference. With this, we proved Vα(r) to be a valuable descriptor to explore SAR in terms of H-bond acidity.
These observations were further substantiated using 1 H-NMR characterization. By way of example, we selected the triplet 2-8-14. In particular, H-bond donating abilities of the core N-Ha/C-Ha and the carboxamide N-Hb of the three compounds were investigated by means of NMR spectroscopy.
In NMR experiments, proton chemical shifts depend upon many factors, among which there are H-bond interactions. In general, a proton involved in such interaction exhibits a decreased shielding and therefore tends to move downfield in the NMR spectrum.
First, the complete assignment of the 1 H-NMR spectra of 2, 8 and 14 compounds was carried out in CDCl 3 by means of the simultaneous 1 H-1 H COSY, 1 H-1 H NOESY, 1 H-13 C HSQC, and 1 H-13 C HMBC 2D NMR analyses (see Supplementary Materials), to unequivocally assign the Ha and Hb resonances.
Subsequently, the 1 H-NMR spectra were acquired in both CDCl 3 and DMSO-d 6 as "not good" and "good" H-bond acceptor solvents, respectively [41]. The variation of chemical shifts of Ha and Hb protons in the two solvents was calculated as ∆δ (δ DMSO − δ CDCl3 ), using the TMS signal as a reference at 0 ppm (Table 5 and Figures S5-S7). The target concentration of the compounds was selected by a serial dilution method (reaching a 0.3 mM final concentration), until a stable chemical shift difference between the two solvents was observed, and the solute-solute interactions were expected as low as possible. The large chemical shift difference of 3.4 ppm observed for proton Ha in compound 2 going from DMSO-d 6 to CDCl 3 suggested that this proton was involved in H-bond with the former solvent, excluding the formation of an intramolecular H-bond. On the contrary, the decrease in ∆δ observed for proton Ha of compounds 8 and 14 suggested a lower acidity of the IMID 1 and IMID 2 core scaffolds. This is in agreement with the pharmacological results where potency reduces moving from the INDZ compounds to the new series, whereas IMID 1 and IMID 2 derivatives are similarly double-digit nanomolar GSK-3β inhibitors. On the other hand, no large differences were observed for the ∆δ of Hb within the triplet (Table 5), suggesting comparable acidity of the carboxamide N-Hb in the three different scaffolds. Moreover, the sharper and better resolved triplet observed for the Hb signals moving from CDCl 3 to DMSO-d 6 ( Figures S5-S7) confirmed the involvement of the protons in H-bond interactions with the latter solvent, which prevented them from engaging in exchange processes with the residual amounts of water in solution.
Overall, preliminary NMR data suggested, as expected, that the INDZ core has the strongest H-bond potential playing a crucial role in the affinity for GSK-3β binding pocket. Moreover, Vα(r) was in accordance with ∆δ, as a significant acidity drop of Ha was observed moving from INDZ to IMID 1 and IMID 2 in both calculated and experimental parameters.
The proton spectra of the 2, 8 and 14 were also acquired in deuterated DMSO at five different temperatures (298, 303, 308, 313, and 318 K). The temperature coefficients (T coeff ) were calculated for N-Hb protons by linear fitting of the chemical shift as function of temperature curves (as exemplified in Figure 8 for compound 14).
in accordance with Δδ, as a significant acidity drop of Ha was observed moving from INDZ to IMID 1 and IMID 2 in both calculated and experimental parameters.
The proton spectra of the 2, 8 and 14 were also acquired in deuterated DMSO at five different temperatures (298, 303, 308, 313, and 318 K). The temperature coefficients (Tcoeff) were calculated for N-Hb protons by linear fitting of the chemical shift as function of temperature curves (as exemplified in Figure 8 for compound 14). In general, amide proton chemical shifts show a strong temperature dependence and this is mainly due to the occurrence of H-bonds [42]. Chemical shifts depend on the inverse third power of the distance between amide proton and H-bond acceptor; the average H-bond length increases with temperature because of the increased thermal motions [43]. Therefore, the deshielding effect induced by an acceptor decreases at higher temperature and amide protons result in being shifted upfield. In line with this general trend, N-Ha and N-Hb chemical shifts followed this behavior and showed a linear dependence with temperature, allowing the determination of Tcoeff values from the slope of the curves (Table 6).  Hb atoms have shown similar Tcoeff, more negative than that of the Ha proton in the INDZ core. Indeed, for weaker H-bonds, structural perturbations are more pronounced and N-Hb chemical shifts are more temperature dependent (more negative). In previous studies, Tcoeff values in the range of -12.3 to −0.5 ppb/K were reported for -OH groups which were solvent exposed and involved in In general, amide proton chemical shifts show a strong temperature dependence and this is mainly due to the occurrence of H-bonds [42]. Chemical shifts depend on the inverse third power of the distance between amide proton and H-bond acceptor; the average H-bond length increases with temperature because of the increased thermal motions [43]. Therefore, the deshielding effect induced by an acceptor decreases at higher temperature and amide protons result in being shifted upfield. In line with this general trend, N-Ha and N-Hb chemical shifts followed this behavior and showed a linear dependence with temperature, allowing the determination of T coeff values from the slope of the curves (Table 6).  Hb atoms have shown similar T coeff , more negative than that of the Ha proton in the INDZ core. Indeed, for weaker H-bonds, structural perturbations are more pronounced and N-Hb chemical shifts are more temperature dependent (more negative). In previous studies, T coeff values in the range of -12.3 to −0.5 ppb/K were reported for -OH groups which were solvent exposed and involved in intermolecular H-bonds [44], while a mean value equal to −7.1 ± 2.5 ppb/K was determined for the solvent-exposed amide protons of the protein backbone [42]. In our case, the obtained values reflected the presence of weak intermolecular H-bonds in DMSO-d 6 ; the more negative T coeff values for the Hb atoms suggested that, by increasing the temperature, intermolecular solute-solvent H-bonds are broken more easily for N-Hb than for N-Ha within the INDZ core, in agreement with the conclusions drawn in light of the ∆δ values.
From the ∆δ downfield shifts reported above, it was possible to estimate the H-bond energies E HB . The most universal approach is the empirical estimation of the H-bond strength by use of the Equation (3) (see Materials and Methods) [45]. The equation is only valid for weak and moderate values of the H-bond energy. The energy varies within a wide range from weak (1-4 kcal/mol) to medium (4-15 kcal/mol) to strong H-bonds (15-40 kcal/mol).
Based on the calculated values ( Table 7), Ha of compound 2 showed the highest energy, while comparable and lower values were observed concerning Hb, regardless of the central core. Accordingly, we can assert that Ha of INDZ and Hb of all core scaffolds are involved in intermolecular H-bonds that can be classified as weak/moderate to highly weak, respectively. Again, this outcome enables to rationalize SAR confirming the key role of the INDZ N-Ha for the high inhibition activity.

Chemistry
Reagents were purchased from Sigma-Aldrich (Madrid, Spain) and Fluorochem (Hadfield, UK) and used as received.
Final compound purity (always >95%) was determined by HPLC-UV, by comparing the area of the target peak and total peak areas of minor impurities. HPLC-UV analyses were performed with a pump/autosampler (Waters 2695, Alliance model), a UV photo diode array detector (Waters 2996 model), and a Waters system data management (Empower 2) (Etten-Leur, The Netherlands). The column used was generally X-Bridge C18 (150 × 4.6 mm, 5 µm).
UPLC-QTOF analyses for HRMS determination of final compounds were carried out by means of a SYNAPT MS-ACQUITY UPLC system (Waters). The system was operated in positive ion mode in the "V-Optics" configuration. Leucine-enkephalin (200 pg/µL) was employed as the lock mass in order to provide authenticated exact mass measurement in MS and MS/MS modes within 5 ppm RMS mass accuracy. The column was an Acquity BEH C18 (50 × 2.1 mm, 1.7 µm).
Nuclear Magnetic Resonance Spectroscopy ( 1 H-NMR and 13 C-NMR) were obtained using Bruker Advance systems (Milan, Italy), operating at 300, 400 and 500 MHz. Spectra were acquired at 300 K, using DMSO-d 6 , CD 3 OD, and CDCl 3 as solvents. All resonance bands were referenced to tetramethylsilane (internal standard). Compounds were purified with one of the following techniques: (i) flash chromatography on silica gel (Grace Reveleris flash chromatography system with a 40 µM silica cartridge, flow = 60 mL/min), with an appropriate gradient of mixtures of DCM/MeOH or Cyclohexane/EtOAc as eluents; (ii) preparative HPLC-MS system (sample was dissolved at 50 mg/mL in DMSO/CH 3 CN in 1:1 ratio), using an appropriate gradient of CH 3 CN and H 2 O + 0.1% NH 4 HCO 3 ; flow =40 mL/min. The preparative HPLC-MS system consisted of a Waters 2767 Sample manager, a Waters 2478 dual λ absorbance detector, and a Waters Micromass ZQ single quadrupole mass spectrometer with an electrospray ionization (ESI) source. The column used was a X-Bridge Prep C18 5 µm OBD with 19 × 10 mm (Waters) precolumn. Fraction collection was available from the system software MassLynx v. 4.1 (Etten-Leur, The Netherlands). Detection wavelength was set to 230 nm and temperature to 25 • C.

General Procedures (A)-(D):
General amide coupling procedure (A) for the synthesis of compounds 21, 22 and 13-18. To a solution of carboxylic acid/lithium carboxylate salt (1.0 eq.) in DMF (0.07 M) HATU (1.2 eq.), DIPEA (4.5 eq.) and the proper amine (1.0-1.2 eq.) were added. The solution was stirred at rt until full conversion was observed by HPLC-MS (24-72 h). The solvent was evaporated under reduced pressure, the residue was taken up with EtOAc and washed with 1N aqueous NaOH solution (×2) and brine. The combined organic phases were dried over anhydrous MgSO 4 , filtered and concentrated in vacuo. The crude was purified by flash chromatography (SiO 2 , DCM/MeOH for 21 and 22) or reverse phase chromatography (C18, NH 4 HCO 3 /ACN for 13-18) to give the target compound.
General boronate synthesis procedure (B) for the synthesis of compounds 23 and 25. A mixture of bromide (1.0 eq.), bis(pinacolato)diboron (2.0 eq.) and potassium acetate (2.0 eq.) in 1,4-dioxane (0. 12 M) was placed in a Schlenk tube. After degassing with N 2 for 10 min, Pd(dppf)Cl 2 (0.15 eq.) was added and the mixture was heated at 100 • C overnight. Conversion was checked by HPLC-MS. The mixture was diluted with EtOAc/MeOH, the insoluble solids were removed by filtration through Celite and the solvent was evaporated under reduced pressure. The crude was purified by flash chromatography (SiO 2 , DCM/MeOH) to give the target compound.
General ester hydrolysis procedure (D) for the synthesis of compounds 30-33.
To a solution of ethyl carboxylate ester (1.0 eq.) in MeOH and water (1:1, 0.03 M), lithium hydroxide hydrate (2.0 eq.) was added. Conversion was checked by HPLC-MS. After stirring at 50 • C overnight, the reaction mixture was concentrated to dryness under reduced pressure to give the target compound that was used in the next step without further purification.

MDCK-MDR1 Cell Line Assay
MDCK-MDR1 cells were obtained from the NIH (Rockville, MD, USA). Cells were seeded onto Millipore Multiscreen Transwell plates at 3.4 × 105 cells/cm2 and cultured in DMEM. Cell culture and assay incubations were carried out at 37 • C in an atmosphere of 5% CO2 with a relative humidity of 95%. The monolayers were prepared by rinsing both apical and basolateral surfaces twice with Hanks Balanced Salt Solution (HBSS) at the desired pH warmed to 37 • C. Cells were then incubated with HBSS at the desired pH in both apical and basolateral compartments for 40 min to stabilize physiological parameters. The dosing solutions were prepared by diluting test compound with assay buffer to give a final test compound concentration of 10 µM (final DMSO concentration of 1% v/v). Analytical standards were prepared from test compound DMSO dilutions and transferred to buffer, maintaining a 1% v/v DMSO concentration. Assay buffer was composed of supplemented HBSS pH 7.4. For assessment of A-B permeability, HBSS was removed from the apical compartment and replaced with test compound dosing solution. The apical compartment insert was then placed into a companion plate containing fresh buffer (containing 1% v/v DMSO). For assessment of B-A permeability, HBSS was removed from the companion plate and replaced with test compound dosing solution. Fresh buffer (containing 1% v/v DMSO) was added to the apical compartment insert, which was then placed into the companion plate. At 60 min, the apical compartment inserts and the companion plates were separated and apical and basolateral samples diluted for analysis. Test compound permeability was assessed in duplicate. Two control compounds were screened alongside the test compounds, propranolol (highly permeable), and prazosin (a substrate for P-gp). Test and control compounds were quantified by LC-MS/MS cassette analysis using an 8-point calibration with appropriate dilution of the samples. Cyprotex generic analytical conditions were used. The starting concentration (C 0 ) was determined from the dosing solution and the experimental recovery calculated from C 0 and both apical and basolateral compartment concentrations. The permeability coefficient (P app ) for each compound was calculated from the following Equation (1): where dQ/dt is the rate of permeation of the drug across the cells, C 0 is the donor compartment concentration at time zero, and A is the area of the cell monolayer. C 0 was obtained from analysis of the dosing solution. For bi-directional experiments, ER was calculated from mean A-B and B-A data. This was derived from Equation (2): ER value greater than two indicates that drug efflux is occurring.

Ligand Docking
The Schrödinger suite was used for all calculations. Ligands were converted to 3D structures and prepared with LigPrep tool using default setting [50]. Overall, the most likely tautomers and protonation states were considered for each ligand. Next, a conformational search was carried out using ConfGen resulting in the generation of 10 conformers per ligand [51]. This represented the input dataset for Glide SP docking [30,31,52]. The coordinates of the crystal structure of GSK-3β in complex with 2 (PDB entry 6Y9R) were used. The protein was prepared using Maestro Protein Preparation Wizard [53]. Hydrogen atoms were added, and water molecules were removed from the protein structure. The default setting was used for receptor grid preparation. The binding site was specified as a box of 10 × 10 × 10 Å 3 centered on the centroid of the co-crystallized inhibitor. In addition, we treated thiol and hydroxyl groups of Tyr134, Thr138, and Cys199 as rotatable (flexible) groups. Three H-bond constraints were defined being the carbonyl oxygen of the hinge residues Asp133 and Val135, and the hydrogen of the N-H group of Val135. In Glide SP docking, we required that at least one such H-bond crucial for compound activity was formed. We also used the enhanced sampling option to enhance the sampling of conformational space by adding variations on the input structure to the conformational search and included aromatic ligand hydrogen atoms as donors. Visual inspection of the docking results was carried out enabling the picking of one pose per ligand as a starting point for further studies.
A rigid-body superposition of INDZ vs. IMID 1 and IMID 2 analogues within each triplets was performed using Maestro software [54], and RMSD was computed ( Table 2). The common substructure was defined considering all heavy atoms with the exception of the ones of the bicyclic system that are different moving from INDZ to IMID 1 and IMID 2 analogues. RMSD was computed without changing structure positions. Figures 3 and 5 were prepared with PyMOL 1.8.6.2 [55]. Crystal structure of 16 in complex of GSK-3β ( Figure 5) was prepared using a default setting of the Maestro Protein Preparation Wizard tool [53].

MD Simulations
The selected binding poses of 1-18 within the ATP binding site of GSK-3β were processed in MD simulation using the graphics processing unit (GPU) tool of the Desmond MD package as supplied with the Schrödinger 2019-1 suite [32,56]. The Desmond System Builder was used to create a fully solvated system with simple point charge (SPC) waters and a cubic solvation box with a 10 Å buffer around the original complex. The system was electrostatically neutralized by adding Na + ions in defined positions by using the ion placement option. Default parameters of OPLS-3e force field were used. The prepared systems included approximately 58,900 atoms. The MD simulation protocol followed an extended relaxation schedule consisting of the following steps:

NMR Study
NMR experiments were performed in deuterated DMSO and CDCl 3 solutions and recorded with a Bruker Avance spectrometer operating at 500 MHz for 1 H and 125 MHz for 13 C, equipped with a 5 mm inverse broad-band probe and z-axis gradients.
The complete NMR spectra assignment of the three compounds was performed in CDCl 3 solutions at 298 K. COSY experiments were acquired with a spectral window of 15 × 15 ppm (carrier frequency at 4.5 ppm) using 2048 × 256 data points, 16 transients, with a relaxation delays of 2 s. NOESY experiments were acquired with a spectral window of 15 × 15 ppm (carrier frequency at 4.5 ppm) using 4096 × 256 data points, 32 transients, with a relaxation delay of 2 s and a mixing time of 600 ms. 1 H-13 C HSQC experiments were acquired with a spectral window of 15 ppm × 200 ppm (carrier frequencies at 6.5 and 100 ppm) using 2048 × 256 data points, 16 transients with relaxation delays of 2 s. 1 H-13 C HMBC experiments were acquired with a spectral window of 15 ppm × 220 ppm (carrier frequency at 6.5 and 110 ppm) using 4096 × 512 data points, 32 transients with relaxation delays of 2 s. All data were processed with TopSpin using a sine squared window function for each dimension of the 2D data.
The variations of chemical shift (∆δ) for protons of interest were calculated from the difference of the chemical shift in DMSO-d6 and in CDCl 3 , detected in the 1 H-NMR spectra at 298 K.
For T coeff determination, 1 H-NMR spectra of the three compounds were performed in deuterated DMSO at five different temperatures: 298 K, 303 K, 308 K, 313 K, and 318 K.

Conclusions
In this study, novel derivatives around the lead series core scaffold were investigated by means of ligand docking and MD simulations. The novel IMID derivatives turned out to be GSK-3β inhibitors in the nanomolar scale with a potency drop compared to the previous INDZ series but improved CNS permeability. Interestingly, the in silico-derived models of the binding mode generated during the design process well suited the crystallographic data (PDB entry 6Y9S). In order to rationalize SAR, computational tools were also employed retrospectively. In particular, MD and QM-based calculations shed light on the differences of the H-bond interactions between the core scaffolds and the hinge residues moving from INDZ to IMID 1 and IMID 2 derivatives. The pharmacological profile of IMID 1 and IMID 2 derivatives was associated to the lower acidity of the C-Ha versus N-Ha group of the core scaffolds. NMR-based studies were also carried out to substantiate this hypothesis. Overall, our findings provided adequate evidence that low C-Ha acidity and consequent unfavorable geometric rearrangement affected H-bond donating abilities so reducing the affinity for GSK-3β enzyme. This represented significant structural activity insight into our GSK-3β inhibitor design and demonstrated the pivotal role of computational approaches properly supported by experimental proofs in the different steps of the drug discovery program. In the future, prospective application of the aforementioned in silico tools may lead to further optimized classes of GSK-3β inhibitors.
Supplementary Materials: The following are available online, Figure S1: Ligand-Protein interactions monitored throughout 100 ns (a) and 500 ns (b) simulations, Figure S2: Ligand RMSF from each MD simulation, Figure  S3: representative snapshots of MD simulations, Figure S4 Table S1: Percentage of interaction of the outer and central H-bond across each MD trajectory, Table S2: Data collection and model statistics of crystal structure of 2 in complex with GSK-3β enzyme (PDB entry 6Y9R), Table S3