Selected arylsulphonyl pyrazole derivatives as potential Chk1 kinase ligands-computational investigations.

Protein kinases control diversity of biochemical processes in human organism. Checkpoint 1 kinase (Chk1) is an important element of the checkpoint signalling pathways and is responsible for DNA damage repair. Hence, this kinase plays an essential role in cancer cells survival and has become an important target for anticancer agents. Our previous investigations showed that some arylsulphonyl indazole derivatives displayed anticancer effect in vitro. In the present study, in order to verify possibility of interactions of pyrazole and indazole derivatives with Chk1, we focused on the docking of selected tosyl derivatives of indazole and condensed pyrazole 1-7 to the Chk1 pocket, analysis of interactions involving optimized ligand-protein system using DFT formalism, and estimation of the interaction enthalpy of the ligand-protein complex by applying the PM7 method. The estimation of binding affinity seems to indicate that the indazole 5-substituted with 3,5-dimethylpyrazole 4 and condensed pyrazoloquinoline derivative 7 fit the best to the Chk1-binding pocket. The values of the energy of interaction, i.e. the enthalpy change (ΔHint), were between - 85.06 and - 124.04 kcal mol-1 for the optimized ligand-Chk1 complexes. The relaxation of the ligands within the complexes azole-protein as well as the distribution of hydrogen contacts between the ligands and kinase pocket amino acids was also analysed using molecular dynamics as a supporting method. Graphical Abstract Presentation of methods used to describe the interactions between arylsulphonyl pyrazole derivatives and Chk1 kinase.

Introduction

Chemotherapy is one of the major treatments of cancer and often relies on DNA damage of fast proliferating cancer cells. The cancer cell response to such damage (DNA damage response, DDR) is activation of the S- and G2/M-Phase checkpoints that leads to effective repair of DNA and cell survival. The above checkpoints are controlled by checkpoint kinase 1 (Chk1), a serine/threonine kinase with self-protection function [1]. Tumour cells have developed a strong dependence on Chk1 for survival, so this kinase can be a promising target for anticancer agents. Inhibition of Chk1 results in the checkpoint abrogation and suppression of DNA repair. Therefore, the Chk1 inhibitors can increase the therapeutic effect of DNA-targeting anticancer agents and are typically used in combination with radiotherapy or DNA-damaging drugs [2].

Overexpression of Chk1 has been found to relate to pulmonary arterial hypertension (PAH), a progressive vascular remodelling of distal arteries, resulting in severe elevation of pulmonary artery pressure [3]. It has been shown that smooth muscle cells of these arteries exhibit cancer-like properties, i.e. fast proliferation and resistance to apoptosis. Hence, Chk1 might be an attractive target for the therapy of patients with PAH.

Chk1 is also involved in DNA repair by targeting repair kinases, like DNA-pyruvate kinase (DNA-PK), as well as some repair pathways, e.g. BRCA-mediated DNA repair [1]. Moreover, the ATR/Chk1 signalling pathway suppresses both a caspase-3-dependent and caspase-2-dependent apoptotic response; thus, there is supposedly a direct link between Chk1 and apoptosis [1].

Pyrazole, indazole, and condensed pyrazole derivatives constitute an important class of kinase inhibitors [4-12] including Chk1 blocking agents [6, 13, 14]. Our previous investigations showed that some indazole derivatives with an arylsulphonyl substituent at the position 3 of indazole displayed anticancer effect against the colon cancer cell line HT29 as well as breast cancer cell lines MCF7 and MDA-MB-231 [15, 16]. These derivatives can interact with magnesium ions, important for carcinogenesis [17, 18], as well as with several amino acids, present in the kinases binding site [19]. We suggested that the activity mechanism might involve either stimulation of proapoptotic proteins or inhibition of signalling pathways proteins, but we were unable to specify the target. To verify the above hypotheses, we decided to study interactions of indazole derivatives 1–7 with Chk1. The study focused on the docking of selected tosyl derivatives of indazole and condensed pyrazole to the Chk1 pocket with the use of AutoDock Vina suite, analysis of interactions involving optimized ligand–protein system with the help of DFT formalism, and estimation of the interaction enthalpy of the ligand–protein complex (PM7 method). For the analysis of the ligand relaxation within the azole–protein complexes, we used the molecular dynamics method only as a supporting technique.

Considering the potential affinity of pyrazole and indazole derivatives to kinases, we decided to investigate the interactions of azoles 1–7 with Chk1 using a Protein Data Bank deposit 2e9n.pdb [20, 21] in complex with A767085 ligand.

The clinically approved kinase inhibitors bind to the catalytic kinase domain—the ATP-binding site. All protein kinases, including Chk1, share the same catalytic domain that consists of an N-terminal lobe, constructed of a five-stranded β-sheet and a single α-helix, and a C-terminal lobe, mainly α-helical [6, 7]. The ATP-binding site forms a cleft between these two lobes and is composed of five regions, important for small molecule inhibitor binding, namely a linker (hinge) region for adenine, ribose pocket, phosphate binding loop (P-loop)-catalytic aspartate region, back hydrophobic (water) pocket, and front specificity pocket.

The linker region, a short, mostly hydrophobic, strand connecting C and N lobes, interacts with the adenine ring of ATP through the key site residues, i.e. glutamic acid E85, tyrosine Y86, and cysteine C87 (the numbering refers to the Chk1 sequence).

The P-loop interacts with the phosphate group of ATP through a glycine rich motif. The catalytic aspartate fragment at the active site gate contains a conserved Asp-Ph-Gly (DFG) motif at its N-terminal edge. The DFG motif adopts normally two conformations, namely DFG-in and DFG-out. In the first conformation, the aspartic acid side chain Asp148 is directed towards the active site and coordinates magnesium. As this orientation allows catalysis to proceed, it is called the active conformation, as opposed to the inactive DFG-out conformation in which the Asp148 side chain is pointed away from the active site. Most inhibitors interact with the active DFG-in orientation.

The ribose pocket in the neighbourhood of the linker region contains glutamic acid E91 that forms important contacts with the ribose hydroxylic groups. The back hydrophobic pocket is usually occupied by water molecules. The entry to this pocket is composed of the gatekeeper residue L84. The front specificity pocket is a relatively small hydrophobic region between the linker site and a hydrophilic, solvent-exposed sector of the protein [6].

Computational methods

For the initial preparation of the analysed ligands, we obtained 1000 conformations of azoles 1–7 (Scheme 1) using the Gabedit 2.4.7 program [22]. The following parameters were applied for energy minimization: heating—0.5–1000 ps, 0–1000 K, and Andersen thermostat [23] for conformation creation and the Verlet velocity algorithm for MD Trajectory [24]. For the conjugate gradient, we applied the Fletcher–Reeves method with maximum line searches of 25 [25]. In the next step, we employed minimization with molecular mechanics (MM), Amber force field with terms for the bond stretch, angle bend, torsion, nonbonded, and H-bonded. Then, all the resulting conformations were optimized with PM7 (Mopac 2016) [26, 27]. Each of the most energetically stable conformers of hetarenes 1–7, i.e. with the lowest final heat of formation (HOF), was refined using density functional theory formalism [28] in the gaseous phase. The DFT calculations were executed, and geometries of each previously pre-optimized conformers of 1–7 (Scheme 1) were further refined using the Gaussian 16 A.03 program [29] at the B3LYP/6-31G(d,p) level of theory.

Scheme 1

The investigated azoles as potential Chk1 kinase ligands

The human Chk1 kinase protein in complex with the A767085 ligand, acquired from the Protein Data Bank base (PDB entry: 2e9n.pdb), was selected as a biological target [20, 21]. An initial target for further optimization was prepared by removing the internal A767085 ligand, and all water molecules from the 2e9n.pdb file but the internal coordinates were kept unchanged. The genetic algorithm (GA) method implemented in the AutoDock Vina program [30] was employed to provide the appropriate binding orientations and conformations of the compounds in the Chk1-binding pocket. Polar hydrogen atoms were added, and partial charges were assigned to the protein. Then, the internal ligand was replaced by the optimized structure hetarenes 1–7, and additionally, the residues were saturated with hydrogen atoms (an example of configuration file used for docking protocol of azole 1 is given in Table S1 in the Supplementary material; Cartesian coordinates of the lowest energy poses of all docked azoles 1–7 are given in Figure S1 in the Supplementary material). A grid box of 10-Å size (centre _x = − 1.0, centre_y = 10.0, centre_z = − 19.0) was defined to carry out the docking simulation. The outputs (*.pdbqt files) after docking procedure were visualized with the Chimera 1.13.1 package [31]. The projections of the first poses of azoles 1–7 docked to the Chk1 pocket (Fig. 1) were visualized with MGLTools 1-5-6 program [32]. The Chemcraft 1.7 software was utilized for the visualization of all poses shown in Figure S1 [33]. For the semiempirical calculations with the use of PM7 method [34], we employed the Mopac 2016 software [26]. For the initial optimization of hydrogen atoms, we used keywords “PM7 XYZ GNORM=5 T=14D DUMP=1200 NOOPT OPT-H MOZYME”; then, we carried out the refinement of complexes using “PM7 XYZ GNORM=1 T=14D DUMP=1200 MOZYME” keywords. For the molecular dynamics (MD) calculations, the GROMACS 2016.4 [35, 36] was employed to simulate the solvated complexes. The Amber99SB-ILDN force field [37] was used to parameterize the protein and counter ions. The general GAFF force field [38] was utilized to represent the ligands and their topology with the help of Topolbuild 1.2.1 [36]. Finally, the complexes were inserted into the cubic water boxes using the TIP3P water model [39] (10 × 10 × 10 nm). The soluble complex consisted of one molecule of Chk1 kinase, one ligand molecule, approximately 31,397 water molecules, and about six Na+ ions depending on the charge of the ligand. The soluble complexes were first minimized using the steepest descent scheme. Then, the minimized configurations were relaxed in NVT and NPT ensembles with 500-ps MD length per simulations. The complexes were restrained by NVT simulations using a small harmonic force. For the complexes free of restraints, we adopted NPT MD simulations. The relaxed system was then used as an initial conformation for 20-ns MD simulations. The time step used throughout the MD calculations was 2 fs.

Fig. 1

Superimposition of docked azoles 1–7 (1st poses, MGLTools 1-5-6 program)

Results and discussion

Ligand docking

Indazole derivatives 1–5, as well as condensed pyrazoles 6 and 7, were initially optimized by molecular mechanics; then, a series of 1000 conformers for each azole derivatives were further refined by the semiempirical method PM7. Next, for each of the azole derivatives 1–7, a conformer with the smallest heat of formation (HOF) was selected and optimized in the SCF procedure (DFT method, B3LYP functional, Gaussian 16 A.03 program [29]). The optimized ligands 1–7 were docked to the 2e9n.pdb protein using the AutoDock Vina suite [30]. Nine poses were obtained for each of azole 1–7 from which the first poses (Fig. 1; Figures S1–S22, Supplementary material) have the lowest negative value of binding affinity (Table S1, Supplementary material). It should be noted that azoles 1–4 and 7 adopt similar docking poses, whereas for azole 6, only the tosyl group assumes a similar orientation in the binding site. Compound 5 differs from the remaining azoles—its first pose is not superimposable on the poses for azoles 1–4 and 6–7.

Next, the molecular electrostatic potential (MEP) was determined by the B3LYP/6-311++G(2d,3p) approach for the conformers of azoles 1–7 (1st poses) with geometry previously optimized at B3LYP/6-31G(d,p) level of theory in gaseous phase (Gaussian 16 A.03 program [29], Figures S1–S7, Supplementary material). In our investigations, involving the multilevel approach to the conformational rotamers search, the results were refined using a basis set enriched with the higher-level polarization functions.

The obtained results show that the pyrrolic nitrogen of the pyrazole ring, an H-bond donor, and the electron-withdrawing tosyl substituent are the most important for the azole–protein interactions. Moreover, the pyridinic nitrogen of the pyrazole ring as well as the same type of atom in the quinoline ring (compound 7) may be of significance for each interaction. However, the contribution of these heteroatoms to the interaction strength with polar amino acids in the kinase pocket can particularly be determined by conformational factors. The stereochemical factors influence also contacts and interaction energy of the dimeric heterocycles containing an additional pyrrole, pyrazole, or indole ring. Crucial differences in the geometry of docked conformations can be expected for the fused azoles 6 and 7 as well.

Fitting of the first poses of azoles 1–7 in the Chk1 domain resulted in the formation of several hydrogen bonds between the ligands and the kinase amino acids (Fig. 2; Table S1, Figures S8–S15, Supplementary material). Apart from agents 5 and 6, all other indazole derivatives form an H-bond between the pyrrolic nitrogen atoms and the carbonyl of the backbone Glu85. These results generally agree with those obtained in the previous studies, e.g. the interactions of carbazole NH atom with Glu85 [40] albeit the NH…O=CGlu85 contact for indazole derivatives seems to be shorter and stronger than those reported for dihydroindeno[2,3-c]pyrazoles, i.e. 2.03–2.267 vs 2.7 Å [20].

Fig. 2

Docking poses of azoles 1–7 (1st poses)

Let us now consider the role of the sulphone group for the ligand–Chk1 interactions. The sulphone or sulphonamide functionality is frequently present in the anticancer molecules where they act as hydrogen bond acceptors from amino acids with pendent amino or hydroxy groups often in the kinases domains [41-43]. Moreover, the sulphone functionality often works as an important spacer for kinase inhibitors [44].

Our studies have shown that the sulphone group in all azole derivatives, except indole 5, form a hydrogen bond with the hydroxy group of Ser147, although for compound 6, this polar interaction is of a long-range nature.

Another hydrogen bond is established between nitrogen N2 of the pyrazole ring in compound 6 and the hydroxy group of Ser147. The HB distance was observed to be shorter than 3 Å, thus within the limits for fairly strong hydrogen bonding. Such values agree with the previously reported HB for the kinase ligands containing pyrazole system [45, 46].

The major hydrogen bonding for ligand 5 is formed between the indazole nitrogen N2 and hydroxy group of Tyr20, whereas the interaction involving the sulphonyl group as an acceptor and hydrogen donating ammonium group of Lys38 is considerably weaker. This docking result agrees with the reported observation concerning interactions of Tyr20 or Lys38 with heterocyclic nitrogen atoms [45, 47, 48]. We also found that the hydrogen bond between the pyrrolic nitrogen NH and carbonyl group of Asn135 significantly stabilized ligand 6 in the kinase domain. Additional support for this finding comes from the literature data [4].

The HBs and hydrophobic interactions of docked azoles 1–7 with the kinase Chk1 pocket amino acids are given in Figures S16–S22 (Supplementary material). Moreover, π–π stacking interactions involving Tyr89 and ligands 1–4 and 7 as well as π–cation interactions involving Tyr86 and ligands 1–7 were observed for the best poses.

The above discussion leads to a conclusion that the 3-tosylindazole ligands with chlorine or small molecular azole substituent on position 5 of indazole form similar overlapping poses in the kinase Chk1 pocket. On the other hand, the presence of a larger substituent like indole (compound 5), or an entirely different heterocyclic system like the pyrazolopyrazole condensed ring (compound 6), results in a different pose and dissimilar distribution of hydrogen bonding network in the enzyme cavity.

Optimization of ligands and residues involved in hydrogen bond formation

To further explore the geometry of docked ligands 1–7 (the best poses) in the Chk1-binding pocket and their interactions with the pocket amino acids, we decide to apply the DFT formalism and the B3LYP/6-31G(d,p) approach. To this end, we optimized only the ligands (Table 1; Cartesian coordinates in Supplementary material) or, additionally, the amino acid functionalities that form hydrogen contacts with the docked ligand (Table 2), while the remaining amino acids were treated as frozen. The computations were limited to a 4-Å sphere around the ligand. Furthermore, the C- and N-terminal amino acids were treated as close-shell systems, i.e. the carboxyl and amino groups were unionized. Analogously to the docking procedure, the protein coordinates were frozen for the calculations of energy changes.

Table 1

Ligand-amino acid distances (Å) between azoles 1–7 and amino acid residues within a 4-Å sphere around the ligand after ligand optimization with B3LYP/6-31G(d,p) method; RMSDcomplex = 0.351 (1), 0.382 (2), 0.505 (3), 0.485 (4), 0.470 (5), 0.743 (6), 0.893 (7) Å, respectively; the corresponding contact distances for the ligands docked to 2e9n.pdbqt were given in brackets

Contacts	Distances calculated for optimized azoles 1–7
	1	2	3	4	5	6	7
N-H…O=CGlu85	2.841	2.854	3.077	2.826	☓	☓	2.630
SO2…H-OSer147	2.517	2.500	3.266	2.817	☓	☓	2.961
N2indol…O-HTyr20	☓	☓	☓	☓	3.063	☓	☓
N-H…OOCAsp148	☓	☓	☓	☓	3.787 (7.115)	☓	☓
N-H…OSer147	☓	☓	☓	☓	☓	3.514	☓
N-H…OGlu17	☓	☓	☓	☓	3.787 (3.309)	☓	☓

Table 2

Ligand-amino acid distances (Å) between azoles 1–7 and amino acid residues within a 4-Å sphere around the ligand after ligand optimization with B3LYP/6-31G(d,p) method; RMSDcomplex = 0.337 (1), 0.368 (2), 0.423 (3), 0.484 (4), 0.337 (5), 0.748 (6), 1.070 (7) Å; the contacts resulted from the docking procedure are given in brackets

Contacts	Distances calculated for optimized azoles 1–7
	1	2	3	4	5	6	7
N-H…O=CGlu85	2.866	2.926	2.661	2.871	☓	☓	2.654
SO2…H-OSer147	2.644	2.593	3.377	2.791	☓	☓	2.846
SO2…H-NLys38	☓	☓	☓	☓	3.060	☓	☓
N2indol…O-HTyr20	☓	☓	☓	☓	2.402	☓	☓
N-H…OOCAsp148	☓	☓	☓	☓	3.330 (7.115)	☓	☓
N-H…OSer147	☓	☓	☓	☓	☓	3.503	☓
N-H…OGlu17	☓	☓	☓	☓	3.330 (3.309)	☓	☓

The HB distances NH…O=CGlu85 for ligands 1–7 were elongated in the geometry optimization by DFT approach in relation to the docking results. For these contacts, the elongation was the most noticeable for pyrazole 3 (Δ = 0.81 Å) and the least significant for quinoline 7 (Δ = 0.57 Å), whereas for azoles 1, 2, and 4, this lengthening was of a similar value (0.75–0.79 Å).

The H-bond distance SO2…HOSer147 was significantly elongated for pyrazole derivative 3 (Δ = 0.68 Å) and the least considerably lengthened for chlorine derivative 1 (Δ = 0.068 Å) while for the remaining azoles, this distance was elongated of 0.251 (2), 0.234 (4), or 0.374 Å (7).

The above data suggest that the stereoelectronic properties of Glu85 that influence energy and geometry changes of the optimized ligands were the most significant for azoles 1–4 and 7. This effect involving Ser147 was negligible for azole 1. The oxygen atom of the Ser 147 hydroxy group was closer to the hydrogen atom of pyrrolic NH in azole 6, but still, it was a weak long-range interaction. The bond OH…N-pyridinic atom (ligand 6) underwent a significant elongation during the refinement; thus, we regarded this bond as a negligible contact. Although the NH…OOCAsp148 contact involving indazole 5 was shortened, it was still within the long-range interactions. The polar interaction of pyrrolic NH…OGlu17 involving the same ligand 5 was considerably elongated from 3.309 to 3.787 Å in comparison with the docking results. If we consider in the optimization the residual functionalities of the amino acids involved in HBs together with the ligand geometry, this would affect the contacts arrangement (Table 2) in the kinase pocket. Similarly to the previous model, all contacts were elongated.

For the NH…O=CGlu85 interaction, the largest change was observed for 3,5-dimethylpyrazole derivative 4 (Δ = 0.840 Å), whereas the smallest alteration was noticed for pyrazole derivative 3 (Δ = 0.394 Å). For the remaining azoles 1, 2, and 7, this contact was stretched out of 0.784, 0.677, and 0.594 Å, respectively.

The strongest elongation effect of the SO2…HOSer147 contact was observed for pyrazole 3 (Δ = 0.791 Å), while the weakest effect was detected for chlorine derivative 1 (Δ = 0.059 Å), thus similarly to the previous interaction model. This HB was lengthened of 0.344, 0.208, or 0.259 Å for azoles 2, 4, or 7, respectively.

Analogously to the above model, the stereoelectronic impact of Glu85 upon the geometry changes of the optimized ligand was significant for azoles 1–4 and 7. However, this effect concerning Ser147 was negligible for azole 1.

The distance between the oxygen atom of Ser147 and the pyrrolic nitrogen atom of condensed azole 6 (3.503 Å) was within the limits of weak interactions. The contacts involving the pyridinic nitrogen of azole 6 and Ser147 hydroxy group as well as the Asp148 carboxyl and indazole 5 pyridinic nitrogen were practically imperceptible. However, the hydrogen bonding including the pyrrolic nitrogen in ligand 5 and carboxyl oxygen of Glu17 was almost unchanged in comparison with the docking result: 3.330 Å vs 3.309 Å, respectively (Table 2).

Next, we applied the semiempirical method PM7 for verification of the interactions between azole ligands 1–7 and the amino acids present in the kinase 2e9n.pdb pocket. The optimization was carried out for the functional amino acid groups that could form contacts with the best poses of ligands 1–7 within a 4-Å sphere leaving other coordinates of the kinase frozen. Then, we inspected distribution of the ligand–amino acid contacts in the optimized complexes (Table 3) and compared them with the previously obtained docking data (Table S1, Supplementary material). The ligand geometry in such optimized complexes did not differ from the corresponding best poses of azoles 1–7 obtained after docking procedure to the Chk1 pocket. The root mean square deviation (RMSD) value was 0.13–0.16 Å for the protein–ligand complexes 1–4, 6, and 7, while for indole derivative 5, it equalled 0.30 Å. In comparison with the docking results, the contacts involving pyrrolic atom NH in ligands 2–4 were virtually unchanged, whereas for chlorine and quinoline derivatives 1 and 7, these bondings considerably weakened and were beyond the strong interactions (d ≥ 2.2 Å). This seems to indicate that the NH…O=CGlu85 contact is important for the structure of complexes involving ligands 2–4 but less significant for the other complexes. On the other hand, the interaction between the Ser147 hydroxy group and azoles sulphonyl function seems to be unimportant because it disappears upon computations. The contact between the Ser147 hydroxy group and pyridinic nitrogen atom N2 of the pyrazole moiety was elongated significantly for condensed azole 6, but the contact NH…O=CAsn135 practically did not reshape in comparison with the docking results.

Table 3

Ligand-amino acid distances (Å) between azoles 1–7 and residues within a 4-Å sphere around the ligand after ligand optimization with PM7 method; RMSDcomplex = 0.130 (1), 0.127 (2), 0.138 (3), 0.119 (4), 0.292 (5), 0.143 (6), 0.157 (7) Å, respectively; the contacts resulted from the docking procedure are given in brackets

Contacts	Distances calculated for optimized azoles 1–7
	1	2	3	4	5	6	7
N-H…O=CGlu85	2.432	1.974	2.060	2.047	☓	☓	2.493
SO2…H-OSer147	5.774	5.409	5.371	5.089	☓	☓	5.386
SO2…H-NLys38	2.947 (4.037)	3.478 (4.054)	3.352 (3.937)	3.492 (4.053)	1.699	1.591 (4.303)	2.953 (3.970)
N2indol…H-NLys38	☓	☓	☓	☓	3.157 (5.320)	☓	☓
N2indol…O-HTyr20	☓	☓	☓	☓	5.624	☓	☓
N-H…OOCAsp148	☓	☓	☓	☓	1.861 (7.115)	☓	☓
N-H…O=CAsn135	☓	☓	☓	☓	☓	2.584	☓
N2…H-OSer147	☓	☓	☓	☓	☓	3.556	☓

The sulphonyl group on azoles 1–4 and 6–7 participated in a long-range polar interaction with the ammonium group of Lys38 or even could form a hydrogen bond (azole 6) with this function as opposed to the docking results. A complete loss of interaction was observed for the pyrazole pyridinic nitrogen in ligand 5 and Tyr20 hydroxy group. However, a new hydrogen bond was spotted between the pyrazole pyridinic nitrogen and Lys38 ammonium group. The above described contact involving indazole pyrrolic nitrogen (indazole 5) and Glu17 carboxyl disappeared.

Estimation of interaction energy

In the next step, we focused on the assessment of enthalpy changes of the interactions of azole ligands 1–7 (ΔHint) in the Chk1 pocket. In this evaluation, we considered values of HOF under standard conditions using Mopac 2016 program and its implemented module, Mozyme [26]. To study the interactions between ligand and kinase pocket, the binding sphere was limited to 4 Å from the best pose. The pocket amino acids were correctly protonated, and the C- and N-terminal amino acids were ionized to obtain COO− or NH3+. Then, the hydrogen atoms of ligand–protein complex were optimized as well as the ligand environment leaving the COO− or NH3+ groups frozen. The resulted hydrogen bonds and polar interactions in such optimized complexes are shown in Table 4.

Table 4

Ligand-amino acid distances (including hydrogen bonds) (Å) between azoles 1–7 and residues within a 4-Å sphere around the ligand after ligand optimization with PM7 method; RMSDcomplex = 0.604 (1), 1.085 (2), 0.622 (3), 1.438 (4), 0.790 (5), 0.058 (6), 0.689 (7) Å, respectively; the contacts resulted from the docking procedure are given in brackets

Contacts	Contacts length calculated for optimized azoles 1–7
	1	2	3	4	5	6	7
N-H…O=CGlu85	1.841	1.913	1.855	1.793	☓	☓	3.244
SO2…H-OSer147	2.919	4.406	3.257	3.754	☓	☓	1.792
SO2…H-NLys38	☓	☓	☓	☓	2.475	☓	4.155 (3.970)
N2indol…H-NLys38	☓	☓	☓	☓	4.289 (5.320)	☓	☓
N2indol…O-HTyr20	☓	☓	☓	☓	3.758	☓	☓
N-H…OO=CAsp148	☓	☓	☓	☓	4.870 (7.115)	☓	☓
N-H…O=CAsn135	☓	☓	☓	☓	☓	1.750	☓
N2…H-OSer147	☓	☓	☓	☓	☓	4.565	☓
N-H…OGlu17	☓	☓	☓	☓	3.486 (3.309)	☓	☓

For the interaction energy calculations, we adopted an approach based on the thermodynamic cycle of Raha and Merz [34]:ΔHint = ΔHf(PL) − [ΔHfcomplex(P) + ΔHfcomplex(L)], where ΔHf(X) were the heats of formation in vacuo of the protein–ligand complex, free ligand (L) or free protein (P), and the ΔHfcomplex(X) parameter corresponds to the enthalpy of the protein or ligand molecule in the complex conformation. An important note must be made here: this procedure does not account for such processes as desolvation of the binding pocket. However, we compare energies of the structures docked into the same binding pocket; therefore, the desolvation energy should be comparable among the 1–7 series.

The application of the above equation to the complexes of ligands 1–7 with Chk1 led to the values shown in Tables 4 and 5. These values provide evidence that 3,5-dimethylpyrazole derivative 4 and indole derivative 5 are the ligands best fitted in the Chk1 pocket. For compound 4, this conclusion is in accordance with the docking result (Table S1, Supplementary material), whereas for compound 5, it seems that the estimation of binding affinity by molecular modelling is rather inconclusive. Moreover, this result is consistent with the above discussed models indicating that ligand 4 forms indeed the strongest interactions with Glu85 and Ser147 in the Chk1 pocket. Noteworthy is the fact that the NH…O=CGlu85 contact is shortened for ligands 1–4 but is elongated for quinoline derivative 7 in comparison with the docking results. The ligand 1 SO2…HOSer147 distance has almost similar value, but the same contact was significantly lengthened for ligands 2–4 or shortened for ligand 7. The interaction energy for azoles for azoles 2, 3, and 6 was on the similar level, whereas azole 1 with chlorine substituent was energetically disfavoured (Table 5). Subsequently, the SO2…HNLys38 contact was virtually unchanged for indazole 5, whereas the pyridinic N…HOTyr20 interaction was markedly stretched out and weakened. The interactions concerning the pyridinic nitrogen atom of indazole 5 with Lys 38 and Asp148, described in the previous model involving PM7 optimization of the best docking pose, practically disappeared. The ligand 5 pyrrolic NH…O=CGlu17 contact (3.486 Å; Table 4) was comparable to the docking results (3.309 Å; Table S1, Supplementary material) but Glu17 was in the periphery of the kinase pocket. We can conclude that the share of Tyr20 in the relatively high interaction energy (ΔHint = − 111.88 kcal mol−1; Table 5) is marginal. The PM7 results concerning azoles 5 and 7 differed from the estimated binding affinity obtained in the docking protocol.

Table 5

Calculated heats of formations (kcal mol−1) for free ligands (ΔHfcomplex(L)), free protein (ΔHfcomplex(P)), ligand–protein complex (ΔHf(PL)), as well as ligand–protein interaction energy (ΔHint)

Compound	HOF of ligand (ΔHfcomplex(L))	HOF of protein (ΔHfcomplex(P))	HOF of complex (ΔHf(PL))	ΔHint
1	− 1.11	− 1164.89	− 1251.06	− 85.06
2	43.41	− 2046.66	− 2105.93	− 102.67
3	62.88	− 1370.37	− 1405.87	− 98.39
4	46.76	− 2130.38	− 2207.67	− 124.04
5	55.48	− 1844.07	− 1900.47	− 111.88
6	49.99	− 1673.68	− 1723.83	− 100.14
7	34.23	− 1212.30	− 1276.91	− 92.18

Molecular dynamics calculations

Next, we applied molecular dynamics (MD) to explore stability of ligands 1–7 in the Chk1 pocket. For this purpose, the GROMACS 2016.4 [35, 36] was employed to simulate the solvated ligand–protein complexes (cubic water boxes). The RMSDs of backbone atoms (with or without the presence of ligands 1–7) were calculated with respect to the initial configuration. The time evolution of RMSD values of the backbone in the ligand–protein complexes is shown in Fig. 3.

Fig. 3

The RMSD plot for the backbone within ligand–protein complex during the productive phase calculated for free kinase (black), as well as its complex with the following: 1 (red), 2 (green), 3 (blue), 4 (yellow), 5 (grey), 6 (light grey), and 7 (purple)

The ligand RMSD plot shows that the docking poses of all azole ligands, apart from chlorine derivative 1, are stable inside the kinase pocket, especially if we compare this plot with the RMSD simulation relative to the empty kinase pocket. Generally, ligands 2–7 remained the backbone stable (regarding their influence on protein) for ca 5 ns in their positions with RMSD oscillating at 1.50–2.00 Å; then, the RMSD curve moved upward and oscillated at 2.00–2.50 Å. In most cases, the RMSD values were smaller than 3.0 Å, which clearly indicates that these systems were quite stable during the last 12 ns. Examining the RMSD plot, we can conclude that apparently ligands 2 and 5 make the backbone most stable.

Hydrogen bonds are important interactions within an enzyme active site, especially these formed when an enzyme substrate reaches the transition state [49]. The phosphorylation process, characteristic for kinases, does not require strong HBs between the enzyme substrate and amino acids present in the active site.

However, HBs and short electrostatic contacts cannot be distinguished reliably by the RDF analysis. As this study was intended to initiate further computational research by us in the field of azole ligand affinity, we chose to proceed with a more detailed analysis of the hydrogen bonds along the MD trajectory.

The geometric criteria proposed for HBs indicate that the donor–acceptor distance D should be less than 3.5 Å and the Ѳ angle more than 130° [50]. However, these criteria for proteins are less restricted (D < 3.9 Å and Ѳ >  90° [51]) due to the large structural flexibility of the enzyme active site during the catalysis process. These structural alterations cause fluctuation of the hydrogen bonds array. Considering the dynamic nature of the contacts formed, we finally chose to use rather conservative values to define the occurrence of hydrogen bonding: D < 3.3 Å and Ѳ > 135°. Therefore, we focused on the occupancy parameter relative to HBs formed by azoles 1–7 with the above-discussed amino acids.

The occupancy parameter involving Glu85 has the following values: 12.66 (1), 23.72 (2), 14.58 (3), 0.11 (4), and 0.26% (7). We can conclude that the impact of Glu85 on the possible formation of hydrogen contacts is significant for azoles 1–3. This matches with the previously discussed models (Tables 1, 2, 3, and 4). On the other hand, we did not observe significant interactions between Ser147 and azoles 1–3, or the occupancy parameter had a small value of 0.03 or 0.04% for azoles 4 or 7. This relatively minor role of Ser147 in the formation of HBs is confirmed by the previously obtained data shown in Tables 1, 2, 3, and 4. A similar value of the occupancy parameter (0.004%) was obtained for azole 5 and Tyr20. For azole 5 and Glu17, however, this parameter reached 9.79%. Therefore, Glu17 participates considerably in the electrostatic interaction with the NH atom of azole 5 that is in accordance with the data presented in Tables 1, 2, and 4.

For azole 6 and Ser147 or Asn135, we did not observe the occupancy parameter relative to the formation of HBs.

The half-life time for the hydrogen contacts involving azoles 1–4 and 7 and Glu85 was the longest in the simulation (Figure S23 in the Supplementary material). However, these contacts half-life time was shortened considerably in relation to the Ser147 interactions (Figure S24 in the Supplementary material). Therefore, Ser147 influence on the formation of HBs with azoles 4 and 7 seems to be less critical, which corresponds with the occupancy parameters and conclusions drawn from the previous models.

The corresponding half-life time involving azole 5 and Tyr20 was relatively short and did not exceed 1.45 ps (Figure S25 in the Supplementary material). Consequently, Tyr20 impact on the formation of HB with the azole 7 pyridinic nitrogen was insignificant.

In contrast to the previous case, the half-life time involving Glu17 and azole 5 was markedly longer at 21 ps (Figure S25 in the Supplementary material). Therefore, this interaction should be regarded as necessary for the ligand stabilization in the kinase pocket. Both the above conclusions concerning azole 5 follow the occupancy parameter and conclusions drawn from the previous models. A relatively high value of interaction enthalpy (Table 5) may be connected with the additional π–cation contact with Lys38.

The MD simulation for condensed pyrazole 6 did not show the formation of HBs with Ser147 and Asn135 that would fulfil the adopted cut-off criteria. This observation agrees with the occupancy parameters and data obtained from the DFT formalism and PM7 method.

Conclusions

Checkpoint kinase 1 is a pivotal element of the checkpoint signalling involving DNA damage response. As the DDR occurs in reaction to DNA destruction by cytotoxic agents and radiation, Chk1 is a potential target for kinase inhibitors [52]. Our previous investigations showed that some indazole derivatives with a toluene sulphonyl substituent at the position 3 of indazole displayed anticancer effect against colon cancer cell line HT29 as well as breast cancer cell lines MCF7 and MDA-MB-231. In the present paper, we provided evidences that such derivatives of indazole and condensed pyrazole can represent a valuable template hit for the future Chk1 inhibitors. The molecular docking to the kinase ATP pocket has shown that the refined azole conformers 1–4 and 7 access the bonding region in Chk1 and form relatively strong hydrogen bonds with the key amino acids, i.e. Glu85 and Ser147. The docking has shown that that the 3-arylsulphonylindazole ligands with chlorine or small molecular azole substituent on position 5 of indazole form similar overlapping poses in the kinase Chk1 pocket. On the other hand, the presence of a larger substituent like indole (compound 5), or an entirely different heterocyclic system like the pyrazolopyrazole condensed ring (compound 6), results in a different pose and dissimilar distribution of hydrogen bondings in the enzyme cavity.

The above observation has been confirmed by the DFT and semiempirical PM7 methods although these computations have given somewhat inconsistent results for the interactions between the ligands and Ser147, particularly for compounds 1, 5, and 6. The calculated interaction energy favours azoles 2, 3, 5, and 6, whereas the remaining azoles, particularly compound 1, are relatively energetically disfavoured.

The molecular dynamics simulations indicate that almost all complexes ligand–kinase Chk1, apart from compound 1, presented stable MD trajectories. We can conclude that the impact of Glu85 on the possible formation of hydrogen contacts is significant especially for azoles 1–3. This matches the previously discussed models (Tables 1, 2, 3, and 4). However, the occupancy parameters for the interactions of Ser147, Tyr20, and Asn135 with ligands 1–7 are usually of small values (apart from azole 5 and Glu17). This may indicate that formation of hydrogen bonds involving the above amino acids is less probable.

(DOC 41358 kb)