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

Abstract Using the density functional theory (DFT) formalism, we have investigated the properties of some arylsulphonyl indazole derivatives that we studied previously for their biological activity and susceptibility to interactions of azoles. This study includes the following physicochemical properties of these derivatives: electronegativity and polarisability (Mulliken charges, adjusted charge partitioning, and iterative-adjusted charge partitioning approaches); free energy of solvation (solvation model based on density model and M062X functional); highest occupied molecular orbital (HOMO)–lowest occupied molecular orbital (LUMO) gap together with the corresponding condensed Fukui functions, time-dependent DFT along with the UV spectra simulations using B3LYP, CAM-B3LYP, MPW1PW91, and WB97XD functionals, as well as linear response polarisable continuum model; and estimation of global chemical reactivity descriptors, particularly the chemical hardness factor. The charges on pyrrolic and pyridinic nitrogen (the latter one in the quinolone ring of compound 8, as well as condensed Fukui functions) reveal a significant role of these atoms in potential interactions of azole ligand–protein binding pocket. The lowest negative value of free energy of solvation can be attributed to carbazole 6, whereas pyrazole 7 has the least negative value of this energy. Moreover, the HOMO–LUMO gap and chemical hardness show that carbazole 6 and indole 5 exist as soft molecules, while fused pyrazole 7 has hard character.


Introduction
Pyrazole and indazole derivatives can act as bioisosteric mimics for some important natural compounds, namely, pyrrole, indole, and purine. They display a broad spectrum of pharmacological activities and are present in many currently available drugs [1,2], including anticancer agents such as axitinib, crizotinib, pazopanib, ruxolitinib, ibrutinib, and niraparib [3,4]. Axitinib is an orally available multitargeted inhibitor of VEGFR, plateletderived growth factor (PDGFR), and c Kit (CD117) tyrosine kinase receptors. It has been approved for the treatment of renal cell carcinoma [5,6]. A pyrazole derivative, crizotinib, is a potent and selective dual inhibitor of anaplastic lymphoma (ALK) and c Met kinases. As the latter enzyme is a downstream effector involved in the pathogenesis of breast cancer, this drug is also in clinical trials for triplenegative breast cancer therapy [7,8]. Pazopanib is an indazolylpyrimidine derivative active orally. It is a potent and selective multitargeted inhibitor of VEGFR 1, VEGFR 2, VEGFR 3, PDGFR, and c Kit kinase that blocks tumour growth and angiogenesis [9]. Another pyrazole derivative, ruxolitinib, is a Janus kinase (JAK2i) inhibitor that was approved in 2011 (USA) for the treatment of intermediate and high-risk myelofibrosis and was in clinical trials for the use in breast cancer cure [10]. Ibrutinib was approved in 2014 (USA) for the treatment of B-cell cancers and is used in mantle cell lymphoma, leukaemia, and macroglobulinemia [11]. Niraparib, an indazole derivative, was approved in 2017 in USA and European Union for the maintenance treatment of adult patients with recurrent epithelial ovarian cancer. This drug is a novel poly(ADP-Ribose) polymerase inhibitor and is currently undergoing clinical trials for the treatment of triplenegative breast and prostate cancers [12,13].
The sulphone or sulphonamide moiety is also frequently present in the anticancer molecules where they operate as hydrogen bond acceptors from amino acids with pendant amino groups, e.g. lysine or asparagine, often in the affinity region of kinases [14][15][16][17]. Moreover, sulphone functionality often works as an important spacer for kinase inhibitors [18].
Our preliminary investigations showed that some indazole derivatives with the toluene sulphonyl (tosyl) substituent at position 3 displayed anticancer effect against colon cancer cell line HT29 as well as breast cancer cell lines MCF7 and MDA-MB-231. We suggested then that their mechanism of activity might involve interaction with magnesium ions, stimulation of proapoptotic proteins, or inhibition of signalling pathway proteins [19][20][21][22].
The purpose of the present study is to examine the electronic properties of azoles 1-8 (Scheme 1) with emphasis on the distribution of charges localised on the selected nitrogen atoms. We also aim to evaluate their condensed Fukui functions as well as global chemical reactivity descriptors (GCRDs). One key question relative to the computations was the estimation of the highest occupied molecular orbital (HOMO)-highest occupied molecular orbital (LUMO) gap of the aforementioned azoles. Such an approach involving hard-soft acid-base (HSAB) principles has been much discussed in the recent literature [23,24]. Our results are supported by analysis of dipole moments of the N-H bonds [25], quantum theory of atoms in molecules (QTAIM) [26], and calculations of free energy of solvation [27]. We believe that such broad and multidimensional comparative analysis concerning sulphonyl derivatives of azoles has not been undertaken before.

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.
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.

Azole
Type of nitrogen atom  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.

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 0 (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 0 (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.

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) 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].
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 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 twocentre bonding orbital (σ C-N ) between the indazole carbon 5 and the pyrrolic nitrogen of azoles linked to this carbon. 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 Table 4: QTAIM analysis of N-H bonds within the structure of azoles 1-8 (B3LYP/6-311++G(2d,3p)//B3LYP-631G(d,p); gaseous phase); electron density at the BCP (ρ BCP /a.u.), Laplacian of electron density (∇ 2 ρ BCP /a.u.), bond ellipticity (ε), the Hamiltonian form of kinetic energy density (K) 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)

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 (∇ 2 ρ 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 ∇ 2 ρ 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 ρ

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 moleculeshaped, 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 nonelectrostatic 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   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,1dioxide 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.
For hetarenes 1-6, the HOMOs were located beyond the tosyl substituent, and for compound 2, the HOMO was placed mainly on the pyrazole ring. On the other hand, the LUMOs had the highest coefficients for the indazole and tosyl rings. The conjugation present in azoles 7 and 8 affected the electron distribution and resulted in the placement of the HOMOs on the fused heterocyclic rings and sulphone group with low electron population on the peripheral toluene ring. The LUMOs of compounds 7 and 8 were different: for the fused pyrazole derivative 7, the LUMO covered almost the whole molecule, whereas for quinoline 8 practically only the quinoline moiety. The details concerning the TD-DFT analysis, including the HOMO and LUMO energies, HOMO-LUMO gaps, λ 1 related to the first computed excited state, and λ max are listed in Tables 6-9 (Tables  S25-S75 in the Supplementary material).
Next, we studied the vertical excited states of azoles 1-8 in the context of solvatochromic effect considering the B3LYP, CAM-B3LYP, MPW1PW91, and WB97XD functionals in the reaction field of 13 solvents. The details concerning the TD-DFT analysis including energy characteristics of the HOMOs and LUMOs, HOMO-LUMO gaps, λ 1 related to first computed excited state and   Tables 6-9). We analysed the estimated maxima of absorption λ max and arranged them according to the lowest and highest values for each of the ligand. Using the B3LYP functional, we obtained the lowest λ max values for azoles 1, 4, 7, and 8 in n-hexane, azoles 2, 5, and 6 in methanol, and azole 3 in water. The highest λ max values were observed for azoles 1 and 8 in DMF, azole 2 in toluene, azole 3 in n-octanol, azoles 4 and 6 in chlorobenzene, and azole 7 in DMSO. Using the CAM-B3LYP functional, the lowest values of λ max were found for azoles 1 and 8 in n-hexane, azole 3 in toluene, azoles 2 and 4 in methanol, azoles 5 and 6 in water, and azole 7 in THF, whereas the highest values of λ max were obtained for azoles 1, 7, and 8 in DMF, azoles 4-6 in toluene, azole 2 in chlorobenzene, and azole 3 in CCl 4 . The use of the MPW1PW91 functional led to the following lowest λ max relative to solvents: azoles 1-4, 7, 8 (n-hexane), azole 5 (methanol), and azole 6 (acetone), while the highest λ max values were detected for azoles 1 and 8 in DMF, azoles 2, 4, and 7 in DMSO, azoles 3 and 6 in toluene, and azole 5 in chlorobenzene. The DFT  computations using the WB97XD functional resulted in the lowest λ max values for azoles 1-4 in methanol, azoles 5 and 7 in water, and azoles 6-8 in n-hexane, while the highest λ max values were found for azole 1 in acetonitrile, azole 2 in n-hexane, azoles 3-5 in toluene, azole 6 in chlorobenzene, and azoles 7-8 in DMF. The above data show that the use of various functionals for the calculations of λ max for azoles 1-8 in a large spectrum of solvents that differ in polarity results in dissimilar values of λ max depending on the reaction field of solvent. However, we noticed that for each of the functionals, the lowest λ max in the quinoline 8 UV spectrum was detected for n-hexane, whereas the highest one was observed in DMF. Moreover, for some of the studied hetarenes, the lowest values of λ max were observed in methanol, while the highest λ max values were detected for toluene or chlorobenzene. Next, we compared the λ max values in the theoretical UV spectra of azoles 1-8 that were obtained by the PCM approach with those in the gaseous phase. The results are listed in Tables 10-13, where the negative Δλ max values correspond to a hypsochromic shift, whereas the positive values correspond to a bathochromic shift in relation to λ max generated in vacuo.
Using the B3LYP functional, the largest hypsochromic shift was observed for chlorine derivative 1 in CCl 4 , toluene, and n-hexane, whereas the largest bathochromic shift was detected for pyrazolopyrazole 7 in all solvents excluding CCl 4 , toluene, and n-hexane (Table 10).
By applying the CAM-B3LYP functional, the largest hypsochromic shift was noticed for carbazole 6 in all solvents, while the bathochromic shift was detected for compound 7 in acetonitrile, acetone, methanol, DMF, DMSO, and water.
The use of the MPW1PW91 functional resulted in a significant hypsochromic shift for pyrrole 2 in all solvents, whereas a weak bathochromic shift was observed for azoles 1, 3, 4, and 8 also in all solvents.
The computations with the use of the WB97XD functional resulted in the largest hypsochromic shift for chlorine derivative 1 in all solvents apart from acetonitrile and THF, while only small bathochromic shifts were observed for the remaining azoles.
In the last step of work, we focused our attention on the GCRD factors as important parameters that give insights into chemical reactivity and stability of the molecules based on the DFT formalism and Koopmans theorem [54,55]. From these parameters, electronegativity (χ), chemical hardness (η), first electronic potential, and electronic potential were especially interesting for us [56,57]. The aforementioned parameters are defined by the following equations: μ = −(I + A)/2 and η = (I − A)/2, and electronegativity χ = (I + A)/2, where I is the first ionisation potential (I = −E HOMO ) and A is the electron affinity (A = −E LUMO ). The chemical hardness is especially extensively used for estimation of the chemical behaviour of the molecules and could be useful for characterisation of organic compounds [58]. The obtained electronegativity values, and consequently the chemical hardness, were comparable with the literature data for similar azole moieties [23]. Basically, in the case of "hard" molecules the HOMO−LUMO gap remains large. For the so-called "soft" derivatives, the character of such gap is counterdirected. From the quantum mechanics standpoint, molecules with a small HOMO−LUMO gap and low excitation energies will have their electron density changed more easily than hard analogues.
The details concerning the analysis of GCRD descriptors for compounds 1-8 (in the gaseous phase, as well in the PCM solvation model for all solvents used) are listed in Tables S17-S73 (Supplementary material). When we arranged the results involving all functionals and solvents in increasing order, we observed small discrepancies between them. The analysis of the matrix covering all values led to a conclusion that carbazole 6 had the lowest chemical hardness η for all functionals and almost all solvents used. On the other hand, the fused pyrazole derivative 7 had the highest value of hardness. Considering hardness, the investigated compounds can be ordered as follows: 6 < 5 < 2 < 4 < 3 < 1 < 8 < 7. For the MPW1PW91 functional and n-hexane as a solvent, compounds 6 and 5 interchanged in the above list. Considering the HOMO-LUMO gap and hardness η, the aforementioned computations showed that carbazole 6 and indole 5 are relatively soft chemical species, thus more reactive than the other studied compounds, whereas pyrazolopyrazole 7 displayed the highest hardness.

Conclusions
In the above discussion of results, we show that the pyrrolic nitrogen atoms in all studied compounds as well as the pyridinic nitrogen atom in the quinoline ring within the structure of 8 are potential centres of interactions in the protein binding pockets. The NBO analysis proves that the pyrrolic nitrogen atom has the main share in the σ N-H bond. This nitrogen atom is sp 2.31-2.33 hybridised, thus it has more than 69% p character.
The highest stabilisation energy E(2) in the relationship NBO donor (σ N-H )-NBO acceptor corresponds to the bond between the pyridinic nitrogen atom and neighbouring carbon in the indazole ring.
The estimation of free energy of solvation (SMD model, 13 different solvents, functionals: B3LYP, CAM-B3LYP, MPW1PW91, and WBXD91) shows that for azoles 1, 2, 6, and 7, this parameter has the lowest negative value in the reaction field of acetonitrile, whereas for compounds 3-5 and 8 these values are the lowest using acetone as a solvent. Considering all solvents used for the investigations, the lowest free energy of solvation in the SMD model can be attributed to carbazole 6, while the highest one to fused compound 7. The observed differences substantiate the dissimilar abilities of carbazole 6 and fused pyrazole 7 for interactions within the protein pocket that were reported in our previous study [22]. Examining the solvatochromic effect for azoles 1-8 with the use of the PCM solvation model and B3LYP, CAM-B3LYP, MPW1PW91, and WBXD91 functionals in 13 solvents, we observed that the UV spectrum of quinoline 8 had the lowest λ max in n-hexane, whereas the highest value was noticed for DMF as a solvent.
The investigations of the HOMO-LUMO gap and η have led to a conclusion that carbazole 6 and indole 5 can be regarded as soft molecules, while fused diazole 7 is a hard molecule.

Supporting information
Cartesian coordinates of all discussed ligands and adducts, MEP analysis visualisation, charge computation, free energy of solvation, and TD-DFT calculation results are given in the Supplementary file. This information is available via the Internet or upon the request from the corresponding author.