Quantum topological molecular similarity (QTMS)

Substituent effects can be gauged by Hammett substituent constants,[1, 5, 47, 148] that is, the change in the pK_a of benzoic acid on substitution in the aromatic ring by the substituent S (Scheme 1):

(12)

Popelier developed a method termed “quantum topological molecular similarity” (QTMS) that has been tested extensively for its capability of reproducing the trends predicted from Hammett constants accurately.[123] The QTMS descriptors have been shown to convey essentially the same information contained in the Hammett substituent constants and may constitute an alternative in QSAR modeling since it is generally easier to compute QTMS descriptors than it is to determine the experimental values of new nontabulated Hammett constants.[124-126, 129, 131, 140, 142]

In the QTMS approach, ρ(r) is sampled at the BCPs, a sampling that entails several advantages such as the speed of the calculation and the emphasis on the chemically relevant bonding regions with only an indirect perturbing influence of the chemically irrelevant nuclear regions. Further, QTMS bypasses the molecular alignment problem, a considerable time-saver compared to methods requiring maximization with respect to the parameter Θ in Eq. (7). We note in passing that unlike Eq. (7), however, QTMS cannot distinguish between enantiomers, which does not constitute a limitation as long as the collected experimental data include those obtained using the correct enantiomer(s). Recapping, the QTMS method has a number of advantages with respect to its applicability in QSAR-type drug design: (i) It is entirely based on the topology of the electron density in the chemically relevant bonding region and not dominated by nuclear maxima and (ii) its is considerably fast given that it does not require molecular prealignment and/or space integration as the electron density is only sampled at a set of point (the BCPs) rather than being considered in its entirety.

The QTMS method starts with the construction of a multidimensional mathematical space using bond properties determined at the BCPs. Each bond property is represented by a Cartesian axis after converting each to a dimensionless relative scale to avoid dimensional nonhomogeneity. The dimensionless Euclidian distance between two molecules in this mathematical space is then a measure of their dissimilarity, the larger the distance the less similar they are.

Examples of bond descriptors include, but are not limited to, the electron density at the BCP (ρ_BCP), the Laplacian of ρ_BCP (▿²ρ_BCP), the bond ellipticity (ε), the three principal curvatures of the electron density evaluated at the BCP (λ₁, λ₂, and λ₃), the kinetic energy density at the BCP (K_BCP), the total energy density at the BCP (H_BCP), and the equilibrium bond length (R_e).

If we suppose the BCP space is constructed from ρ_BCP, ▿²ρ_BCP, and ε, then the distance between molecules A and B is:

(13)

where the sum runs over the i common bonds between the two molecules A and B, and where standardised BCP properties are used, which means that property x is replaced by {[x – mean(x)]/standard deviation of x}, to avoid dimensional inhomogeneity. In this way, the smaller d(A,B) the more similar are A and B with a lower bound of d = 0 when A and B are one and the same molecule.

By construction, QTMS is an electronic structure fingerprinting tool, but it is devoid of a built-in capability to represent molecular size or hydrophobicity (as recognized by Popelier on a number of occasions). This is so since the sampling is of intensive properties (densities) rather than extensive properties (integrated densities). This feature can actually be a strength of QTMS rather than a limitation, because whenever QTMS alone (i.e., without the inclusion of extraneous steric or log P information) is capable of yielding a strongly predictive QSAR model, then this implies that the properties being modeled are determined primarily by electronic structure rather than by steric bulk, molecular size, or solvent partitioning tendencies. Conversely, failure of QTMS alone to predict activity is a diagnostic that steric bulk and/or hydrophobicity are important in determining activity. In the latter cases, QTMS can be used in conjunction with traditional steric bulk and/or hydrophobicity measures to yield powerful predictive combined models.

The exclusive sensitivity of QTMS to electronic structure parallels that of the Hammett σ constants. Both Hammett constants and QTMS must be supplemented with steric bulk and/or hydrophobicity descriptors when these are important in QSAR models. This parallelism between QTMS and Hammett constants lends merit to Popelier's proposal of the eventual replacement of the latter by the former given the relative ease with which QTMS descriptors may be obtained.[130, 132]

A remarkable feature of QTMS is its apparent inherent capability to pinpoint the “active center” automatically. By this, we mean that the experimental Hammett constants (σ)-sequence relative to p-aminobenzoic acid (Table 1) is reproduced if the BCPs of the “active center,” and only the active center (here the bonds in the carboxylic group: O()H, C()O, and C(—)OH), are included in the distance calculation. The inclusion of any additional bond(s) destroys the agreement with the experimental sequence, and the more “irrelevant” bonds are incorporated in the Euclidian space the more the disagreement with experiment[123] (Table 1, especially the last column). As a corollary and by construction, the Euclidian QTMS molecular comparison space need not include (and probably is better not to include) all the common bonds of a series of molecules. In other words, one can include in the QTMS distance calculation only those bonds existing in a common “subgraph” that maximizes the predictability of the model. That is to say, given a well-characterized active center in otherwise unrelated molecules (i.e., molecules differing in their molecular skeletons), the relative activity can be predicted from the distances computed from the properties of the bonds in this active center alone. Figure 4a supports this hypothesis as can be seen from the plot of pK_a values predicted on the basis of the Euclidian distance including only the bonds of the carboxylic groups of 40 unrelated organic acids[129] after a selection of the strongest predictors (from the set of BCP properties) on the basis of a partial least square (PLS) regression procedure. A less stringent test is presented in the plots of Figures 4b and 4c that demonstrate the ability of QTMS to predict the pK_a's of related molecules with accuracy: Substituted anilines in 4b and substituted phenols in 4c. The possibility of chance correlations has been eliminated in this study through the “leave one out” cross validation[129] leading to final statistics that are of higher quality than those obtained from traditional descriptors such as Hammett constants.[129]

Experimental sequence		NH2	OCH3	CH3	H	F	Cl	CN	NO2	# of disagreements/8
QTMS sequences	No. of BCPs
aAdapted from Ref. 123 with permission © 1999 American Chemical Society. bS denotes a substituent attached to the ring.
OH	1	NH2	OCH3	CH3	H	F	Cl	CN	NO2	0
CO	1	NH2	OCH3	CH3	H	F	Cl	CN	NO2	0
COOH	3	NH2	OCH3	CH3	H	F	Cl	CN	NO2	0
CCOOH	4	NH2	OCH3	F	CH3	Cl	H	CN	NO2	4
C6	6	NH2	CN	CH3	OCH3	H	Cl	NO2	F	5
C6H4	10	NH2	CH3	OCH3	CN	H	Cl	NO2	F	6
C6COOH	10	NH2	OCH3	CH3	Cl	CN	H	NO2	F	5
C6H4COOH	14	NH2	OCH3	CH3	Cl	H	CN	NO2	F	5
SC6H4COOHb	15	NH2	CN	H	OCH3	CH3	NO2	F	Cl	7

The versatility of QTMS has been demonstrated by its ability to predict a wide range of physicochemical and pharmacological properties. A few highlights from this work are given next for the sake of illustration.

Examples of predictive QSAR using QTMS only (without incorporating steric or log P information into the model) include the prediction of the proton affinities and the pK_a's of 125 pyridine derivatives,[140] the relative bond dissociation enthalpies of phenolic antioxidants,[131] and the interaction energies of substituted cytosine–guanine (CG) Watson–Crick (WC) base pairs.[149] The predicted properties in this selection of examples are all primarily dependent on electronic (and much less on steric or solvent partitioning) properties. The latter study will now be briefly discussed.

Xue and Popelier (XP)[149] calculated the effect of 37 substituents at position 8 of guanine on the binding energy of the CG WC base pair. The position of the substituent R, attached to C8 of guanine, is indicated in the following structure with the most “predictive bonds” (see below) enclosed in boxes (Scheme 2):

The binding energies of the unsubstituted and substituted dimers were all calculated at the B3LYP/6-311+G(2d,p) level of density functional theory (DFT) with the inclusion of counterpoise basis set superposition error corrections.[150, 151] If ΔE represents the interaction energy of the unsubstituted GC base pair and ΔE_S that of a guanine(C8)-substituted GC pair, then ΔΔE ≡ ΔE_S − ΔE represents the change in ΔE on substitution. Defined in this manner, a positive value of ΔΔE implies a substitution that destabilizes the base pair and a negative value is a substitution that favors the hydrogen bonding of the WC pair. XP then compare a multiple linear regression model with a model constructed solely using Hammett σ_m constants. The QSAR models these authors report are:

(14)

where RMSD is root mean square deviation, using traditional Hammett constants, and:

(15)

using the QTAIM properties,[17] where ε (≡λ₁/λ₂ − 1) is the bond ellipticity and G the gradient kinetic energy density ( ) evaluated at the BCP, and where the bonds are those enclosed in boxes in the chemical structure above.

As can be seen from the regression quality indicator lines below Eq. (15), the molecular set has been divided into to subsets: A training set used to construct the model and a set external to the training set to test the performance of the constructed model outside of its bounds (this is necessary to ensure the predictivity of the QSAR model rather than just its ability to model the given data accurately). Clearly, the QSAR model based on QTAIM descriptors is characterized by significantly superior statistical quality indicators that the corresponding traditional Hammett constants-based model. This example, while not strictly QTMS since Eq. (13) has not been used, is nevertheless closely associated to QTMS since the descriptors used in Eq. (15) are typically incorporated into QTMS models. These results add weight to the proposal that QTMS can be an inexpensive alternative to the Hammett constants-based QSAR approach.

The utility of QTMS has been extended beyond the prediction of static molecular properties. For example, this approach has been shown capable of reproducing the rate constants (k) of a model reaction, namely, the hydrolysis of esters:[127]

The alkaline hydrolysis rate constant at 25°C of 40 different esters (16 ethyl-containing and 24 other esters) were calculated from a QTMS QSAR model and compared with experimental values. The QSAR model yielded r² = 0.930 and q² = 0.863. Further, an examination of the “variable importance in the projection (VIP)” obtained from the PLS regression has (correctly) identified the bonds in molecular fragment enclosed in square brackets, R[(CO)O]R′, as the most important in determining activity,[127] that is, essentially automatically locating the active center.

It is customary to regard the accurate prediction of pK_a's as one of the initial validation steps of a new drug design approach. QTMS has moved beyond that stage and has been applied to the direct prediction of biological properties. For example, QTMS[152] has been shown to have comparable predictivity as the more traditional comparative molecular field analysis (CoMFA)[153] modeling of the toxicological indices of a series of health-hazardous lipid-soluble environmental pollutants (polyhalogenated dibenzo-p-dioxins).[154] A regression of r² = 0.91 and q² = 0.86 was obtained in a QTMS modeling of the cell-growth inhibitory activities of 15 substituted (E)−1-phenylbut-1-ene-3-ones.[155] The method has also delivered the automatic identification of the region in these congeners responsible for their biological activity, namely, the site of the Michael addition.[155] A study of the cytotoxicity of o-alkyl substituted 4-X-phenol (X being the substituent)[133] has shown that this biological activity is accurately predicted by QTMS alone (without the inclusion of any steric information), which indicates that cytotoxicity is primarily governed electronically, a result that does not support a previous proposal by the Hansch group that the latter is more important in determining activity.[156] The reverse situation is encountered in a study of the hepatotoxicity of a series of phenols, which demonstrates that, in this case, QTMS must be combined with external hydrophobicity descriptors such as log P, (see the section entitled “Partitioning Between Immiscible Liquids and log P” below) for a successful modeling,[134] liposolubility being a known key factor in determining the hepatoxicity of phenols. More recently and similarly, it has been shown that QSAR models for the toxicity of aldehydes, another class of environmental pollutants, to a model organism (Tetrahymena pyriformis) can be constructed from a proper mix of QTMS descriptors and log K_o/w[141] [defined in Eq. (42) and discussed in the section entitled “Partitioning Between Immiscible Liquids and log P”]. An earlier study has also shown that log K_o/w and E_LUMO are essential extraneous ingredients in the QTMS modeling of the toxicity of nitroaromatics, another class of pollutants, to Saccharomyces cerevisae.[136]

Bond properties as predictors of spectroscopic transitions and NMR proton chemical shifts

The calculation of UV excitation spectra accurately remains a challenge for computational quantum chemistry to this day and usually necessitate an expensive high level of explicit treatment of electron correlation. The fast empirical modeling of such key analytical characteristics of molecules is, thus, of particular significance from the practical point of view.

The first HK theorem[39] [relation (5)] has been proven for nondegenerate ground states, as discussed above in the section entitled “The Electron Density ρ(r) and QSAR Modeling”. The mapping expressed in expression (5) shows, however, that the ground-state density, ρ(r), does specify completely the Hamiltonian operator (Fig. 2). By specifying and through the time-independent many-particle Schrödinger equation, the ground-state density ρ(r) in principle also determines the excited states and their properties. In other words, the energies of the excited states are also functionals of the ground-state density, even though they are not functionals of the corresponding excited state densities since there exists different external potentials that can yield the same excited state density (i.e., there is no HK theorem for excited states).[157]

It is not entirely surprising, thus, to realize that (approximations to) the ground-state density and ground-state geometry, both of which are mapped to the ground-state external potential, contain information about excited states. In this section, we show that the ground-state bond properties do indeed encode information that can be used in the QSAR modeling of UV excitation energies.

About seven years ago, Buttingsrud, Alsberg, and Åstrand (BAA) have been able to accurately predict λ_max and the corresponding excitation energies (ΔE_hν) of 191 substituted azobenzene dyes from two sets of ground-state QSAR models: The first is based on optimized ground-state bond lengths, and the second on QTAIM BCP descriptors of the ground state.[146] The molecular set studied has the following general structure (Scheme 3):

where some of the R's may belong to condensed cyclic systems.

Several PLS models based only on optimized ground-state bond lengths (including different subsets of bonds in the molecules) have been tested in modeling the excitation properties. Even the simplest model, the one that includes only the NN bond length, yields already an elevated r² of 0.90 (λ_max) and a corresponding r² = 0.86 (ΔE_hν); indicating that this bond appears essential in describing the excitation[146] and hence represents the chromophore.

The cross validated statistics of the one-bond length model are q² = 0.89 for the training set and q²(test) = 0.91 for the test set with a root mean square error of cross validation (RMSECV) of 11.2 nm, and a root mean square error of prediction (RMSEP) of 12.6 nm. The equivalent one-bond (NN) QTAIM bond property model (based on two descriptors: ▿²ρ_BCP and ε) yields q² = 0.91 for the training set and q²(test) = 0.92 for the test set with a RMSECV of 9.9 nm and a RMSEP of 12.3 nm.

The best QSAR of all in both sets of models are those that includes all bonds, in which case the statistics are (bond length model/QTAIM model): q² = 0.94/0.97, q²(test) = 0.96/0.98, RMSECV = 8.0/5.4 nm, and RMSEP = 7.6/5.4 nm. The above observations are consistent with the trends obtained from an inspection of Table 2 of Ref. 146 that (i) generally, QTAIM modeling yields more robust statistics than the modeling based on bond lengths only, and that (ii) there is an improvement, albeit modest, from the simplest QTAIM model (based only on the NN bond) to the most complex based on all bonds, a small gain at the cost of a significant increase in the complexity of the model by the empirical fine-tuning of the perturbation of the chromophore by its surroundings beyond the perturbation captured in the mere change in the NN bond length and properties.

In another study, BAA use QTAIM bond properties to predict 170 distinct NMR proton isotropic shieldings in 60 substituted benzenes obtained from the direct electronic structure calculation. The 60 molecules were grouped into four separate groups, with some overlap (some molecules are members of more than one group), according to the nature of the substituents: group I, II, III, and IV, which include 12, 11, 38, and 28 compounds, respectively.[145] The purpose of the study is not to calculate the NMR shifts themselves since these are calculated directly (and more accurately) by the authors, but rather to use these data as a test-case example of a local property to be predicted by QTAIM bond properties.

Models of increasing complexity ranging from the simplest model, which includes the QTAIM bond properties of only the bond shared by the proton of interest (model a) to the most complex model that includes the properties of all bonds (model e) are tested for their predictivity. The models are further classified in three categories: Those based on bond lengths only (models B); those based on three “standard” QTAIM bond properties, namely, ρ_BCP, ▿²ρ_BCP, and ε, (model S); and finally those that include more QTAIM bond properties, a group of models termed “extended” (model E). An examination of Table 4 of Ref. 145 reveals a general and significant improvement in the S models over the B models, but only a marginal improvement over the S models on passing to the E models. This observation is consistent across the table, hence we quote here only the extreme values for the simplest model a (one-bond only) and the most complex model e (all bonds) applied to group I of 12 molecules exhibiting 20 distinct proton shieldings in the order B/S/E: model Ia: q² = 0.6/0.98/0.98, RMSECV = 0.12/0.027/0.025 ppm, and model Ie: q² = 0.995/0.997/0.996, RMSECV = 0.014/0.011/0.012 ppm.[145]

In some cases, the simplest model a fails completely to capture the proton shielding (group II of chlorobenzenes). The modeling in this case fails presumably due to the size and relatively large electron population of chlorine (compared to fluorine, for example) and which may cause “through space” shielding that invalidates the BCP modeling. Again, as in QTMS, bond properties can be used in electronic fingerprinting but fail to capture size effects. Even in these failed cases, the accuracy of the modeling is quickly recovered in models of the complexity c (seven bonds) and beyond for this group of molecules. For a comparison (B/S/E): model IIa: q² = 0.1/0.5/0.0, RMSECV = 0.11/0.11/0.08 ppm, and model IIe: q² = 0.97/0.96/0.96, RMSECV = 0.019/0.022/0.020 ppm.[145]

Atomic energy and charge of the acidic hydrogen as predictors of pK_a

Adam demonstrated that the pK_a of a weak acid is inversely proportional to the QTAIM energy of the dissociating (acidic) hydrogen in the parent undissociated acid molecule at a given set of thermodynamic conditions.[158] The milestones of his demonstration are summarized here.

The pK_a of an organic acid is related to the Gibbs energy of the acid dissociation reaction ( ) and is expressed:

(16)

Re-expressing ΔG⁰ in terms of the molar dissociation energy (U⁰) and the standard molar partition functions [q⁰; not to be confused with atomic charge, which is given the symbol q(Ω)], we get:

(17)

where

(18)

and where

(19)

The last term in Eq. (17), defined in Eq. (18), depends on the characteristics of H₂O and H₃O⁺ and the temperature and, hence, is independent of the nature of the acid. The assumption implicit in Eq. (19) is expected to be more valid the larger the acid molecule.[158] Substituting approximation (19) into Eq. (17):

(20)

The molecular dissociation energy is given by:

(21)

where E is the bottom-of-the-well energy of dissociation, E^ZPE the zero-point vibrational energy, E_i is the electronic energy of the separated ith ground-state atom. The numerator of the first term of the R.H.S. of Eq. (20) can then be re-expressed as:

(22)

in which and E_HA are the vibrationless total energies of the anion and the undissociated acid, respectively. The difference is expected to be small and approximately constant for a homologous series, and E_H is the energy of free ground-state hydrogen atoms, which is another constant, leading to a re-expression of Eq. (20) as:

(23)

where C is a new constant.

Assuming that the atomic energy of the dissociating hydrogen atom can be approximated as the difference between the total energies of the anion and the neutral undissociated acid, especially for a larger acid molecule:[158]

(24)

then,

(25)

It is important to realize that this result is independent of the particular model chemistry, which also may or may not incorporate solvation effects. The model embodied in Eq. (25) is validated by an excellent agreement with experiment. For instance, the following equation is the result of a statistical fitting of the pK_a's of 53 unrelated organic acids to Eq. (25):[158]

(26)

which yields an r² of 0.991.

The correlation of the experimental and calculated pK_a values and the correlation of pK_a's with the QTAIM atomic energies of the acidic hydrogen are displayed in Figure 5. The statistical correlation indicators reported by Adam[158] based on Eqs. (25) and (26) for different sets of molecules are collected in Table 2.

Adam uses DFT [PW91/6-311+G(d,p)] model chemistry from which he obtains the atomic energies of the acidic hydrogen,[158] raising the questions of (i) the physical meaning of QTAIM atomic energies obtained from DFT calculations and (ii) the effect of using a different level of theory with the possible result of obtaining significantly different atomic energies. The first question is addressed in detail in the Appendix of Ref. 159. As for the second question, QTAIM atomic energies that are numerically obtained by the same integration procedure but from different underlying electronic structure theories (such as Hartree–Fock, Møller–Plesset perturbation theory, or DFT) are very strongly correlated statistically (r² typically = 0.99) and hence, from the point of view of statistical modeling, atomic energies contain essentially the same information no matter what underlying computational level of theory has been used and despite of their possibly significantly differing magnitudes.[32] References 30 and 31 contain additional discussions of the effects of underlying electronic structure calculation levels of theory on QTAIM-derived properties.

Simultaneous to the appearance of Adam's paper, Hollingsworth, Seybold, and Hadad (HSH)[160] report QSAR models to predict the pK_a values of a series of substituted benzoic acids from models based on the atomic charges of the carboxylic group and of the acidic hydrogen, using several different definitions of atomic charges including that of QTAIM. Various data from the HSH paper are assembled in Table 3, which lists the pK_a's against the Hammet constants [defined in Eq. (12)], and the total QTAIM charges of the acidic hydrogens, q(H), and of the carboxylic acid group, q(COOH). As can be seen from the table, Hammett constants are the better predictors in this case (r² = 1.000), which is expected since they are constructed to deliver pK_a's of substituted benzoic acids, but the QTAIM charges of the acidic hydrogens are also excellent predictors of the respective acidities (r² = 0.913):[160]

(27)

which is in line with the simple intuition that an acid with a more electropositive acidic hydrogen (greater positive charge, q(H) > 0) is more acidic and, thus, has a lower pK_a value, as reflected in the negative coefficient of the second term.

	V0 (Exptl.)	V0 (Calc)	V(vdW)	CSI
aExperimental values taken from Ref. 162. V0 values are in cm3 mol−1, van der Waals volume and CSI are in a.u. The calculated V0 are obtained from Eqs. (35) and (36).
Amino acid
gly	43.3	41.9	47.25	0.009
ala	60.5	57.8	212.22	0.221
ser	60.6	62.0	277.93	2.707
cys	73.4	76.3	406.30	0.741
asp(−)	73.8	79.0	473.10	5.083
thr	76.9	77.3	435.30	2.805
asn	78.0	80.5	493.10	5.487
pro	82.8	81.4	461.70	1.070
glu(−)	85.9	93.6	625.23	5.350
val	90.8	87.3	519.70	0.737
gln	93.9	95.2	644.87	5.665
his	98.8	96.2	666.83	6.982
met	105.4	106.9	715.79	0.275
ile	105.8	101.9	670.83	1.010
leu	107.8	102.2	674.48	1.019
lys(+)	108.5	107.8	760.30	4.175
phe	121.5	122.5	879.93	0.698
tyr	123.6	127.3	947.65	2.826
arg(+)	127.3	119.1	926.27	9.679
trp	143.9	146.3	1146.63	3.238
Group
NH	11.6	10.7	126.83	1.826
CO	13.1	12.7	167.78	3.119
NH2	15.4	16.4	175.63	2.013
CH2	15.9	17.4	150.49	0.309
COOH	25.8	25.5	306.26	5.077
CH3	26.5	25.5	213.79	0.277
CH2OH	28.2	27.7	277.93	2.707
CONH2	28.8	29.6	344.49	5.363

There is a striking consistency and similarity in the form of Eqs. (27) and (26), which deserves a comment. Atomic energies are always negative (E(Ω) < 0, always), and they usually become more negative with an increase in the electron population [N(Ω) = Z_Ω − q(Ω)] of a given atom. This is especially true for electropositive hydrogen atoms where the intra-atomic electron–electron repulsion is nonexistent. It is, thus, reasonable that the intercept in Eq. (26) is also negative leading to a lower pK_a for an acid with more electropositive acidic hydrogen since it has a higher (less negative) energy, just as HSH result, Eq. (27), suggests.

QTAIM descriptors of bulk and charge-separation for the prediction of partial molar volumes (PMV)

The partial molar volume (PMV, or V⁰) is the rate of change of the volume of a solution with respect to the number of moles n of solute, extrapolated to infinite dilution. This quantity has the units (L³ mol⁻¹) and for a two-component system is defined:

(28)

where V is the volume of the solution at a given temperature (T), pressure (P), and number (n) of moles of solvent. Equivalently, V⁰ is the first derivative of the chemical potential of the solute with respect to the pressure.[161] Experimentally, the “apparent” PMV ( ) is measured, which only equals V⁰ at infinite dilution (to eliminate solute–solute association) since[162]

(29)

where the concentration of the solute is expressed in molality (m), and thus, the second term of Eq. (29) vanishes only at infinite dilution (when m = 0). All reported experimental PMVs are to be understood, thus, as extrapolated values at infinite dilution.

PMVs are additive and can be estimated from group additivity schemes,[162-164] paralleling the behavior of QTAIM atomic properties expressed in Eq. (11). Further, the observed V⁰ can be considered as the resultant of two principal opposing contributions: A positive intrinsic steric bulk contribution due to the physical space occupied by each solute molecule, which displaces the surrounding solvent molecules ( ), and a generally negative electrostriction contribution ( ) due to the shrinking of the solvent shells surrounding the solute when both the solute and the solvent are polar ( can at times be positive if the solvent–solvent attraction is stronger than the solute–solvent attraction). Ignoring other contributions such as the temperature-dependent solute molecular vibration contribution (positive), we can write[163]

(30)

which can be expressed in terms of QTAIM atomic contributions (after substituting the symbol for approximate equality by that of equality for simplicity):

(31)

leading to the definition of an atomic contribution to the PMV as

(32)

where the atomic intrinsic volume is identified with the QTAIM atomic volume up to the ρ= 0.001 a.u. isodensity envelope, namely, the van der Waals volume , leading to the equivalent definition for a group (or a molecule) R:

(33)

where the sum runs over all the atoms in R.

The atomic electrostriction contribution in water is equated in this modeling to the magnitude of the atomic charge since a charge of either sign will pull water molecules toward the solute (with reversed orientation), ignoring the solvent–solute- and electric charge sign-dependent difference in the magnitude of electrostriction. The justification of this assumption is empirical, that is, the ability of the resultant model(s) to fit experimental data with accuracy. For a molecule or a group of several atoms R, we identify the charge separation index (CSI),[165] namely the sum of the magnitudes of the atomic charges, as its electrostriction term:

(34)

where multiplication of the CSI by a constant of unit magnitude and of dimensions [volume/charge] is implied.

Figure 6 displays the atomic charges on the side chains of a nonpolar amino acid (alanine) and a polar amino amino acid (serine) of roughly similar molecular size. The total charge on the two side chains (R) is +0.08 indicating very little charge transfer to the C_αH_α(NH₂)COOH group. The CSI clearly distinguishes the polar and nonpolar side chain being 0.58 and 2.70 a.u., respectively.

One can construct a QSAR model on the basis of the two contributions described above to model the experimental PMVs of the hydrogen-capped genetically encoded amino acids side chains[162] leading to[80]

(35)

where the subscript AA refers to an entire amino acid and R refers to its side chain.

The agreement between the experimental and calculated values according to this equation can be gleaned from the upper section of Table 4 and upper plot of Figure 7. The agreement between the calculated and experimental values is remarkable especially that only two parameters were used in the fitting of 20 PMVs (with a ratio of data points to parameters = 10:1). This agreement between the QSAR model and experiment extends to the PMVs of the chemical groups composing the amino acids:[80]

(36)

where is the empirical additive group contributions to the PMV,[162] and the subscript G is a short for “group.” The bottom part of Table 4 and bottom plot in Figure 7 illustrate the agreement between the empirical and calculated PMV of the groups composing the naturally occurring amino acids according to model (36).

The coefficients of the last two terms of the R.H.S.' of Eqs. (35) and (36) exhibit the expected signs: Positive for the intrinsic volume terms and negative for the electrostriction terms. There is, however, a subtle but important difference between these two equations: The former has a significant intercept while the latter has practically none. The intercept of about 37.3 cm³ mol⁻¹ in Eq. (35) represents the (predicted) PMV of the group which is common to all amino acids and which has not been directly modeled in this work. The (experimental) contribution of that group to the PMV can be estimated using the additivity property of PMVs and the experimental values listed in Table 4:

(37)

The experimentally estimated value of 34.0 cm³ mol⁻¹ is gratifyingly close to the intercept of Eq. (35). Further, the estimated PMV of the zwitter-ionic group, 37.3 cm³ mol⁻¹, is significantly smaller than that of a hypothetical nonzwitter ionic group estimated to be 52.7 cm³ mol⁻¹ using Eq. (36) and the group values listed in Table 4, and this is also expected due to electrostriction shrinkage accompanying the charge-separated zwitter-ionic form.

Another important class of biomolecules are the five nucleic acid bases: Adenine (A), guanine (G), cytosine (C), thymine (T), and uracil (U). Lee and Chalikian[166] determined accurately the PMVs of the five corresponding nucleosides as well as those of free A, C, T, and U bases at 18 and 55°C (the PMV of the free G base was not determined due to experimental impediments). In the following model, it is shown that the QSAR modeling using again an intrinsic volume term and an electrostriction term is sufficiently sensitive to capture not only the PMVs but also their small changes as a function of a temperature change of 37°.

The contribution of the sugar moiety to the PMV of the nucleosides can be evaluated using the additivity and transferability of groups contributions by assuming:[167]

(38)

The entries in Table 5 demonstrate the simultaneous constancy (transferability) and additivity of the contribution of the sugar moiety to the PMV (80.4 ± 1.0 cm³ mol⁻¹), Eq. (38), a value unaffected by the temperature change from 18 to 55°C within the combined experimental and statistical precision. The empirical PMV of guanine that has not been directly measured can be estimated from the difference between the PMV of guanosine minus the constant contribution of the sugar moiety, that is[167]

(39)

With a complete data set of PMV of all five nucleic acid bases, we can now proceed as before and model their PMV with an intrinsic and an electrostriction term to obtain[167]

(40)

and

(41)

The excellent agreement of the modeling at both temperatures can be seen from the listings in Table 5 and the plots of Figure 8. The increase in the average molecular kinetic energy and population shifts to higher vibrational levels increase the effective molecular volume at the higher temperature, which is reflected in the larger effective PMV at 55°C, an effect well-reproduced by the fitting.

Partitioning between immiscible liquids and log P

A solute in excess at equilibrium with two immiscible liquid phases sharing an interface will result in saturated solutions containing different concentrations of the solute in each of the two phases. In the absence of an excess of solute, the solute will be partitioned between the two phases in the same ratio if the activity coefficient remains approximately constant.[168] The ratio of the solute concentrations in two touching phases in equilibrium, C₁ and C₂, respectively, namely the equilibrium constant K is also known as the distribution ratio, the distribution coefficient, or more commonly as the partition coefficient (P_o/w or simply P). Often the two solvents include an organic hydrophobic/oily phase (often octanol) and an aqueous phase. Thus, the partition coefficient, defined as

(42)

measures the tendency of a compound to be preferentially distributed in the organic or in the aqueous phase. The concentrations (often molarities) rather than activities are used in the definition of P since at equilibrium the activities of the solute in both the organic and the aqueous phase are equal and in this case K would always be unity, which is not useful.

Due to the form of thermodynamic equations, it is more common to use log P in QSAR modeling rather than P itself. The value of log P is smaller than zero when the concentration of the solute is higher in the aqueous phase, that is, in the case of hydrophilic solutes, and conversely it is positive for hydrophobic solutes.

The partition coefficient plays a central role in pharmaceutical sciences and technology as it governs the absorption, distribution, and binding of drugs and determines a diversity of phenomena from the activity of narcotic drugs to the efficiency of solvent extraction of a compound from its solution in another solvent.[168] For these reasons, extensive compilations of log P obtained from experimental measurement and/or estimated from group contributions have long been available.[2] The importance of log P has been underscored, more recently, in the Lipinski's rule of five obeyed by the majority of orally active drugs. The rule states that orally active drugs must satisfy at least three of the following criteria:[169]

• A molecular mass less than 500 D.
• No more than five hydrogen bond donor groups.
• No more than 10 hydrogen bond acceptor groups.
• An octanol–water partition coefficient of less than 5 (log P < 5).

Drugs that significantly violate these rules, except compounds that are substrates for biological transporters, are poorly absorbed and have to be administered by other routes, for example by injection.[169] (The rule acquired its name apparently because all the numbers that appear in it are multiples of five).

Partitioning between two phases across an interface such as a biomembrane is common in the living system. For example, the serum of blood in contact with the neuronal membranes constitute the blood–brain-barrier that must be crossed by any drug intended to act inside the brain tissue. Another example is the necessary ability of a molecule to cross the mucosa of the small intestine before reaching the systemic circulation. This passage through the intestinal mucosa is frequently accomplished by passive diffusion through partitioning of the drug between the aqueous phase in the lumen of the intestine, the oily phase of the mucosa, and finally again an aqueous phase on the other side of the mucosa, that is, the blood serum.

Soluble serum proteins are often exploited as drug reservoirs for the slow and steady release of sparingly soluble drugs, the steady plasma concentration being maintained due to the partitioning with the protein reservoir. The equilibrium constant of binding of a molecules to a protein, K, can be predicted accurately from P by the equation:[2]

(43)

where k₂ and k₄ are dimensionless fitting parameters, which have the values 0.751 ± 0.007 and 2.301 ± 0.15, respectively, when the protein is bovine serum albumin (n = 42 compounds, the standard error = 0.159 log K units, and an r² = 0.960).[11]

A pharmacological endpoint is often expressed at the L.H.S. of equations similar in form to Eq. (43) [Eq. (1)] with a R.H.S. consisting of f(P). The linearity of Eq. (43) is deceptive since it is merely an approximation to a linear tail of a more general parabolic functional relationship between log (1/C) and log P when the range of values is enlarged. The more general relation has the form:[1-3, 5]

(44)

where k₁ and k₃ are fitting parameters, σ is the Hammett substituent constant defined in Eq. (12), and where the second and fourth terms are the R.H.S. of Eq. (43). The third term, k₃σ, is of importance only for molecules that possess a common skeleton with different substituents but looses its relevance when the molecular set used to construct the QSAR model do not share a common skeleton, in which case the equation simplifies to[5]

(45)

The parabolic dependence of log (1/C) on log P in Eq. (45) expresses the requirement that a molecule must possess an optimal log P, referred to as log P⁰,[170] for maximal pharmacological activity—higher or lower log P results in a reduction in activity. The optimal log P ensures that the solubility in the oily phase is just enough to cross biological membrane barriers but not elevated enough to hinder the eventual release of the drug from the membrane to be redistributed on the other (aqueous) side where the site of action is located. Consistent with these remarks is the observation that a log P⁰ ≈ 2.3 seems to be a universal requirement for a molecule to be a general anesthetic,[170] which is the value obtained from the fitting of the activities of ethers for which k₁ = 0.22, k₂ = 1.04, and k₄ = 2.16.[170]

The partition coefficient [Eq. (42)] may be expressed as a ratio of the equilibrium mole fractions (rather than molarities) in which case it is labeled with a prime (P^′). The reason for this modified definition is the appearance of P^′ in the thermodynamic expression of the standard free energy of transfer of a solute from solvent 1 to solvent 2:[2]

(46)

which embodies the experimental procedure of measuring the Gibbs energy of transfer from one solute to the other from partitioning experiments.

Conceptually, the most drastic partitioning between two phases is the partitioning of a solute between its gas-phase and its aqueous solution phase. The Gibbs energy of transfer from the gas-phase to the aqueous medium measures the affinity of a given solute to water, a property crucial in determining proteins folding in water.[171] The molar free energy of transfer from the gas phase to the aqueous phase, termed the molar free energy of hydration or “hydration potential” (ΔG_hydr.), has been accurately determined by Wolfenden et al.[172] at 25°C after correction to pH 7 for each of the hydrogen-capped amino acid side chains ((R); since the amino acids themselves have very small vapor pressures due to their zwitter-ionization in water). These hydration potentials are accurately modeled using a single fitting parameter that measures the hydrophilicity/hydrophobicity, namely, the CSI[80] [defined in Eq. (34)]:

(47)

The experimental and calculated ΔG_hydr. are listed along with CSI_R in Table 6, which is sorted in order of decreasing CSI (hydrophilicity). The correlation of the experimental and calculated values obtained from Eq. (47) is displayed in the top of Figure 9.

AA	Descriptors		Modeled properties
			mRNA 2nd letter	Baseb
	CSIR	V(vdW)R			exptl.	calc.	exptl.	calc.	exptl.	calc.
aCSI and van der Waals volumes are in a.u. and all energies are in kcal/mol. Experimental values were determined by Wolfenden et al.[172] bChemical nature of the base: I = purine, II = pyrimidine. cSerine is the only amino acid with a degenerate middle letter of the triplet code.
arg+	9.679	926.27	G	I	−19.92	−19.28	−14.92	−12.65	−1.32	−0.58
his+	6.982	666.83	A	I	−10.27	−13.42	−4.66	−9.13	0.95	−0.42
gln	5.665	644.87	A	I	−9.38	−10.57	−5.54	−6.70	−0.07	−0.06
asn	5.487	493.10	A	I	−9.68	−10.18	−6.64	−7.39	−0.01	−0.42
glu–	5.350	625.23	A	I	−10.20	−9.88	−6.81	−6.22	−0.79	−0.01
asp–	5.083	473.10	A	I	−10.95	−9.30	−8.72	−6.73		−0.34
lys+	4.175	760.30	A	I	−9.52	−7.33	−5.55	−2.99	0.08	0.73
trp	3.238	1146.63	G	I	−5.88	−5.30	2.33	1.50	2.51	2.08
tyr	2.826	947.65	A	I	−6.11	−4.41	−0.14	0.95	1.63	1.67
thr	2.805	435.30	C	II	−4.88	−4.36	−2.57	−2.52	0.27	0.29
ser	2.707	277.93	G/C	I/IIc	−5.06	−4.15	−3.40	−3.41	0.04	−0.11
pro	1.070	461.70	C	II		−0.59		1.06		0.91
leu	1.019	674.48	U	II	2.28	−0.48	4.92	2.62	1.76	1.51
ile	1.010	670.83	U	II	2.15	−0.46	4.92	2.61	2.04	1.50
cys	0.741	406.30	G	I	−1.24	0.12	1.28	1.33		0.87
val	0.737	519.70	U	II	1.99	0.13	4.04	2.11	1.18	1.18
phe	0.698	879.93	U	II	−0.76	0.21	2.98	4.66	2.09	2.17
met	0.275	715.79	U	II	−1.48	1.13	2.35	4.36	1.32	1.86
ala	0.221	212.22	C	II	1.94	1.25	1.81	1.02	0.52	0.51
gly	0.009	47.25	G	I	2.39	1.71	0.94	0.31	0.00	0.13

It has long been recognized that the second position of the triplet genetic codon of an amino acid is the most important in determining the physical properties of that amino acid.[173] By sorting the amino acids on the basis of the hydration potentials of the side chains, Wolfenden et al.[172] noted that in the mRNA genetic code, hydrophilic amino acids have a purine base in the second position (A or G) while hydrophobic amino acids tend to have a pyrimidine base (U or C) in the second position. Since the hydrophilicity/hydrophobicity of the amino acid side chains is highly correlated to the CSI, one can glean from Table 6, a correspondingly strong correlation between the chemical nature of the middle letter of the triplet genetic codon and the CSI of the side chain of the encoded amino acid. All amino acid side chains with 9.68 > CSI_R ≥ 2.83 a.u., the hydrophilic amino acids according to CSI, possess a purine base in the second position without exceptions. Conversely, most hydrophobic amino acids with CSI_R ≤ 2.81 a.u. possess a pyrimidine base in the second position, with a few exceptions: (i) Serine, which is the only amino acid with a degenerate middle letter in its codons (middle letter can be either “G” or “C”), (ii) glycine, which is an exceptional amino acids in that it does not have a side chain and is the only amino acid with an achiral α-carbon, and (iii) cysteine, which is the only amino acid that dimerizes through its side chain to form SS bridges of cystine. To the author's knowledge, this is probably the most direct link between the genetic code and the electron density.

The experimental Gibbs energies of transfer of the 20 hydrogen-capped amino acid side chains from cyclohexane (cyc) and from octanol (oct) to water (w) were reported by Radzicka and Wolfenden.[174] These experimental can be modeled with an electrostriction term [CSI_R, Eq.(34)] and an intrinsic volume term [V(vdW)_R, Eq.(33)]:[15, 80, 167]

(48)

and

(49)

where energies are in kcal mol⁻¹, and CSI values and volumes are in a.u. Table 6 provides a listing of the experimental and calculated values corresponding to Eqs. (48) and (49) along with a listing of the amino acid descriptors CSI_R and V(vdW)_R used in the modeling, while Figure 9 (middle, and bottom) displays the agreement between calculated and experimental values graphically. A number of other remarkable models were reported where CSI_R and V(vdW)_R could be used in the accurate prediction of the change in the stability of proteins on mutation, that is, the change in Gibbs energy upon denaturation of the mutant compared to this change for the wild type.[80, 167]

The versatility of the CSI index to model widely differing properties is quite notable. Recently, Sánchez-Floresy et al.[175] report accurate modeling of excitation energies in simple molecules using this single descriptor. For example, these workers show that the excitation energies for the low-lying singlet and triplet π→π* states of CO can be closely fitted to the following linear regression equation:[175]

(50)

while transition energies between low-lying singlet and triplet π→π* states of benzene can be fitted to:

(51)

again using the ground-state CSI, which reinforces the arguments presented above in the opening of the section entitled “Bond Properties as Predictors of Spectroscopic Transitions and NMR Proton Chemical Shifts”.

A general conclusion that can be drawn is that the QTAIM ground-state CSI is capable of capturing a good deal of important molecular physics, in the QSAR sense, and is an important and valuable descriptor worthy of note and further exploration.

π-Stacking of antineoplastic agents (Casiopenas^®) and charge transfer from flanking DNA base pairs

Casiopeinas^® are a promising class of antineoplastic cytotoxic agents of uncertain mechanism of action.[176] These complexes are denoted [Cu(NN)(NO)]NO₃ and [Cu(NN)(OO)]NO₃ where NN symbolizes a substituted bipyridine or phenanthroline aromatic ring system, NO an α-aminoacidate or peptide, and OO acetylacetonate or salicylaldehyde. These compounds consist of tertradentate complexes of Cu²⁺ where the metal is chelated to a substituted bipyridine or phenanthroline ring from one side (top moiety in the structure below) and to an α-aminoacidate, a peptide, an acetylacetonate, or a salicylaldehyde moiety from the other side (bottom moiety). The general structure of these molecules can be represented by Scheme 4,

where the light gray part is variable, existing in the phenanthroline congeners but nonexisting in the bipyridine congeners, and where R₁ and R₂ are parts of a closed ring system where the two atoms bonded to the central Cu²⁺ may be both oxygen atoms or one can be an oxygen and the other a nitrogen.

Galindo-Murillo et al.[177] have recently suggested that these complexes can intercalate between DNA base pairs through their aromatic moiety and that the π-stacking interaction is driven by an electron depletion of the planar ligand (the substituted bipyridine or phenanthroline ring) due to the transfer of charge to the metal center which, in turn, drives charge transfer from the flanking DNA bases to the intercalating ligand. Using the integrated QTAIM charges summed over entire molecular fragments, these workers found a simple but strong statistical correlation between the complex stabilization energy of the adenine–Casiopeinas^® complex and the net electron population transferred from adenine to the aromatic ligand:[177]

(52)

which lends numerical support to the charge-transfer assisted π-stacking hypothesis advanced by these workers and, simultaneously, enhance the plausibility of their proposed mechanism of action initiated by stacking intercalation of Casiopeinas^® between DNA base pairs. What is desirable for a future study is, perhaps, to correlate energies of stacking interaction with a direct measure of biological antineoplastic activity.

Before leaving this section, we note that there exists a number of similar empirical correlations of π-stacking interactions modeled with QTAIM (bond) properties, we cite as examples those of Zhikol et al.,[178] Platts and co-workers,[179] and Parthasarathi and Subramanian.[180]

Definition of LDM (ζ), DM, and DWCM

Among the properties that can be associated with a bounded region of space is the number of electrons localized within that region. When this region is an atom in a molecule defined by QTAIM, Ω, the average number of electrons localized within this atom, that is, not shared with any other atom in the molecule, is termed the “localization index” (LI, for short) and is given the symbol Λ(Ω). Since not all electrons are localized, in general Λ(Ω) ≤ N(Ω), the equality is approached at the closed-shell limit as in ionic bonding and is exact only for infinitely separated isolated atoms or ions. For an illustration, in the case of the ionic molecule Li^+0.935F^−0.935 molecule (the superscripts are the atomic charges): [Λ(F) = 9.845] < [N(F) = 9.935] with a difference of N(F) − Λ(F) = 0.090. In contrast, for the covalent F₂, the corresponding values are: [Λ(F) = 8.549] < [N(F) = 9.000] with a significantly larger difference of N(F) − Λ(F) = 0.451. (QCISD/6–311++G(3df,2pd) level of theory).

The difference between N(Ω) and Λ(Ω) accounts for the number of electrons belonging to Ω but shared with other atoms in the molecules. In the case of Li^+0.935F^−0.935, 0.090 e^– are contributed by F and an equal number is shared from the Li with a total of 0.180 e^– shared between the two basins [denoted by δ(Li,F) and which is termed the “delocalization index” (DI), in QTAIM parlance]. In this manner, all electrons in the LiF molecule are accounted for: [Λ(Li) = 1.975] + [Λ(F) = 9.845] + [δ(Li,F) = 0.180] = 12.000 e^–.

Thus, QTAIM defines a DI that counts the number of electrons shared between any two atomic basins in a molecule and which is generally denoted by δ(Ω₁,Ω₂). The formal definition and methods of calculation of the LI and DI are discussed elsewhere.[118, 181-187]

Two important points must be emphasized before moving to examples of applications. The first is that, as mentioned above, LI and DI account for the whereabouts of all electrons in a molecule composed of n atoms [Ω_i, i = 1,2,.n], their general relation to an atomic electron population being:

(53)

Given expression (53), the total molecular electron population can be expressed as the sum of two subpopulations:

(54)

where

(55)

and

(56)

The second point, which directly follows from Eq. (53), is that the DI measures the number of electrons (and not the number of pairs of electrons), an erroroneous characterization that has propagated in many publications (including some by the present author). A source of the confusion may be that when two atoms Ω_i and Ω_j are bonded the number of electrons shared between them, δ(Ω_i,Ω_j), is numerically close to the number of shared pairs in the Lewis bonded structure and hence can be misconstrued as a bond order. One can assert, though, that δ is generally proportional to the classical bond order when two atoms share a bond path. A simple example helps settle this confusion. For the H₂ molecule, δ(H,H^′) ≈ 1.0 and Λ(H) = Λ(H^′) ≈ 0.5, which means that N_loc = 2 × 0.5, leaving only 1 e to be delocalized (shared).

It is also worth reminding the reader that δ(Ω_i,Ω_j) is not only defined between bonded atoms but also between any two atoms in a molecule, no matter how distant and regardless of the presence or absence of a bond path linkng their nuclei. Obviously the DI cannot be related to any bond order if the atoms do not share a bond path, that is, when they are not bonded in the first place.

Ferro-Costas, Vila, and Mosquera (F-CVM)[188] have recently demonstrated that the anomeric effect in halogen-substituted methanols cannot be explained by hyperconjugation arguments based on the behavior of atomic populations and the QTAIM localization/delocalization indices. These authors have presented lower-triangular matrix-like tabulations of the delocalization indices and have used differences between these matrices in their argumentation.[188]

Given that the δ indices were found to be excellent predictors of experimental quantities such as NMR JJ coupling constants between protons[185] and fluorine atoms,[187] it is tempting to construct matrices similar to those of F-CVM and use them as tools in QSAR-type studies. Initial explorations of these ideas are presented here for the first time, meanwhile more extensive studies in this direction are being undertook in the author's laboratory.

Disregarding the nature of the matrix elements, the form of the matrices used by (F-CVM)[188] is identical to the oft-used graph-theoretical connectivity (adjacency) matrices. Connectivity matrices have been brought from the field of graph theory, a branch of pure mathematics, to mathematical chemistry by scientists such as Balaban, Gutman, Hosoya, Nicolić, Randić, Trinajstić, and several others.[6, 7, 51, 189-193] In a typical graph-theoretic connectivity matrix, one places 1 or 0 depending on whether two atoms are bonded or not in the hydrogen-suppressed graph. Atom-connectivity matrices have been proposed that do not suppress the hydrogen atoms and that include as off-diagonal entries properties that depend on atomic pairs such as interatomic distances, vibrational force constants, or the bond dissociation energies while atomic symbols are introduced along the diagonal elements.[51]

A natural extension of the graph-theoretical approach is to merge it with the topological analysis of the electron density. The idea does not appear to be new, already in 1981 Dmitriev, in a 155 page book on chemical graph theory, devotes pp. 110–115 to a section on Bader's “Topology of Molecular Charge Distributions,” but makes no connections between it and the main theme of the book.[193] In 1994, Balasubramanian has suggested the “[i]ntegration of graph theory and quantum chemistry”[194] in a paper that concludes with a review of the basic concepts of the topological analysis of the electron density according to the theory now known as QTAIM. In his paper, Balasubramanian has emphasized (i) the remarkable correspondence between the graph-theoretical graph of a molecule and the set of bond paths constituting its QTAIM molecular graph[17] and (ii) the role of the Laplacian of the electron density, ▿²ρ(r), in determining the sites of nucleophilic and electrophilic attack in a molecule. The present author is unaware of other proposals to merge or integrate QTAIM and the graph-theoretical approach and presents below some thoughts on how this might possibly be achieved in the future.

Balasubramanian's integration is intuitively implemented through the construction and the use of “electron density weighted adjacency matrices (EDWAM)” (quoting the words and the idea of Massa (L. Massa, Personal communication, 2014)). An EDWAM places the value of the electron density at the BCP, ρ_BCP, as the off-diagonal matrix element whenever a bond path exists between two atoms instead of the “1” of the graph-theoretic connectivity matrices. There is no hydrogen suppression in the proposed EDWAM. Figure 10 displays a EDWAM for chloroacetic acid. The respective research groups of Massa and of the present author are jointly exploring the use of EDWAMs in QSAR-type studies at the time of writing.

We now introduce a new matrix representation of molecules that encapsulates one- and two-electron information simultaneously and which will be referred to as “localization–delocalization matrix (LDM).” An early version of the LDM was introduced in 2001 in the documentation manual of the program AIMDELOC[195] and bears some similarity to the matrices used here and those used recently by F-CVM.[188]

The diagonal elements of LDMs proposed here are the atomic localization indices Λ(Ω), as originally proposed in 2001,[195] but the off-diagonal elements are the delocalization indices divided by 2, that is, δ(Ω_i,Ω_j)/2, a feature distinguishing the LDMs in this work from the 2001 proposal and from that of F-CVM. Further, we preserve all the elements of these symmetric matrices (not just the lower or upper triangular part). Thus, we define the LDM that we denote by ζ as

(57)

In this manner, and according to Eq. (53), the sum of any row or column of ζ equals to the given atomic electron population N(Ω), while the sum of the column sums or row sums yields the total number of electrons in the molecule [Eq. (54)]. Further, the trace of ζ is the total number of localized electrons (N_loc) in the molecule [Eq. (55)], and the number of delocalized electrons is, then, given by difference:

(58)

As known, δ(Ω_i,Ω_j) decays exponentially with internuclear separation,[196] which implies that the ζ matrix, by including all the delocalization indices in a molecule, does implicitly contain interatomic distance matrix information. The ζ matrix is, thus, rich in coded information about not only the electron density distribution but also about the pair density and even contains some geometric distance information. For these reasons, this matrix can be expected to constitute the raw material for a new class of QSAR correlations.

All matrix representations of molecules have common and well-known drawbacks, which we briefly discuss below with the particular reference to LDMs.

Problem 1: Nonuniqueness due to ambiguity of atom labelling (i.e., a family of matrices related by permutations all code for the same structure)

The matrices are not unique since they are labelling-dependent and an interchange of a pair of atomic labels will result in the permutation of a pair of rows and the simultaneous permutation of a pair of columns. There are n! ways to label n atoms and correspondingly, there are as many different matrices describing the same system. This problem has prompted the graph theoretic community to search for “matrix invariants,” that is, properties associated with these matrices that are labelling-independent. Common invariants include the characteristic polynomial and its roots (the eigenvalues), the trace of the matrix, the determinant of the matrix and so forth. The labelling problem is of course independent of the nature of the matrix elements and hence is also relevant to the LDMs proposed here. We will hence follow the lead of the chemical graph-theoretic community and extract invariants from the LDMs at the cost of complicating the physical interpretation of the final (invariant) descriptors, which are, nevertheless, derived from well-characterized measures of electron localization and delocalization compactly collected in the LDMs. Since an LDM is by construction real and symmetric, then it is diagonalizable to a similar matrix D (the “diagonalized localization–delocalization matrix,” or DLDM) through a similarity transformation:

(59)

The matrix D is indispensible when the atomic numbering scheme cannot be made consistent across a molecular series that lacks, for example, a common skeleton. In these instances, the original ζ matrices of different molecules cannot be compared directly. At times, however, the ζ matrices themselves can be compared directly when a consistent labelling scheme can be imposed on all molecules (e.g., molecules with the same skeleton as in the example of substituted acetic acids discussed in the section entitled “The LDM ζ and the pK_a of Halogenated Acetic Acids” below). Matrices comparison can be accomplished, for example, by calculating the distance between them (e.g., the Frobenius distance, see below).

Problem 2: Unequal matrix dimensions for differently sized molecules

Unless the number of atoms in the compared molecules is the same, the matrix dimensions of ζ will differ from one molecule to another rendering the direct comparison such as calculating the Frobenius distance impossible. A solution to this problem has been proposed by White and Wilson[197] (WW) in the context of image statistical pattern recognition where different sized (sub)graphs are treated on equal footing. The approach of WW consists in determining first the matrix representative of the largest graph (n×n) and enlarge all other smaller matrices to n×n by filling all new empty matrix elements with zeros after placing the smaller matrix in the upper left block of the enlarged matrix, that is, “padding” all matrices other than the largest with zeros to the same common n×n size.[197]

The WW approach is applicable to LDMs since a zero diagonal element means that there are zero electrons localized within a given atomic basin while the zeros in the other (off-diagonal) matrix elements of the row and column that correspond to that atom imply that there are zero electrons shared between that basin and all other atoms in the molecule. This is what computational chemists refer to as a “ghost atom”. Thus, if one wishes to compare the first four saturated aliphatic hydrocarbons, the matrix size will be that of butane (14×14), and the LDMs of methane, ethane, and propane will be enlarged to that size by “zeros padding”. At the DFT-B3LYP/6–311G(d,p) level of theory, and given the atomic numbering schemes in the Supporting Information, to two decimal places, an LDM of butane is:

(60)

while a corresponding (ghost atom-padded) LDM of ethane is:

(61)

Problem 3: The mapping of the set of molecular graphs to the set of adjacency matrix representative can be surjective not injective (nonbijective), i.e., one and the same matrix invariant can code for more than one distinct graph

This problem is illustrated by the example recognized by Živković (Quoted through Ref. 192 due to unavailability of the original reference: T. Živković, Report on Summer School, Repino, Leningrad, 1973) and whereby two chemically distinct molecules (Scheme 5)

coded by their respective adjacency matrices yield one and the same characteristic polynomial (x¹⁰ – 10x⁸ + 33x⁶ – 44x⁴ + 24x² – 4). Naturally, this problem will never occur when the molecules are coded by their respective LDMs or electron density-weighted adjacency matrices with matrix elements that are very different than simple ones and zeros.

The LDM ζ and transferability in aliphatic alkanes (methane, ethane, propane, and butane)

The LDM matrices for the first four members of the aliphatic alkanes homologous series have been obtained at the B3LYP/6-311G** (the complete set of the ζ matrices and their eigenvalues are available in the Supporting Information). We note in passing that while the meaning of the LI and DI is strictly defined within Hartree–Fock theory, the indices obtained from DFT calculations are often used as they yield numerical values of LIs and DIs that are close to the Hartree–Fock values.

First, we investigate one way by which the similarity or dissimilarity of these four molecules can be quantified using their DLDM representations. It is possible to define a “distance” between two matrices (of the same dimensions), even though such a definition is not unique. The definition used in this work is the Frobenius distance, that is, the “Frobenius norm” of the difference matrix:

(62)

There are n(n − 1)/2 different distances d(A,B) between n molecules that may be arranged into a lower or upper triangular intermolecular distance matrix. The Frobenius distance matrix between the DLDMs (D) representing the first four members of the saturated alkane series is given in Table 7.

An examination of Table 7 reveals that the “contribution of a methylene group” to the distance of one member compared to the next in the homologous series follows expectations in the sense that the first addition of a CH₂ to CH₄ to obtain C₂H₆ results in a somewhat anomalous change compared to the subsequent additions of CH₂ to the growing alkane chain, and where this change appears to quickly reach a constant value (as emphasized by Bader in relation, for example, to the contribution of a methylene group to the heats of formation of aliphatic alkanes).[17] Thus, to two decimals, d(CH₄,C₂H₆) = 3.29, d(C₂H₆,C₃H₈) = 3.04, and d(C₃H₈,C₄H₁₀) = 2.92, fast converging to an apparent asymptotic value around 3.

Since this is a relatively uncharted territory, it may be of interest to explore the geometry of this 14-D comparison space. One criterion worthy of consideration is the satisfaction of the triangle inequality as emphasized by Muskulus.[198] The inequality applied to intermolecular metrics of three molecules A, B, and C is:

(63)

A violation of this inequality indicates that the geometric representation is inadequate as a proximity measure of the studied objects[198] and hence would invalidate the use of d(X,Y) as a similarity or dissimilarity criterion of molecules X and Y. Figure 11 displays the triangular distance and angular relationships between the three possible triads of molecules in this group of four aliphatic alkanes. The figure demonstrates that the inequality is clearly satisfied given the intermolecular distance matrix in Table 7. Further, the angular relationships displayed in the figure exhibit a curious regularity: The angles of the triangles formed by three consecutive members of the homologous series are almost identical (compare the first and third triangles). The almost constancy of the obtuse angle ∼99°–100° in all triangles cannot escape notice either. The significance of these observations is not evident to the author at the time of writing.

Another LDM invariant, other than the DLDM representation discussed above, is the trace tr(ζ) = N_loc, Eq. (58). Since in most molecules, the majority of electrons are localized within their respective atomic basins, for a given homologous series N_loc is roughly proportional to the total number of electrons N and is thus an extensive variable that grows with size. In other words, tr(ζ) can be expected to capture molecular properties that depend on steric bulk and/or lipid solubility. Figure 12 displays the correlation of tr(ζ) of the first four aliphatic hydrocarbons with a typical size-extensive property, the calculated total energy E, yielding a correlation with an r² = 1.000 (top), and with a lipid solubility, measured by log P, with an r² = 0.994 (bottom). These preliminary results call for further investigation with many more data points to provide definitive conclusions though.

The LDM ζ and the pK_a of halogenated acetic acids

The LDMs of a set of six homohalogenated acetic acid (CH₂FCOOH, CHF₂COOH, CF₃COOH, CH₂ClCOOH, CHCl₂COOH, and CCl₃COOH) are compared with that of the parent unsubstituted compound, acetic acid. (LDM ζ matrices, their eigenvalues, and the atomic numbering scheme for all members of this molecular set are available in the Supporting Information). In this example, a consistent atomic numbering scheme can be adopted for all seven molecules, all of which have a common skeleton and the same number of atoms, obviating the need to diagonalize their LDMs or pad any matrices with ghost atoms, thus here the distances are calculated directly between the LDMs.

The Frobenius intermolecular distance matrix is given in Table 8. In this table, the first column, for example, lists the distance between every substituted acid from the parent unsubstituted acetic acid. A plot of these distances against the corresponding experimental pK_a's is displayed in Figure 13. The figure shows that the experimental pK_a's of all congeners substituted with one and the same halogen (either F or Cl) are almost perfectly linearly correlated with the Frobenius distances of their respective LDMs from that of unsubstituted acetic acid (r² = 0.98 for both regression lines). The slopes of the two lines are different, however, which signals that the pK_a is not captured in its entirety by the LDM intermolecular Frobenius distance matrix.

CHiX3-i	CH3	CH2F	CHF2	CF3	CH2Cl	CHCl2
CH2F	4.52	0
CHF2	7.78	4.76	0
CF3	11.16	7.89	5.05	0
CH2Cl	9.29	7.47	8.67	10.65	0
CHCl2	14.05	11.91	10.47	11.40	9.17	0
CCl3	19.56	15.83	14.03	12.72	13.87	9.17

The LDM is generally “biased” somewhat by its diagonal elements Λ(Ω_i), which typically have larger magnitudes than the off-diagonal ½δ(Ω_i,Ω_j). In a homologous series, both N_loc [Eq. (55)] and N_deloc [Eq. (56)] increase roughly in proportion to the total number of electrons N [Eq. (54)]. The situation is different for a series of congeners molecules RX differing only in the identity of a substituent, say with X = F, Cl, Br, or I. In this case, the traces of the LDMs will be significantly different since these substituents will contribute a large fraction of their atomic numbers of electrons to the count yet the delocalization indices in these congeners will, at most, be slightly perturbed by the substitution. For these four congeners RX, then, the diagonal elements are the sole carrier of size-extensive information in their respective LDMs. For a property primarily driven by electronic structure, such LDMs, biased by steric bulk, may not work just as we have seen in the case of QTMS (see the section entitled “Quantum Topological Molecular Similarity (QTMS)” above). The inability to fit the two sets of congeners of substituted acetic acids (X_n = F, Cl) on one line using the distances between their LDMs proves the point (Fig. 13).

We thus, define a size-independent electronic fingerprinting matrix representation of a molecule as the diagonal-suppressed-LDM or simply the “delocalization matrix” (DM) in which all diagonal elements are suppressed to zeros (no explicit localization information). The similarity can then be measured as the Frobenius distance between DMs when the modeled properties are primarily governed by electronic structure and not by steric bulk.

The experimental pK_a's are plotted in Figure 14 against the distances between the DMs of the substituted and unsubstituted acetic acids, denoted as d_deloc(AA,SAA), and are listed in Table 9. This size-independent electronic fingerprinting distance correlates much better with the experimental pK_a's than the size dependent LDMs-based distances bringing both sets of congeners (X_n = F, Cl) together on a single exponential regression line (r² = 0.979). A similar conclusion is reached within the QTMS framework that models pK_a's accurately on its own without injecting extraneous size-dependent information in the model, since pK_a is inherently governed by electronics rather than by steric bulk. It may be concluded that the full LDMs appear to capture properties that depend simultaneously on the electronic structure and size (steric bulk) while the DMs are better tuned to capture the properties primarily dependent on electronic structure.

CHiX3-i	pKaa	ddeloc(CH3, CHiX3-i)
aExperimental values obtained from Ref. 200.
CH3	4.756	0.0000
CH2Cl	2.87	0.1000
CH2F	2.59	0.1122
CHCl2	1.35	0.1814
CHF2	1.33	0.2229
CCl3	0.66	0.2404
CF3	0.52	0.3470

We return to Table 8 for some final observations. The distances between CH_iF_3–iCOOH and CH_i−1F_2–iCOOH (i = 3, 2, 1), that is, two homologous members differing by one fluorine atom, are 4.52, 4.76, and 5.05 for i = 3, 2, and 1, respectively. A fluorine atom's contribution to the distance is between ∼ 4.5 and 5.1 and it changes by roughly a constant amount of ∼0.25 with i. The corresponding distances between CH_iCl_3–iCOOH and CH_i−1Cl_2–iCOOH in the case of chlorine substitutions are 9.29, 9.17, and 9.17, exhibiting an almost constant contribution to the distance for every addition of a chorine atom. The table exhibits other regularities that we skip for brevity.