Abstract

Quantitative structure-property relationship (QSPR) study on the acid dissociation constant, pKa of various 22 N-base ligands including pyridines, pyrimidines, purines, and quinolines has been carried out using Codessa Pro methodology and software. In addition, the formation constant, Kc of these ligands with Pt(II) (bpy = 2,2′-bipyridine) ion has also been modelled with the same methodology. Linear regression QSPR models of pKa and Kc were established with descriptors derived from AM1 calculations. Among the obtained QSPR models of pKa presented in the study, statistically the most significant one is a four parameters linear equation with the squared correlation coefficient, values of ca. 0.95 and the squared cross-validated correlation coefficient, values of ca. 0.89, and external the squared correlation coefficient, values of ca. 0.97. Statistically the most significant QSPR model of Kc is also a four parameters linear equation with the squared correlation coefficient, values of ca. 0.75 and the squared cross-validated correlation coefficient, values of ca. 0.55, and external the squared correlation coefficient, values of ca. 0.81. An analysis of descriptors that involved in the pKa models indicate that reactivity index and charge distribution related descriptors play major roles to model acid dissociation constant of ligands of N bases.

1. Introduction

The acid dissociation constant, pKa, which describes the extent to which a compound dissociates in an aqueous solution, is a fundamental physical property of a chemical. The pKa value plays a fundamental role in many analytical procedures such as acid-base titrations, solvent extraction, complex formation, and ion transport, and is especially relevant in medicinal chemistry because it affects ADME and activity. The degree of ionization determines permeability and solubility, two properties widely used in pharmaceutical research to predict the pharmacokinetic profile of a compound. Differences in adsorption, toxicology, solubility, bio-concentration, and reactivity are common when comparing the properties of the ionized molecule to its neutral form [1]. All of these reasons make the pKa a common physico-chemical property as a descriptor for QSAR and QSPR (quantitative structure and property relationship) studies. The experimental pKa values of organic compounds have been largely determined through well-established methods such as potentiometry, UV-visible absorption spectrometry, conductimetry, and competitive reactions. A pKa value of a compound may have to be determined, as it may not be available in published literature. Therefore, it is of interest to develop modeling methods for estimating the pKa of compounds, and use these methods to predict the properties of a chemical in a solvent environment. Appropriate modeling methods may minimize cost and time to predict accurately pKa values of unknown compounds.

Numerous computational studies have attempted to predict pKa of organic compounds by applying different theoretical models [213]. Three main approaches have been employed in these theoretical models. First approach relieds upon using linear free energy relationships (LFER). The primary disadvantage of the purely LFER-based approaches is the need to derive a vast number of fragment constants and correction factors which are used in the estimation methods. Some researchers have also criticized the use of LFER approaches on more philosophical grounds, arguing that the prediction of molecular properties by fragment constant methods lacks solid scientific support [14, 15]. A second approach is based on accurate energy calculation of the molecules in their neutral and deprotonated forms in the gas phase and solvent, to account for the solute-solvent interaction. In order to obtain adequately accurate energies to calculate solution phase dissociation constants, one must account for electron correlation at the ab initio or density functional theory (DFT) methods, and consider the effects of solvation on the molecule. When moderate- or relatively large-sized molecules are considered, due to the great computational cost of ab initio or DFT calculations this approach is not feasible, especially for virtual screening applications. Third is the QSPR approach that is a mathematical equation relating chemical structure to a wide variety of physical, chemical, biological, and technological property.

In the study herein, one of the pKa modeling approaches described above, the QSPR method has been used to construct quantitative models to predict pKa values of various 22 N-bases ligands including pyridines, pyrimidines, purines, and quinolines. In addition, we have developed QSPR models to predict the formation constant, Kc of reactions of these ligands with Pt(II) (bpy = 2,2′-bipyridine) ion. After the discovery of the cell division-inhibiting effect of cisplatin by Rosenberg et al. [16], there is considerable interest in platinum chemistry. To date, six platinum complexes (including cisplatin) have been used in drug anticancer chemotherapy [17]. To obtain pKa and Kc QSPR models, program Codessa Pro (comprehensive descriptors for structural and Statistical Analysis), Version 2.7.2 [18] was employed to build a multilinear regression (MLR) method using observed pKa and Kc values taken from Kawanishi et al. [19] in combination with the descriptors which were calculated by the AMPAC [20] semiempirical quantum chemistry code.

2. Results and Discussion

The structures of 22 N-base ligands are shown in Figure 1. Table 1 shows the following information QSPR modeling of pKa of compounds: (i) AM1-based calculated molecular descriptor values involved in the obtained models; (ii) experimental pKa values taken from the original references; (iii) the predicted pK values using Model I, II, and III obtained in this study. Table 2 shows the obtained QSPR equations together with their statistical parameters for pKa modeling. Table 3 shows the intercorrelation of descriptors involved in the pKa models. The predicted values of pKa using Model I is plotted versus the experimental values in Figure 2. Reliability of Model I was tested by the Y-randomization test. 250 random shuffles of the Y (pKa) were chosen, and the modeling process was performed for all cases of Model I. Results are shown in Figure 3. It should be noted that all of the results of the 250 random shuffles of the Y were not included, only the highest 60 data points were taken in the Figure 3. The lower values of and in comparison with the real model’s results support the hypothesis that the good statistical results obtained by the QSPR model are not due to a chance correlation, or structural dependency of the training set.

Table 4 shows the following information for Kc QSPR modeling: (i) AM1-based calculated molecular descriptor values involved in the models; (ii) experimental Kc values taken from the original references; (iii) the predicted Kc values using Model IV, V, and VI obtained in this study. Figure 4 shows plot of experimental versus calculated Kc using Model IV. Table 5 shows the obtained QSPR equations together with their statistical parameters for Kc modeling. Table 6 shows the intercorrelation of descriptors involved in Kc models. At this point, it should be observed that Tables 3 and 6 demonstrate that although some descriptors exhibit relatively high intercorrelation, they are not involved in the same models as seen in Tables 2 and 5. Reliability of the best Kc model (Model IV) was also tested by the Y-randomization test. 250 random shuffles of the Y (Kc) were chosen, and the modeling process was performed for all cases of Model IV from which the results are shown in Figure 5. It should be noted that all of the results of the 250 random shuffles of the Y were not included, only the highest 60 data points were taken in the Figure 5. The lower values of and in comparison with the real model’s results support the hypothesis that the good statistical results obtained by the QSPR model are not due to a chance correlation, or structural dependency of the training set.

2.1. Analysis of pKa Models

A perusal of Table 2 shows that nine types of descriptors are involved in the three pKa models. Among the three models, the best, statistically speaking, is Model I, as evident from the fact that it has a very good statistical fit (, , ) and predictive ability (, , RMSE = 0.51) parameters. By analyzing Model I, it is clear that the most representative descriptor is HNMVF with the positive coefficient and the highest t-test value. The magnitude of HNMVF of a molecule comes from O-H, N-H, or C-H stretching vibrations. Hybridization state (sp, sp2, and sp3) also affects the magnitude of HNMVF of a molecule. An increase in the magnitude of HNMVF favors the exhibitions of the pKa of a N-base ligands. In this model, FNSA3 and FPSA which are charged partial surface area (CPSA) descriptors, also have positive coefficients. They describe the polar interactions between the molecules. MNRI-C has a negative coefficient, highlighting that the lower the magnitude of MNRI-C, the lower the pKa. This descriptor estimates the relative reactivity of the atom (carbon) in the molecule, and is related to the activation energy of the corresponding chemical reaction.

Model II is also a four-parametric model, similar to Model I. This model also has very good statistical fit (, , and ) and good predictive ability (, , RMSE = 1.11) parameters. The descriptors involved in this model are MERI-C, HDSA1/TMSA, WPSA2, and . In this model, only HDSA1/TMSA has a positive coefficient. HDSA1/TMSA is statistically the most significant descriptor, as evident from the fact that it has the highest t-test value. HDSA1 is the hydrogen bonding donor ability of a molecule and TMSA is total molecular surface area of a molecule. An increase in the magnitude of HDSA1/TMSA favours the exhibitions of the pKa of N-base ligands. MERI-C is the the second statistically significant descriptor with a negative coefficient. This is a relative reactivity index, and corresponds to activation energy of a corresponding chemical reaction. WPSA2, and are related to charge distribution of the molecules. In this model, negative coefficients of MERI-C, WPSA2 and indicate that decreases in the magnitudes of MERI-C, WPSA2 and are not favorable for an increase of pKa of the molecules.

Model III is a two-parametric model. HDCA1/TMSA and FNSA descriptors are involved in this model with positive coefficients. These are CPSA descriptors, and describe the polar interactions between the molecules. Although, the goodness of fit of this model is only satisfactory, evident from the fact that it has , , and , its external predictive power is very good (). It is worthy here to mention that model validation is historically a hot topic for discussions within the QSAR community. Some researchers use only internal for model validation, and some others strongly believe that external validation is vital for a reliable QSAR model. If one surveys the QSAR literature, it can be shown that many QSAR models have been published only with internal validation until a study published by A. Golbraikh and A. Tropsha, entitled “Beware of q2!” [21]. In that study, the authors reexamined several published QSAR data sets, and demonstrated that there is lack of any relationship between and . They concluded that the high value of appears to be the necessary, but not sufficient condition, for the model to have a high predictive power and the external validation is the only way to establish a reliable QSAR model. Our results regarding the pKa models support these conclusive remarks. Model II has a higher internal validation value, , but its external validation value, is only 0.85, whereas Model III has a lower internal validation value, , and its external validation value, is very good 0.91.

2.2. Analysis of Kc Models

A perusal of Table 5 shows that nine types of descriptors are involved in the three Kc models. Among the three Kc models, Model IV has good fit (, , ) parameters, whereas its predictive ability only is acceptable (, , RMSE = 3.18). BETA Polarizability, Internal Entropy, ESP-FNS3, and Min e-e REP-N are involved in this model.  BETA Polarizability which has a positive coefficientisthe most representative descriptor,as evident from the fact that it has the highest t-test value. Internal Entropy divided by number of atoms in the molecules which has a negative coefficient in the model, is the least representative descriptor, as evident from the fact that it has the lowest t-test value. Also, the positive sign of coefficient of ESP-FNS3 and Min e-e REP-N indicate that increases in the magnitude of ESP-FNS3, and Min e-e REP-N are favorable for an increase of Kc of Bis(2,2′-bipyridine)platinum(II)-N-base adduct reactions. Mechanistic interpretation of the model is quite complex due to the diverse nature of the involved descriptors. BETA Polarizability reflects information about possible inductive interactions in the molecule. It also characterizes the properties of a molecule as an electron acceptor. ESP-FNS3 is a CPSA descriptor, and describes the polar interactions between the molecules. Min e-e REP-N describes the electron repulsion driven processes in the molecule, and may be related to the atomic reactivity in the molecule [22].

In Model V, BETA Polarizability and Internal Entropy have the same sing of coefficients as Model IV. Internal Entropy is the most representative descriptor in this model. The third descriptor is the Min Re. E-CN which is an energy-related descriptor with a positive sign of coefficient. Statistical parameters of this model are similar to Model III in term of validation parameters. Although this model has the lowest internal validation () value, its external predictive ability, is relatively high. When these six Kc and pKa models are considered together, one could draw a conclusion that there is not any relationship between and .

The final model, Model VI has four descriptors. Internal Entropy is also the most representative descriptor in this model. Internal Entropy and avg. 1-e RI-C have a positive sign of coefficient whereas max n-n Re.-CN and DPSA1D in CPSA have negative sign of coefficient. As observed for Model IV and V, an increase of in the magnitude of Internal Entropy divided by number of atoms in the molecules favours the exhibitions of the Kc for Bis(2,2′-bipyridine)platinum(II)-N-base adduct reactions. It may be concluded that the entropy changes seem to have an important effect on the Bis(2,2′-bipyridine)platinum(II)-N-base adduct reactions. Although, it has not been included as a Model in Table 5, but Internal Entropy divided by number of atoms in the molecules, itself as a monoparametric model has an acceptable statistical parameter, , , ,  , and RMSE . It is worthy to mention that all of the Kc models have one nitrogen-related descriptor (Min e-e REP-N, Min Re. E-CN and Max n-n Re.-CN). When the proposed Bis(2,2′-bipyridine)platinum(II)-N-base adduct by formation of N-Pt(II) bond [19] is considered, the appearance of N-related descriptors in all of the models is not a surprise.

Finally, as mentioned in the introduction section, pKa can be seen as a descriptor for QSAR and QSPR studies in literature. In this study, pKa has been tested as a descriptor for the modeling of Kc, but it has failed during the preselection of the descriptors by the software (Codessa Pro) due to its very low squared correlation coefficient of the one-parameter equation is less than 0.01 by default.

3. Experimental Section

3.1. Quantum Chemical Calculations

All structures were optimized without geometry constraints using the standard AM1 Hamiltonian [23] within the AMPAC quantum chemistry code [20]. Structure fundamental vibrations were also calculated using the same method to check if there were true minima. All computations were carried out for the ground states of these molecules as single states. Output files of the molecules from AMPAC code were used as the input file for Codessa Pro [18] code for descriptor generation.

3.2. Descriptors Generation and Their Definitions

In the present work, more than five hundred descriptors were exploited by using Codessa Pro code, and they were divided into groups such as constitutional, topological, geometrical, electrostatic, quantum chemical, thermodynamic, and contracted. Constitutional descriptors are related to the number of atoms and bonds in each molecule. Topological descriptors include valence and nonvalence molecular connectivity indices calculated from the hydrogen-suppressed formula of the molecule, encoding information about the size, composition, and the degree of branching of a molecule. Geometrical descriptors are calculated from 3D atomic coordinates of the molecule and comprise moments of inertia, shadow indices, molecular volumes, molecular surface areas, and gravitation indices. Electrostatic descriptors reflect characteristics of the charge distribution of the molecule. Quantum chemical descriptors encode the polar interactions between molecules or their chemical reactivity and the activation energy of the corresponding chemical reaction. Thermodynamic descriptors are quantum mechanically calculated on the basis of the total partition function of the molecule, and its electronic, translational, rotational, and vibrational components. Codessa Pro also allows one to construct new descriptors by using the existing descriptors. In this way, the author constructed some common quantum chemical indices, namely, chemical hardness, electronegativity, and electrophilicity from HOMO and LUMO orbital energies. More than five hundred descriptors, andseventeen types of descriptors were involved in the selected pKa and Kc models as shown in Table 1 and Table 4. Their definitions are adopted from the CODESSA Pro Manual [24] as described below.

Highest normal mode vibrational frequency (HNMVF) is an extreme maximum value of the normal mode vibrational frequencies in the molecule. Definition of the normal mode of vibration arises from the quantum mechanical Harmonic Oscillator model of a diatomic molecule.

Fractional PNSA (PNSA3/TMSA) or FNSA3 is one of the electrostatics charged partial surface area (CPSA) descriptors, and was invented by Jurs et al. [25, 26]. Definition of FNSA3 is atomic charge weighted partial negative surface area (PNSA) divided by total molecular surface area (TMSA) which has related features responsible for polar interactions between molecules.

Fractional CPSA (PPSA1/TMSA), FPSA1 is another CPSA descriptor and its definition is given as partial positive surface area (PPSA) divided by TMSA.

Maximum nucleophilic reaction index for a C atom, MNRI-C is one of the quantum-chemically calculated charge distribution-related reactivity indices. These descriptors represent or depend directly on the quantum-chemically calculated charge distribution in the molecules, and describes the polar interactions between molecules or their chemical reactivity. Definition of MNRI-C in a molecule is given as follows: where the summations are performed over all of the atomic orbitals, of a Carbon atom in a molecule, HOMO denotes the th Atomic Orbital (AO) coefficient on the highest occupied molecular orbital (HOMO) and is the energy of HOMO orbital. The reactivity indices estimate the relative reactivity of the atoms in the molecule for a given series of compounds, and refer to the activation energy of the corresponding chemical reaction.

Maximum electrophilic reaction index for a C atom, MERI-C is another quantum-chemically calculated charge distribution-related reactivity index and its definition is given as follows: where the summations are performed over all atomic orbitals, of a carbon atom in a molecule, LUMO denotes the th AO coefficient on the lowest unoccupied molecular orbital (LUMO), and is the energy of LUMO.

ESP-HA dependent HDSA1/TMSA, HDSA1/TMSA is one of the quantum-chemically calculated descriptors.

HDSA1 is the hydrogen bonding donor ability of a molecule and its definition is given as follows: where is the solvent-accessible surface area of H-bonding donor H atoms.

ESP-WPSA2 Weighted PPSA , WPSA2 is the surface-weighted charged partial-positive charged surface area. This is one of the quantum-chemically calculated CPSA descriptors. Definition of WPSA2 is given as follows: where PPSA2 is the total charge weighted by partial positive surface area and TMSA is the total molecular surface area.

Maximum partial charge, is one of the electrostatics descriptors and reflects characteristics of the charge distribution of the molecule. The empirical partial charges in the molecule are calculated using the approach proposed by Zefirov et al. [27, 28]. This method is based on the Sanderson electronegativity scale, and uses the concept which represents the molecular electronegativity as a geometric mean of atomic electronegativities. Definition of is given as follows: where is the positive atomic partial charges in the molecule.

ESP-HA/HDCA Dependent, HDCA1/TMSA is the hydrogen bonding donor ability of the molecule divided by TMSA. Definition of HDCA1 is given as: where is the solvent-accessible surface area of the H-bonding donor H atoms, selected by threshold charge.

DPSA1 Difference in CPSAs (PPSA1-PNSA1) (Zefirov's PC), DPSA1D in CPSA is the difference between partial positively- and negatively-charged surface areas:

ALPHA and BETA Polarizability is the polarization of a molecule by an external electric field, and is given in terms of the nth order susceptibility tensors of the molecular bulk. The first-order term and the second-order term that contain information about possible inductive interactions in a molecule are referred to as and polarizabilitity, respectively, [29, 30]: where is the permanent dipole moment of the molecule, is the induced dipole moment of the molecule and is the external electric field.

Internal Entropy at 300 K/(Number of atoms in the molecules) where the definition of internal entropy is given aswhere are the frequencies of normal vibrations in the molecule, are the principal moments of inertia of the molecule, is the symmetry number of the molecule, is the Planck’s constant, is the Boltzmann’s constant and is theabsolute temperature (K).

ESP-FNSA3 Fractional PNSA (PNSA3/TMSA) [Quantum-Chemical PC], ESP-FNS3 is the fractional atomic-charge weighted partial negative surface area and its definition is given as where PNSA3 is the total charge-weighted partial negatively charged molecular surface area.

Minimum e-e repulsion for an N atom, Min e-e REP-N where the extreme (maximum or minimum) values of the electron-electron repulsion energy are for a given atomic species (N) in the molecule, calculated as follows: where and are the density matrix elements over atomic basis and are the electron repulsion integrals on the atomic basis . The electron-electron repulsion energy describes the electron repulsion driven processes in the molecule, and may be related to the conformational (rotational, inversional) changes or atomic reactivity in the molecule [31].

Minimum resonance energy for a C-H bond, Min Re. E-CN is a quantum mechanical energy-related descriptor and its definition is given as where is a given atomic species, is the another atomic species, is the density matrix elements over atomic basis , and is the resonance integrals on atomic basis .

Average 1-electron reaction index for a atom, Avg. 1-e RI-C is a charge distribution-related reactivity index and it estimates the relative reactivity of the atoms in the molecule for a given series of compounds, and this corresponds to the activation energy of the corresponding chemical reaction. Definition is given as where is the highest occupied molecular orbital MO coefficients, and is the lowest unoccupied molecular orbital MO coefficients.

Maximum n-n repulsion for a C-N bond, Max n-n Re.-CN is a quantum mechanical energy-related descriptor. This energy describes the nuclear repulsion-driven processes in the molecule and may be related to the conformational (rotational, inversional) changes or atomic reactivity in the molecule [27, 28]. Maximum nuclear repulsion energy between two given atomic species (atoms and ) in the molecule is calculated as follows: where and are the nuclear (core) charges of atoms and , respectively, and is the distance between them.

3.3. Statistical Analysis

Codessa Pro was also used for statistical analysis. This code uses diverse statistical structure property/activity correlation techniques for the analysis of experimental data in combination with the calculated molecular descriptors. Heuristic and Best Multi-Linear Regression methods implemented in Codessa Pro were employed for selecting the “best” regression models [24].

In the study herein, the statistical quality of obtained QSPR models were assessed by the statistical parameters , , F, , , RMSE, and the Y-randomization test. , the squared correlation coefficient, is a measure of the fit of the regression equation. , the “leave one out” (LOO) cross-validated squared correlation coefficient, is an internal validation parameter for a model. The LOO approach involves developing a number of models with one sample omitted at a time. F, the Fisher test value, reflects the ratio of the variance explained by the model and the variance due to the error in the model. Higher values of F-test indicate the significance of the equation. is the standard deviation of the regression. is the predicted square of the correlation coefficient for external validation, and is calculated from the test set by applying the equation developed on the training set. Prediction accuracy of the models is also given with root mean square error (RMSE) values. Reliability of the models is indicated by the Y-randomization test.

3.4. Data Set

The experimental acidic dissociation constants pKa data of some N-base ligands and their formation constant, Kc with Bis(2,2′-bipyridine)platinum(II)2+ ion were taken from [19]. The data set consists of 22 N-base ligands. For pKa QSPR modeling, 5 of the ligands were selected as a test set. For Kc QSPR modeling, 4 compounds were selected as a test set. In the Kc QSPR modeling, compound 11 was selected as an outlier due to the fact that its reactivity is several magnitudes higher than the mean value of the set. While selection of training and test set of the compounds, attention was paid to all of the sets for spanning structural and activity diversity for compounds.

4. Conclusions

In the present study herein, the quantum chemical structural descriptors of 22 N-base ligands including pyridines, pyrimidines, purines, and quinolines have been correlated with their experimental pKa using Codessa Pro methodology. In addition, the formation constant (Kc) of reactions of these ligands with Pt(II) (bpy = 2,2′-bipyridine) ion has also been modeled using the methodology. We have introduced three models each for pKa and Kc. pKa models demonstrated relatively better statistics than the Kc models.

The best obtained pKa Model I is the four-parametric regression equation displaying very good statistical fit and predictive power as evident from its , , , and , , RMSE = 0.51 values. An analysis of descriptors that are involved in the pKa models indicate that reactivity index and charge distribution-related descriptors play major roles in the model acid dissociation constant of N-base ligands.

Obtained Kc models exhibit interesting results in terms of validations parameters. If one analyses these models, and conclusion could be drawn that there is not any relationship between and . An analysis of descriptors that are involved in the Kc models indicate that polarity and internal entropy considered with the Nitrogen-related index of N base ligands have dominant effects on the formation constant of reactions of these ligands with the Pt(II) (bpy = 2,2′-bipyridine) ion.

Acknowledgment

E. Eroğlu thankfully acknowledges the BAP Unit of Akdeniz University for the partial financial support.