DFT investigations on arylsulphonyl pyrazole derivatives as potential ligands of selected kinases

4.1 Estimation of atomic charges

Indazole derivatives 1–6 and condensed pyrazoles 7 and 8 (Scheme 1) were initially optimised by molecular mechanics, then a series of 1,000 conformers for each azole derivative were further refined by the semiempirical method PM7. Next, for each of the azole derivatives 1–8, a conformer with the smallest heat of formation was selected and optimised in the SCF (self-consistent field) procedure (density functional theory [DFT] method, B3LYP functional, and Gaussian 16 A.03 program [28]).

To evaluate the interactions of azoles 1–8 with perturbing external field, we carried out an analysis of the electron density distribution among the azole atoms. First, we decided to investigate MEP corresponding to the above charge distribution. The MEP was determined by the B3LYP/6-311++G(2d,3p) approach for the conformers of azoles 1–8 with geometry previously optimised at the B3LYP/6-31G(d,p) level of theory in the gaseous phase using tight criteria of convergence (Gaussian 16 A.03 program [28], Figure 1). We determined both Mulliken and natural bond orbital charges (NBO 3.0 as implemented in the Gaussian 16 A.03 program). We also used the adjusted charge partitioning (ACP) approach [29] as well as a newly developed more accurate scheme, namely the iterative-adjusted charge partitioning (I-ACP) method [30]. In our investigations, involving the multilevel approach to the conformational rotamer search, the results were refined using a basis set enriched with the higher-level polarisation functions. The results are shown in Table 1 and Tables S1–S8 (Supplementary material, corresponding atom numbering is given in Figures S1–S8). The results in Table 1 apply only to the pyrazole pyrrolic nitrogen atoms except for compound 8, where the endocyclic quinoline nitrogen is considered as well.

Figure 1

Electrostatic potential map of investigated azoles 1–8 (B3LYP/6-311++G(2d,3p)//B3LYP-631G(d,p); gaseous phase); isovalue = 0.002 a.u.

Depending on the approach relative to the aforementioned nitrogen atoms, different charge distributions can be observed. The Mulliken population analysis showed that the charges of the pyrrolic nitrogen atoms were arranged from the most to the least electronegative in the following order: 7 < 1 < 6 < 4 < 5 < 2 < 3 < 8. On the other hand, the NBO scheme provided another sequence of the pyrrolic nitrogen atoms, i.e. 7 < 8 < 3 = 4 < 6 < 5 < 1 < 2. Thus, both Mulliken and NBO methods indicated the pyrrolic nitrogens of pyrazolopyrazole 7 as the most electronegative atoms among the analysed heteroaromatic rings.

The charges of the pyrrolic nitrogen atoms of azoles 1–8 determined by the ACP or I-ACP methods did not differ significantly and equalled from −0.248 to −0.246 or from −0.319 to −0.316, respectively. The pyrrolic nitrogen located within the structure of fused azole 7 was indicated again as the most electronegative atom. Note that for compound 8, the most electronegative atom was the quinoline pyridinic nitrogen, regardless of the computation method. The most substantial difference between the pyridinic and pyrrolic nitrogen for this azoloazine was observed in the Mulliken approach.

Despite the differences in the calculated charges for the pyrrolic and pyridinic nitrogen atoms, the obtained results lead to a conclusion that the pyrrolic nitrogen, an H-bond donor, and the electron-withdrawing tosyl substituent are the most important for the azole–protein interactions. Moreover, the pyridinic nitrogen present in the quinoline ring can also participate in such interactions as an H-bond acceptor.

4.2 Estimation of condensed Fukui functions

Reactivity indicators are essential descriptors of the responsiveness of a substance to a perturbing external field. The Fukui function originated from the conceptual DFT [31] and relates to the local electron density changes in a molecule during a reaction. It is known that the Fukui function, f(r), is attributed to the distortion of the electron density in an open system, while the external potential remains constant. Due to the discontinuity of electron density relative to a given number of electrons in the system, in the literature three types of the Fukui functions are defined, i.e. a backward Fukui function f⁺(r), forward Fukui function f⁻(r), or average Fukui function f⁰(r) corresponding to a nucleophilic, electrophilic, or radical attack, respectively.

In our investigations, the electronic densities and charges for the neutral species were evaluated at the B3LYP/6-311++G(2d,3p) level of theory using the geometries of the previously optimised neutral species in the gaseous phase based on the time-dependent DFT analysis (for details regarding HOMO–LUMO energy calculations in terms of global descriptors of reactivity see Section 5) [32] as a quick protocol for the estimation of the susceptibility to different kinds of attack. For computations of condensed Fukui functions, we used software implementation that corresponded to the literature approach [33]. The results are listed in Table 2 and Tables S9–S15 (Supplementary material).

According to the HSAB principle, the functionalities of an analyte with the large value of Fukui function are softer than the regions where the Fukui functions are small. Our computations showed that (f⁺(r)), describing interaction with a nucleophile, had similar values for the pyrrolic nitrogen atoms in azole rings 1 and 2 (Table 2). Analogously, the (f⁺(r)) functions for these atoms in azoles 4 and 5 were almost identical. Considering the values of (f⁺(r)), azoles 1–8 can be arranged in the following ascending order of nucleophilicity: 8 < 3 < 6 < 4 < 5 < 1 < 2 < 7. This order may reflect the susceptibility of these compounds to possible interactions with polar groups of the amino acids present in binding pockets of protein kinases.

The reactivity of the pyrrolic nitrogen atoms in azoles 1–8 towards electrophiles (f⁻(r)) or radicals (f⁰(r)) increased in a similar order, i.e. 2 < 6 < 5 < 4 < 8 < 3 < 1 < 7 or 2 < 6 < 5 < 8 < 4 < 3 < 1 < 7, respectively. This sequence was different only for azoles 4 and 8. For all condensed Fukui functions, the highest values could be attributed to the pyrrolic atom N–H of the fused pyrazolopyrazole 7. For quinoline 8, the Fukui function (f⁻(r)) had the highest concentration on the nitrogen atom within the pyridine ring. This nitrogen atom is the centre of the electrophilic attack. The lowest value of the function f⁻(r) was observed for chlorine derivative 1. The above observations may have profound implications concerning the formation of hydrogen bonding between the azole ligand and amino acids present in kinase pockets.

4.3 Estimation of dipole moment and results of NBO calculations

The molecular dipole moment is an intrinsic property of a molecule and affects its solubility in solvents as well as drug–receptor interactions. The interactions involving the dipole moment include permanent dipole–ion and permanent dipole–permanent dipole interactions, as well as van der Waals or London dispersion forces. The first type of interactions arises when an ionic molecule approaches a compound with a permanent dipole, while the second one appears between molecules with fixed dipoles. The van der Waals forces are very weak interactions and involve attractions between molecules with temporary induced dipoles.

Our computations of molecular dipoles of azoles 1–8 involving the B3LYP/6-311++G(2d,3p)//B3LYP-631G(d,p) (gaseous phase) approach together with ACP and I-ACP implementations led to a conclusion that the dipole moments increased in the following order (Table 3): 4 < 7 < 8 < 3 < 6 < 2 < 5 < 1. The orientation of the dipole moments for azoles 1 and 4 is shown in Figure 2 (Figures S1–S8 in the Supplementary material). Note that the values of these dipole moments correspond with the literature data for analogous sulphones [34,35].

Figure 2

Visualisation of the dipole moment of azoles 1 and 4 (B3LYP/6-311++G(2d,3p)//B3LYP-631G(d,p); gaseous phase). The atoms forming the dihedral angles are marked in light grey colour; the molecular dipole moment is shown with the blue arrow.

The influence of the dihedral angle between indazole and tosyl planes on the molecular symmetry seemed to be negligible because its value was comparable for all azoles 1–8. On the other hand, the dihedral angle between the azole substituent at position 5 of indazole and the indazole ring had significantly different values for the dimeric type of azoles 2–6, i.e.: −35.58 (2), −16.42 (3), 38.36 (4), −45.02 (5), and −55.43° (6) (Figures S2–S6 in the Supplementary material).

The lack of steric strain led to an increase in conformational freedom for compound 1. The presence of fused azaaromatic systems within azoles 7 and 8 significantly affects their dipole moments. The above conclusion is connected with the hybridisation of the nitrogen atoms linked to position 5 of the indazole ring. This relationship can particularly be observed for compounds 2–6. The NBO method (B3LYP/6-311++G(2d,3p)//B3LYP-631G(d,p); gaseous phase) demonstrated that these nitrogen atoms had the largest contributions to the two-centre bonding orbital (σ_C–N) between the indazole carbon 5 and the pyrrolic nitrogen of azoles linked to this carbon. These contributions corresponded to the following values: 62.67 (2), 63.10 (3), 62.90 (4), 62.78 (5), and 62.82% (6). The N hybrids that participated in the σ_C–N bond formation had predominant p character, i.e. s (34.41%) p 1.90(65.52%) (2), s (36.14%) p 1.77(63.79%) (3), s (35.49%) p 1.82(64.44%) (4), s (34.71%) p 1.88 (65.23%) (5), and s (34.62%) p 1.89(65.32%) (6).

In the relationship NBO–donor, represented by the bonding orbital BD, and NBO–acceptor, represented by the antibonding orbital BD*, we considered the highest stabilisation energy E(2). For pyrrole 2, the bonds between pyrrolic nitrogen and neighbouring carbon atoms in the pyrrole ring were NBO acceptors with a stabilisation energy of 1.64 kcal/mol. We found that for pyrazole derivative 3, the N_pyridinic–C bond within the pyrazole moiety functioned as an NBO acceptor and had a significant contribution to the stabilisation energy (2.60 kcal/mol). On the other hand, for dimethylpyrazole derivative 4, the N_pyrrolic–C bond was an NBO acceptor with the stabilisation energy of 2.11 kcal/mol. For indazole 5, we observed that the bond between carbons 3a and 7a, forming the indazole rings junction, was an NBO acceptor with the stabilisation energy of 3.32 kcal/mol. For carbazole derivative 6, the stabilisation energy concerning the bond between the pyrrolic nitrogen and adjacent carbon of the pyrrole ring was 3.53 kcal/mol.

Using the NBO method, we analysed the N–H bonds of azoles 1–8 (Figures S1–S8, Supplementary material). The highest contribution to the σ_N–H bond was found for the nitrogen atoms. That contribution was as follows: 71.03 (1), 70.99 (2), 70.99 (3), 70.98 (4), 71.00 (5), 71.02 (6), 71.05 (7), and 70.98% (8). We concluded that the nitrogen atoms in this bond were sp^2.31 (azoles 1–5), sp^2.26 (7), or sp^2.33 (8) hybridised that translated into more than 69% p in character.

By analysing the interactions between NBO donor (σ_N–H) and NBO acceptor (antibonding BD* orbital), and considering the highest stabilisation energy E(2), we demonstrated that the bond between the indazole pyridinic nitrogen and neighbouring carbon was an NBO acceptor in azoles 1–8, and the energy E(2) had the following values: 2.11 (5), 2.13 (6), 2.17 (1–2, 4), 2.18 (3), 2.19 (8), or 2.29 (7) kcal/mol.

4.4 AIM analysis of N–H bonds

QTAIM scheme is a sound way for analysis of various interactions concerning neutral molecules, ions, complexes, molecular aggregates, or even crystals [26,36]. Applying the QTAIM method, we can analyse the distribution of the electron density between atoms in a molecule. Considering mobility of the protons linked to the pyrrolic nitrogen atom in azoles 1–8 as well as the conclusions involving the analysis of electron density distribution and NBO results, we decided to study the N–H bond of these compounds by the QTAIM method (B3LYP/6-311++G(2d,3p)//B3LYP-631G(d,p); gaseous phase). We focused mainly on electron density, i.e. ρ_BCP as a sum of the kinetic and potential density of electrons at the bond critical point (BCP), Laplacian of electron density (∇²ρ_BCP), and bond ellipticity (ε) at the BCP.

The results are listed in Table 4 (Poincare–Hopf relationship was satisfied). A negative value of Laplacian for the electron density at the BCP shows whether the local electron density is concentrated between neighbouring atoms’ characteristic of a covalent bond. The ellipticity is a parameter used to quantify the double bond character of chemical bonds; thus, it allows us to differentiate between single and double bonds.

Our investigations confirmed that the N–H bond had a cylindrical section corresponding to a single bond. Small differences could only be detected in the second digit after the decimal point for azoles 7 and 8. These differences were notably reflected in the most negative value of Laplacian of electron density for the N–H bond in the fused pyrazole derivative 7. For the remaining azoles, the differences in ∇²ρ_BCP were insignificant and were running at ca. –1.87 a.u. They were particularly small (0.35 a.u.) for the N–H bond. The lowest value of ρ was observed for BCP concerning the N–H bond in the chlorine derivative 1, whereas the highest for fused pyrazole 7.

4.5 Estimation of free enthalpy of solvation

The free energy is a key property of molecules that determines many chemical and biological processes. Accurate calculation of free energy has become the central objective for computer simulations of chem- and biorelevant systems [37]. By invoking the continuum approximation for the solvent, the electronic structure problem for a molecule in a liquid is reduced to the size of the solute of interest. Continuum solvation models (CSMs, implicit solvation models and implicit solvent models) represent a solvated molecule at an atomic level of detail inside a molecule-sized, and usually molecule-shaped, electrostatic cavity surrounded by a dielectric medium that represents the solvent. These methods have been parametrised to deliver accurate values of the free energies of solvation, which can be added to accurate values of the free energies in the gas phase, to obtain the corresponding solution-phase free energy [38]. An accurate computationally efficient model of aqueous solvation represents an important challenge to molecular modelling. CSMs such as the PCM [39], the conductor-like screening model [40], and the SMD [41] offer a computational efficient representation of solvation for molecules treated with electronic structure methods [27,41]. Basically, the intrinsic free energy of solvation corresponds to the change in Gibbs free energy of the solute in the gas and solution phases. It is noteworthy that CSM free energy of solvation is not calculated as the difference in solution- and gas-phase free energies. The actual computational protocol of most CSM models (PCM, solvation models, etc.) requires an equation as follows:

where E_s and E_g are the electronic energies of the solute in CSM and in the gas phase, and G_ne are the non-electrostatic contributions (cavitation, dispersion, and exchange–repulsion). Equation (1) assumes that any thermal contributions in the gas phase and in the solution should cancel and that any noncancelling differences are indirectly incorporated into (E_s + G_ne) through parametrisation. Thus, the free energy in solution in the CSM framework should be calculated by using the following equation:

where ΔG_g is treated as the thermal correction to the free energy estimated in the gas phase [41,42]. In our studies, we decided to estimate the free energy of solvation for investigated azoles 1–8 using the Truhlar SMD approach [41] and related Minnesota M062X functional [43]. For this purpose, the initial geometries of hetarenes 1–8 were optimised at the M062X/6-31G(d,p) level of theory (self-consistent reaction field using very tight criteria, conversion from gas-phase standard state to 1 M was 1.893 kcal/mol). The optimisation protocol was carried out both in vacuo and by applying the following 13 solvents (the value of their dielectric constant is given in bracket): n-hexane (1.8819), carbon tetrachloride (2.228), toluene (2.3741), chloroform (4.7113), chlorobenzene (5.6968), tetrahydrofuran (7.4257), n-octanol (9.8629), acetone (20.493), methanol (32.613), acetonitrile (35.688), N,N-dimethylformamide (37.219), dimethylsulfoxide (46.826), and water (78.3553). However, for the energy calculations of optimised ligands, we used the M062X/6-311++G(2d,3p) level of theory (in the gaseous phase, as well as using the SMD model and the above solvents). The results are listed in Table 5. The literature search has revealed a lack of data concerning the use of the SMD approach for the analysis of sulphonyl azole derivatives.

Our results were comparable in values with the energy determined by the SMD approach for some azaheterocyclic systems, such as caffeine or phthalimide, or to the energy calculated for toluene or compounds with an octahydro-1,2-benzisothiazole-1,1-dioxide core by this solvation model [44,45]. The results showed that for azoles 1, 2, 6, and 7, the free energy of solvation parameter has the most negative value in the reaction field of acetonitrile, whereas for azoles 3–5 and 8, these energy values were the lowest for acetone as a solvent. For all solvents, the most negative values of the free energy of solvation in the SMD model were attributed to carbazole 6, while the highest value of this energy could be ascribed to fused pyrazole 7. When we imitated cell environment using water as a solvent, the investigated azoles were ordered in the following way from the least to the most negative free energy of solvation: 6 < 5 < 4 < 3 < 8 < 2 < 1 < 7. The observed differences in that parameter could reflect the azole susceptibility towards interactions with amino acids present in the binding pockets of kinases as these targets exist in an aqueous environment. The above conclusion particularly applies to compounds 6 and 7.