On the Interactions of Fused Pyrazole Derivative with Selected Amino Acids : DFT Calculations

Due to the increasing prevalence of neoplasms, there is a permanent need for new selective cytostatic compounds. Anticancer drugs can act in different ways, affecting protein expression and synthesis, including disruption of signaling pathways within cells. Continuing our previous research aiming at elucidating the mechanism of pyrazole’s anticancer activity, we carried out in silico studies on the interactions of fused pyrazole derivative with alanine, lysine, glutamic acid, and methionine. The objective of the study is to improve our understanding of the possible interactions of pyrazole derivatives with the above-mentioned amino acids. For this purpose, we apply the DFT formalism (optimization using the B3LYP, CAM-B3LYP, PBE0, and M06L functionals) and interaction energy calculations (counterpoise corrected method based on the basis set superposition error, BSSE) together with QTAIM approach and estimation of the HNMR chemical shifts of analyzed pyrazole derivative using different basis sets and DFT functionals in CPCM solvation model (and water used as a solvent).


Introduction
Currently, there is a growing demand for new cytostatic compounds that act selectively.Drugs with anticancer properties exert them through different modes of action [1].These include disruption of DNA replication and thus inhibition of cell proliferation and inhibition of translation of target proteins, such as kinases that regulate signaling pathways within cells.Eventually, these compounds induce apoptosis of neoplastic cells.So far, numerous oncogenes have been described and many of them are dysregulated protein kinases.Decreasing the activity of upregulated kinases results in a stunted tumor growth [2][3][4][5].This is why kinase inhibitors have attracted significant attention for targeted therapy of many cancers.
Heterocycles containing pyrazole are a part of a group of anticancer drugs on which we focused our attention.This type of anticancer compounds includes those which interact with active sites of kinases, for example, through binding with amino acids.Alanine (Ala), lysine (Lys), glutamic acid (Glu), and methionine (Met) existing in these cavities are especially important as they can interact with pyrazole derivatives [2][3][4][5].
According to the results obtained by our group previously, certain condensed pyrazole derivatives exert a cytostatic activity and induce apoptosis in the HT29 human colorectal neoplastic cells [6,7].In our subsequent papers, we hypothesized that they could disrupt the intracellular Mg 2+ ions balance through the action of pyridinic nitrogen within the pyrazole moiety [8,9].An NMR experiment elucidating the interaction of chosen amino acids and the condensed pyrazole derivative showed that this interaction primarily happens between the pyrrolic nitrogen atom (NH) and carbon fragment B (Scheme 1) [7][8][9].
Continuing our studies that aimed at elucidating the potential mode of pyrazole cytostatic activity, we conducted further theoretical studies investigating simulated interactions of fused derivative 7 (synthesis is given in Scheme  A-E: protons of the final product considered for NMR investigations) with certain amino acids (Ala, Lys, Glu, and Met).The objective of the study is to improve our understanding of the tendency regarding mechanism of interactions between pyrazole derivatives and selected amino acids.We focused our attention on the electrostatic contacts, especially hydrogen bond formation between the azole and amino acids.To achieve that, we decided to choose the DFT formalism and interaction energy calculations with QTAIM approach and calculation of the 1 H NMR chemical shifts of compound 7 with several basis sets and DFT functionals in conductor-like polarizable continuum model (CPCM) of solvation, so that the influence of the solvent environment in cell is also accounted for.What is more, we aimed to determine the methods and basis sets that gave the closest results to the experimental shifts in the 1 H NMR spectrum and the most stable and preferable conformation of compound 7.

Materials and Methods
We carried out density functional calculations and optimized the compound geometries at the DFT [10] level of theory using the Gaussian 09 D.01 program [11], B3LYP/6-31G(d,p) [12], (b) CAM-B3LYP/6-31G(d,p) [13], (c) B3LYP/6-311+G(d,p) [14,15], (d) PBE0/6-31G(d,p) [16,17], and (e) M06L/6-31G(d,p) [18] approaches in the CPCM of solvation (water molecules as solvent).The ideal gas, rigid rotor, and harmonic oscillator approximations were applied to calculate the vibrational frequencies and thermodynamic properties.We obtained no imaginary frequencies within the vibrational spectra and therefore confirmed the energy minimum for all conformers.We created the conformers through the rotations of C1-S1-C2 bond of azole and bonds connecting N-H moiety of azole with carboxyl or amine group of the particular amino acid in dihedral angle increments of 20 ∘ .AIMAll Professional package for topological analysis was used for QTAIM calculations [19].For this purpose, the "Proaim" basin integration approach was used using "superfine" interatomic surface mesh and "very high" outer angular quadrature.Thanks to the atomic (integrated) Laplacian values, we could monitor the accuracy of basin integrations.What is more, the Poincaré-Hopf rule was satisfied for all carried out calculations.The  and ∇ 2  parameters at the bond critical points (BCPs) concerning dipole polarizability tensors were calculated based on the optimized structures (taking into consideration the already calculated interaction energy) using the AIM approach and applying B3LYP/6-31G(d,p) (adducts 7AlaB, 7LysB, 7GluB, and 7MetB), B3LYP/6-311+G(d,p) (adducts 7GluB6 and 7MetB6), PBE0/6-31G(d,p) (adducts 7LysP), and M06L/6-31G(d,p) (adduct 7AlaM) levels of theory, respectively.The components of molecular dipole polarizability tensors were calculated as the set of first derivatives of the dipole moment components with respect to the components of an applied electric field.We calculated the value of NMR shielding for proton (H ref ) in TMS at the B3LYP/6-31G(d,p) level of theory (CPCM of solvation and water as solvent).The abovementioned method was also used for calculation of TMS (reference compound) and the investigated 7-amino acid adducts.Then the values obtained for TMS allowed calculating the chemical shifts of 7 and 7-amino acid adducts, in accordance with the following equation:   =  ref −   , with   being chemical shift of -nuclei of 7-amino acid adducts and  ref and   being the calculated isotropic magnetic shielding tensors of TMS and 7-amino acid adducts, respectively [8,[20][21][22].An average of the estimated chemical shifts was calculated for the methyl groups protons A or E and also these of C or D (Scheme 1).We calculated interaction energy of all optimized adducts applying counterpoise corrected method based on the basis set superposition error (BSSE) at B3LYP/6-31G(d,p) and B3LYP/6-311++G(2d,3p) levels of theory, respectively [23,24].
The ChemCraft 1.7 software was utilized for the visualization of all optimized conformers and adducts [25].The calculations were carried out using resources provided by the Polish Grid Infrastructure (PL-Grid) and Wrocław Center for Networking and Supercomputing (Bem cluster, WCSS Grant number 327/2014).

Results and Discussion
The analysis of azole 7 coordination with magnesium ions proved that the pyrrolic nitrogen NH with a labile exchangeable proton is the principal position in this interaction [8].Moreover, the pH alterations are unimportant for small molecules that, like compound 7, contain practically one protonation site.Additionally, we prove that such compounds do not tautomerize (tautomerism involving protonation of the nitrogen atoms can be neglected) [7][8][9].Based on this information, we decided to examine the interaction of the pyrrolic nitrogen with the following amino acids: alanine, lysine, glutamic acid, and methionine.The geometry of the azole-amino acids adducts was optimized using four different functionals, including hybrid, GGA, and meta-GGA functionals, that is, B3LYP, CAM-B3LYP, PBE0, and M06L.The data concerning the lowest (in bold) and highest interaction energies (CPCM of solvation and water as solvent, counterpoise corrected method including the basis set superposition BSSE error, Gaussian G09 D.01 software [11]) for the azole-amino acids adducts are shown in Table 1.
Interaction energy between two atoms or molecules A and B can be generally calculated to be a result of subtraction of the energies of components A and B from the energy of the fused AB structure: The   index reflects the geometry of the fused AB structure, whereas the   index reflects that of separate A and B. However, this equation often results in incorrect estimates, particularly when hydrogen bonding or dispersion interactions are present within the fused structures.The error associated with dispersion interactions can be lessened by using smaller basis sets, which make the AB structure more stable than A and B separately using the basis set superposition error (BSSE).This result is achieved thanks to the expansion of the wavefunction over fewer basis functions for A and B molecules compared to the AB structure [23].Using very big basis sets with the BSSE is impossible for the majority of investigated compounds.Therefore, counterpoise method (CP) can be applied, which can estimate the BSSE size.Then the AB structure is described in the same way, but A and B have basis sets which are of equal size, the same as that of the basis set for the AB structure.Thus, interaction energy is calculated using CP in the following equation: The AB indices in (2) show that absolute basis sets for the calculations concerning all the components are the same and are equal to those of the AB structure [23].Hence, it can be concluded that such methodology significantly better describes the interaction phenomenon within the analyzed herein adducts.
Although the interactions between amino acid carboxylic group and heterocyclic nitrogen are significantly stronger  than the interactions involving amino group of amino acid and heterocyclic nitrogen [20], the adduct 7-Ala was additionally optimized regarding the interaction of the alanine amino group and pyrazole nitrogen atoms (adduct 7AlaB-N, B3LYP/6-31G(d,p) level of theory; Figure S3, Supplementary Material).
The interaction between the amine group of lysine and azole 7 pyrrolic nitrogen (adduct 7LysP-N, Figure S6 The above discussion led to a conclusion that adduct 7LysP (optimized at the PBE0/6-31G(d,p) level of theory) involving azole 7 and lysine has the strongest linking due to the interaction between the amino acid carboxylic group and pyrrolic nitrogen of compound 7.
The optimization of the adduct 7MetB-N at the B3LYP/6-31G(d,p) level of theory (Figure S11, Supplementary Material) resulted in a relatively low interaction energy, that is, −8.42 [interaction energy calculated at the B3LYP/6-31G(d,p) level of theory] or −7.79 kcal mol −1 [interaction energy calculated at the B3LYP/6-311++G(2d,3p) level of theory].The interaction involved the methionine  amine group and the condensed pyrazole system led to two relatively

AIM Analysis.
The application of AIM method includes determination of type of bond (ionic or covalent) and its properties [26] also when it comes to hydrogen bonding.In our AIM calculations, we focused on the interactions of the amino acids carboxylic group with condensed azole 7. Particular attention was devoted to adducts with the highest and lowest interaction energies.
As a principle, there is a relation between the  parameter at the bond critical point (BCP) and the strength of the hydrogen bond.The fact whether the Laplacian ∇ 2  at the BCP is positive or negative determines whether the hydrogen bond dominates the covalent (∇ 2  < 0) or electrostatic (∇ 2  > 0) interactions.The values for the electron density topological parameters  and ∇ 2  (related with a physical basis via local statement of the virial theorem for bonding interactions [26]) for the analysed adducts are given in Tables 2 and S1-S3 (Supplementary Material).It is noteworthy that the AIM parameters in the literature were correlated with the strength of analyzed bond [27]; however, in cited paper only one bond within the structure of the analyte was considered, whereas in our models we have detected even two hydrogen contacts within the structure of the particular azole-amino acid rotamer being analyzed.
In case of all analyzed adducts, the absolute value of bond ellipticity () was ca 0.01-0.04,which indicates cylindrical, directional character of investigated contacts.
The ∇ 2  values were positive but relatively small; thus they indicate that the adducts were connected through strong hydrogen bonds of electrostatic character.Moreover, the hydrogen bonds of adducts azole 7-amino acid with lower (more negative) interaction energies have higher values of  1).This demonstrates that adducts 7AlaB, 7LysB6, LysM, and 7MetB are connected by stronger hydrogen bonds in comparison with adducts 7AlaM, 7LysP, 7GluB6, and 7MetB6, respectively.

NMR Analysis.
The mechanism of interactions of metal ions and organic compounds, such as synthesized drugs and proteins, is of utmost significance from the biological point of view.Various analytical methods can be applied to investigate them, including 1 H NMR spectroscopy.Noteworthy is the fact that NMR spectroscopy has found its use for the determination of interactions of some ligands with proteins and allows assessing quantitatively these interactions if the equilibrium constants and changes in chemical shifts are taken into account [8,20,21,26,[28][29][30][31].As a result, we decided to focus only on 1 H NMR spectroscopy in our study.Albeit compound 7 has four nitrogen atoms within its structure, the low abundance of the 15 N isotope in comparison with the 1 H isotope as well as the significant signals broadening due to the large quadrupole moment of 14 N renders the nitrogen NMR spectroscopy impractical [8,30,31].Considering the low utility of 13 C NMR spectroscopy for the investigations of interactions between metal ions and small molecule organic compounds [8,9,29] as well as the low natural abundance of 13 C, we decided to use 1 H NMR for the studies depicted in our paper.The theoretical 1 H NMR spectrum was generated for all azole-amino acid adducts using the gauge-including atomic orbital (GIAO) method [22] as implemented in Gaussian G09 D.01.In order to maintain clarity, we considered four adducts with the lowest (more negative) energy, previously optimized in CPCM of solvation (water as solvent) at the B3LYP/6-31G(d,p) (Tables S4 The theoretical values of chemical shifts correlate well with the experimental ones for pyrazole derivative 7, which is shown by low mean absolute deviation (MAD) [8].These low MAD values for the theoretical 1 H NMR spectra of the adducts azole 7-amino acids (Tables 3 as well as S4, S16, and S18; Supplementary Material) correspond to the method and basis used for optimization of the adducts 7AlaB (B3LYP/6-31G(d,p) approach), 7LysP (PBE0/6-31G(d,p) approach),  1) and emphasize the importance of hydrogen bonds within these adducts.
From theoretical point of view, fascinating conclusions can be drawn from the calculations involving spin-spin coupling constants between the nuclei of hydrogen bond donor and acceptor as calculated at the B3LYP/6-311++G(2d,3p) level of theory.The conclusions are related with the adducts with the highest (more negative) and lowest (more positive) interaction energies (Figures S12-S19 and Table S23, Supplementary Material).For the interactions of azole 7 and alanine, the coupling constants involving atoms connected by hydrogen bonds O-H⋅ ⋅ ⋅ N have values of ca −2 Hz.The increase in length of this bond in adducts 7AlaB and 7AlaM leads to an increase in both the torsion angle  of the hydrogen bond donor and negative value of the dihedral angle .Conversely, the increase in length of C-H⋅ ⋅ ⋅ O=C bond (adducts 7AlaM and 7AlaB) resulted in an insignificant decrease of the angle  and more negative value of the dihedral angle  as well as a decrease in coupling constant from 2 to 1 Hz.It should be added that the hydrogen contacts in adduct 7AlaB are characterized by greater values of electron density  calculated by the QTAIM method (Table S1, Supplementary Material).
The N-H⋅ ⋅ ⋅ O=C hydrogen contact in azole 7-lysine adduct (adducts 7LysB and 7LysP) can be characterized by comparable values of the bond length, torsion angle, and coupling constant for hydrogen bond acceptor-donor (ca 5 Hz), as well as electron density  (Table 2).However, the dihedral angle  has a negative value for adduct 7LysB.
The higher interaction energy of adduct 7LysP in comparison with complex 7LysB is probably due to the already noted smaller distance of the hydroxyl group (carboxylic function) of lysine from the pyrazole N2 nitrogen atom, as well as the second hydrogen contact O-H⋅ ⋅ ⋅ N.
Comparison of the interaction energy of 7Met adducts 7MetB and 7MetB6 reveals that the higher energy of adduct 7MetB is connected with the presence of a second hydrogen bond O-H⋅ ⋅ ⋅ O 2 S, despite the fact that this contact has lower electron density  and is longer compared to adduct 7MetB6 (Tables S3 and S23, Supplementary Material).However, this is associated with a smaller value of torsion angle  and dihedral angle  as well as coupling constant that decreases to ca 4 Hz (adduct 7MetB) in comparison with 7 Hz (adduct 7MetB6).
The interaction of azole 7 with glutamic acid leads to adducts 7GluB and 7GluB6 that are linked through a N-H⋅ ⋅ ⋅ O=C hydrogen bond.This hydrogen contact has similar values of the bond length, torsion angle , and coupling constant (ca 4 Hz, Table S23, Supplementary Material), although the dihedral angle  (O-H⋅ ⋅ ⋅ N) in adduct 7GluB is shorter and has somewhat higher value of torsion angle  and lower coupling constant, that is, −1.18 Hz versus −9.38 Hz, in comparison with the same bond parameters in adduct 7GluB6.It should be added that the electron densities  in the paths of both hydrogen bonds between the nuclear attractors of adduct 7GluB attain higher values.

Conclusions
The paper has focused attention on the interactions of the cytotoxic condensed derivative of pyrazole 7 with selected amino acids present in the receptor cavities of human protein kinases.The most significant interactions involve the amino acids carboxylic groups with azole 7 pyrrolic nitrogen.Results of interaction energy calculations among the 7-amino acid adducts regarding method and basis set used for initial optimization protocol for each of the analyzed adducts are gathered in Table 1.From this point of view, we observed that for 7Ala, 7Glu, and 7Met adducts the optimization at the B3LYP/6-31G(d,p) level of theory was undoubtedly a reasonable choice.However, for the 7Lys adduct, the PB0 functional resulted with optimized structure, within which more significant interactions were detected.
We have shown that the low MAD values for the theoretical 1 H NMR spectra of the adducts azole 7-amino acids correspond to the method and basis used for optimization of the adducts 7AlaB, 7LysP, 7GluB, and 7MetB6 with the highest interaction energies.
The results presented in the previous sections show that, for adducts with similar geometric parameters and coupling constants, the electron densities  of the paths of bonds between nuclear attractors (adducts 7Lys) have a decisive influence on the interaction strength.For the remaining adducts, that is, 7Ala, 7Glu, and 7Met, a second hydrogen bond significantly affects the interaction.Furthermore, an increase in the hydrogen bond length in these adducts is accompanied by a decrease in the value of donor torsion angle  and coupling constant.

Table 2 :
Electron densities  and Laplacian fields ∇ 2  [e Å] determined for the 7Lys adducts with a hydrogen bond; nd: no data.