Anatomy of the substituent effect in complex N-heterocycles

Abstract

The substituent effect in heterocyclic compounds is a critical concept in chemistry, influencing the properties and reactivity of these molecules. In simple systems, like benzene derivatives, it is well-understood. However, in heterocycles with multiple substituents and fused rings—important nature's building blocks—the situation becomes complex and hard to comprehend. In this work, quantum-chemical calculations (at the DSD-PBEP86-D3BJ/def2-TZVPP level of theory), including substituent effect descriptors and the electron density of delocalized bonds decomposition scheme, were performed on derivatives of cytosine, isocytosine, guanine, isoguanine, thioguanine, hypoxanthine and 5-aza-7-deazaguanine in different tautomeric forms, allowing the determination and quantitative comparison of the strength of different interactions between substituents and heteroatomic fragments in each studied molecule. The comparisons allowed formulating several rules that help understand and predict the strength, and the resonance/inductive nature of the substituent effect in complex N-heterocyclic systems.

---

Introduction

N-heterocycles, molecules containing at least one nitrogen atom in their ring structures, constitute a significant and important class of organic compounds. Due to their unique chemical and biological properties, these compounds have a wide range of applications in pharmaceuticals, agrochemicals, synthetic organic chemistry and materials science.^1–5 For instance, they constitute 82% of FDA-approved pharmaceuticals between 2013 and 2023,⁶ and over 50% of all approved small-molecule drugs are N-heterocycles.⁷ Many include substituents or heteroatoms, among the most popular are F, S, Cl, and CN.⁶ Currently, machine learning models are increasingly used to predict and optimize the reactivity, selectivity, and bioactivity of N-heterocycles.⁸ These models use different descriptors and computational techniques to achieve high accuracy and interpretability, although challenges related to data quality and experimental design remain.⁹ For this purpose, a deep understanding and qualitative description of the substituent effect (SE) is extremely helpful. The SE is a fundamental concept in organic chemistry. Its understanding allows for the targeted modification of molecular properties such as reactivity, optical properties, intermolecular motif formation, ligand binding, biological activity, and many others.^10–13 The most commonly used characteristics of the substituents are the σ constants derived by Hammett, which provide quantitative information on the electron-withdrawing (EW) or electron-donating (ED) properties of the substituents.^14,15 The classical approach to the substituent effect includes a series of X–R–Y derivatives, in which one can distinguish: (i) the reaction site (fixed group) Y, (ii) the set of changing substituents X and (iii) the transmitting moiety (spacer) R.^16,17 The effect of X on Y is a usually studied aspect. For this purpose, the Hammett equation is used: P(X) = ρ·σ(X) + constant, where ρ, called the reaction constant, describes the sensitivity of the property P to the substituent effect, characterized by σ(X), under given conditions (e.g. solvation, temperature). However, the substituent (X) properties differ when R and/or Y is changed. In the Hammett model, this is represented by the constants σ_m, σ_p for meta and para relationships between X and Y, and by σ⁺, σ⁻ for positively and negatively charged Y, respectively. It follows that the properties of the substituent X depend on Y and R—the so-called reverse substituent effect.^18,19 It should be emphasized that the substituent constants were obtained from dissociation constants of substituted benzene derivatives (i.e. benzoic acids, phenols, and anilines). In other systems, they are used only by analogy. Particularly, in complex heterocycles, where we encounter the lone electron pairs of heteroatoms conjugated with the rest of the π-bond system, the situation differs significantly from that in the model benzene derivatives. In such systems, the quantum-chemical descriptors allow studying the effect of substituents X on each of the exocyclic group. In this study, we use the cSAR (charge of the substituent active region)^20,21 parameter, which quantifies charge accumulation on donating (positive partial charge) and withdrawing (negative partial charge) substituents. Moreover, the cSAR can be used to describe the properties of both the reaction site and the substituents, so the slope of the linear equation cSAR(Y) vs. cSAR(X) characterizes the strength of the substituent effect, as shown in Fig. 1a. To determine the π-electronic structure of heterocycles and visualize resonance interactions associated with the substituent effect, we use a recent π-resonance analysis function defined within electron density of delocalized bonds (EDDB) formalism^22–24—the π-bond delocalization function (BDF_π), explained in Fig. 1b.²⁵

Fig. 1 Introduction of computational methods used in this study: (a) determination of the substituent effect strength with cSAR(X) vs. cSAR(Y) correlations (in the presented case, Y = NH₂); (b) definition of π-bond delocalization function (BDF_π), BDF_π is the difference between delocalized and localized π-electron densities from the EDDB (electron density of delocalized bonds) density partitioning method; and (c) definition of ΔBDF_π, it is a difference between BDF_π functions for substituted and unsubstituted molecules, it isolates changes in delocalization resulting from substitution.

Research objects are the derivatives of nucleic acid bases and their analogues: guanine (G), isoguanine (iG), thioguanine (tG), hypoxanthine (Hx), 5-aza-7-deazaguanine (adG), cytosine (Ct), isocytosine (iCt) and 6-aminopyridin-2-one (6AP) in their most stable tautomeric forms, which are presented in Scheme S1 and discussed further in the manuscript; the study includes a wide range of substituents (X = NO₂, CN, Br, Cl, F, H, SH, OH, and NH₂) with properties ranging from electron-withdrawing to electron-donating, attached at the C8 position of purine bases and C5 or C6 positions of pyrimidine bases. The naming scheme includes the acronym for the heterocycle type, substituent position (omitted for purine bases, only C8–X included) and the tautomeric hydrogen positions (according to the atom numbering repeated in each figure). The goal is to characterize how small, incremental changes in a N-heterocycle structure influence the strength of the substituent effect and identify the structural features responsible for the differences between the systems.

Results and discussion

The substituent effect is transmitted by induction, which acts in the σ-cloud and resonance, acting in the π-cloud. The strength of the inductive effect depends on the groups involved and the distance between them.^16,17 Therefore, between heterocycles substituted at the same position, we can expect similar strength of the inductive effect, but, comparing C5–X, C6–X and C8–X substitutions the differences in the inductive effect will play a role. On the other hand, the π-cloud is often completely reorganized when changing the tautomeric form. This reorganization changes positions of double bonds and lone pairs and can drastically influence the resonance effect, by disrupting existing conjugation paths between the groups or by opening new possibilities for conjugation.

Resonance interactions within N-heterocycles

First, let us introduce different resonance interactions between X, heteroatoms and exocyclic groups which we expect in complex N-heterocycles (Fig. 2). Attaching the substituent introduces additional zwitterionic resonance structures (Fig. 2) into the superposition of resonance forms. These zwitterionic forms are typically high in energy, thus their contribution to the wavefunction should be slight. Nevertheless, the greater their contribution, the stronger the resonance substituent effect is expected. As shown in Fig. 2, the electron withdrawing group (EWG) can accept charge from the exocyclic donating group, or from the endocyclic N(H). Electron-donating groups donate their charge onto exocyclic C=O groups or endocylic N atoms.

Fig. 2 Four types of charge-transfer resonance interactions involving electron donating and withdrawing substituents. Curved arrows indicate the direction in which negative charge is transferred; ortho, meta and para indicate the relative position of fragments necessary for the resonance interaction to occur. The red or blue shaded circle under heteroatom/group indicates that it is conjugated with withdrawing (red) or donating (blue) substituent X, i.e. it can transfer charge onto the substituent becoming positively (red) or negatively (blue) charged, as shown by resonance structures. Charge transfer associated with the presented resonance structures is confirmed by the changes in π-electron populations from NPA analysis and by visualization of isosurfaces (isovalue = ±0.005) of ΔBDF_π (Fig. 1c). ΔBDF_π changes in the π-electronic structure upon attaching the X group—localization (ΔBDF_π < 0, yellow) and delocalization (ΔBDF_π > 0, turquoise).

Simple representations with resonance structures are corroborated by wavefunction analysis: ΔBDF_π results and π-electron populations from natural population analysis (NPA). ΔBDF_π shows that substitution increases the delocalized character of π-electrons along the conjugated pathway connecting the substituent and the interacting group. π-NPA results show an increase or a decrease in π-electron populations at interacting atoms upon substitution, compatible with what the resonance structures show.

How structural changes influence the substituent effect

Fig. 3a–c shows that changing the tautomeric form can strongly alter the SE. When tautomerization disrupts conjugation between X and the NH₂ group in the ortho or para relationship, the SE decreases by ∼60%, accompanied by reduced charge accumulation at the electron-withdrawing CN substituent. This conjugation in the 3H tautomer and its lack in 1H might contribute to the relative stabilization of the 3H form (relative Gibbs energies of 3H in kcal mol⁻¹: +3.6 for NO₂, +4.8 for CN vs. +6.0–12.1 for non π-withdrawing X groups, Table S1).

Fig. 3 Properties (cSAR(X) values in boxes) of X = CN (red), Cl (green), NH₂ (blue) and how they change due to small structural changes noted above the black curved arrows. Intensity of box shading is proportional to the strength of withdrawal or donation of electrons by X. Green and blue curved arrows represent the strength of the substituent effect of X substituent series on exocyclic groups (Y): =O/OH(=S/SH) and NH₂/NH, which was evaluated from the slopes of corresponding cSAR(Y) vs. cSAR(X) linear correlations (X = NO₂, CN, Br, Cl, F, H, SH, OH, and NH₂), as illustrated in Fig. 1a. Acronyms: guanine (G), isoguanine (iG), thioguanine (tG), hypoxanthine (Hx), 5-aza-7-deazaguanine (adG), cytosine (Ct), isocytosine (iCt) and 6-aminopyridin-2-one (6AP).

In studied derivatives of purine bases, the SE is transmitted across two fused rings. In Hx, G and tG systems, breaking conjugation between X and OH/SH, for example by changing tautomer from 9H (Fig. 3d) to 7H (Fig. 3e), reduces the SE by about 10%. In G and tG, the same tautomeric change additionally creates the conjugation between X and NH₂—in non-conjugated 9H the SE between these groups is weaker by 27%. This exemplifies how π-electronic structure reorganization during tautomerization affects the interactions between the functional groups.

Fig. 3d illustrates how adding an NH₂ group in the second ring (transition from Hx 9H,C6–OH to G 9H,C6–OH) strengthens the electron-withdrawing effect of X = CN but weakens the electron-donating effect of X = NH₂ as represented by the cSAR(X) values. These changes are much more pronounced in the 9H tautomer (Fig. 3d) than in 7H (Fig. 3e), because only in 9H is the X → C2–NH₂ pathway conjugated, which is also reflected in the SE strength (−0.177 in G 9H,C6–OH vs. −0.241 in G 7H,C6–OH). This confirms the presence of long-distance resonance interactions between C8–X and C2–NH₂. Disrupting this interaction weakens the SE, as illustrated by transitions tG 1H,7H → tG 1H,9H (Fig. 3h) and G 1H,7H → G 1H,9H (Fig. 3f). The latter transition also breaks the C8–X → C6=O conjugation but creates the C8–X → C2–NH₂ one. Both effects are reflected in the SE strength—C8–X → C6=O is stronger in G 1H,7H (−0.187 vs. −0.168), while C8–X → C2–NH₂ is stronger in G 1H,9H (−0.209 vs. −0.140).

The properties of the C8–X groups change accordingly—electron-withdrawing C8–CN accumulates more negative charge in G 1H,9H while electron-donating C8–NH₂ more positive charge in G 1H,7H. Notably, C8–X → C6=O and C8–X → C2–NH₂ resonance effects operate despite a complicated conjugation pathway between the groups, which in both cases passes through the central C4–C5 bond.

If positions of the NH₂ and =O groups are swapped (G 1H,9H → iG 1H,9H transition in Fig. 3f), the SE of C8–X on C6–NH₂ is enhanced by 50%, compared to C2–NH₂ in G 1H,9H. Accordingly, the distance between X and C=O in this transition increases, and the SE of X on C=O weakens (from −0.168 in C6=O to −0.129 in C2=O, Fig. 3f). These differences are caused by weakening of both inductive and resonance effects with distance.

Resonance substituent effect visualized with ΔBDF_π

Bond delocalization function (BDF_π, Fig. 1b) is a recently introduced approach which allows visualizing the regions of π-electron localization and delocalization, based on the first-principles EDDB electron density partitioning scheme.²⁵ The BDF_π maps in Fig. 4 illustrate the π-electronic structure of unsubstituted heterocycles, while the ΔBDF_π maps illustrate the changes in BDF_π due to the substituent effect—changes in the π-electronic structure associated with the introduction of charge-transfer resonance structures involving X substituents.

Fig. 4 Selected structural changes from Fig. 3 analyzed in terms of the π-electronic structure. BDF_π visualization represents the π-electronic structure of unsubstituted system, while ΔBDF_π shows how it is changed due to CN or NH₂ substitution. Detailed explanation of BDF_π and ΔBDF_π is presented in Fig. 1.

In Fig. 4a, ΔBDF_π reveals the strong resonance interaction between C5–X and C4–NH₂ in the Ct C5–X 3H tautomer—there is a large increase in NH₂ lone pair delocalization upon X = CN substitution. Contrarily, in the Ct C5–X 1H tautomer, C5–X → C4–NH₂ conjugation is absent. In that case, the CN substituent instead interacts via resonance with the endocyclic N1(H) atom, increasing the delocalization along the X → N1(H) conjugation pathway, while leaving the C4–NH₂ lone pair unaffected.

In iCt C5–X 3H (Fig. 4b), ΔBDF_π reveals a large increase in delocalization along the para conjugated C5–X → C2–NH₂ upon X = CN substitution. In contrast, in iCt C5–X 1H this conjugation is disrupted, and the ΔBDF_π map shows no increase in delocalization between C5–CN and C2–NH₂ groups. This explains the changes in the relative SE strength between 1H and 3H (Fig. 4b) in terms of resonance interactions—the positions of double bonds and lone pairs forbid the charge transfer between the groups in the π-cloud, making the SE over two times weaker. Notably, ΔBDF_π allows us to identify when resonance interactions occur, visualize the conjugation pathways involved, and rationalize the observed differences in terms of resonance structures.

Fig. 4c shows a comparison of the SE in the G 9H,C6–OH and G 7H,C6–OH tautomers. In G 9H, ΔBDF_π reveals increased delocalization along two conjugated pathways, C8–X → C6–OH and C8–X → C2–NH₂, upon X = CN. This is accompanied by an increase in the SE strength, ∼7% for OH, and a larger, ∼36% for NH₂.

Interestingly, ΔBDF_π reveals that in 7H and 9H tautomers the changes in delocalization (a uniform increase) occur mainly in the five-membered ring, meaning that the CN group is highly engaged in resonance interactions with the five-membered ring π-system.

In Fig. 4d and e, ΔBDF_π shows that C8–NH₂ → C6=O conjugation in Hx/G 1H,7H tautomers slightly increases the delocalization along the pathway which passes through the central C4–C5 bond. In contrast, in 1H,9H tautomers, in which these two groups are nonconjugated, the same bond is localized upon NH₂ substitution.

In iG 1H,9H (Fig. 4e), ΔBDF_π indicates that X = CN substitution increases delocalization of the C6–NH₂ lone pair and two bonds along the conjugation pathway: N7–C8 and C5–N7. In G 1H,9H, ΔBDF_π for X = CN uncovers a long range C8–X → C2–NH₂ conjugation which increases delocalization throughout the conjugation pathway and C2–NH₂ lone pair. In conclusion, using ΔBDF_π we can see the resonance interactions appearing within the π-electron wavefunction, which nicely corroborates the analysis results of the changes in the SE strength.

Identifying resonance forms in heterocycles with BDFπ

In turn, BDF_π allows determining the dominant resonance forms in the π-electronic structure; apart from the examples presented in Fig. 4, and the BDF_π maps for all studied heterocycles are shown in Fig. S1–S8. BDF_π shows that the π-electronic structure is determined mostly by the tautomeric form.

In guanine (and its thio-, iso- forms and hypoxanthine), the tautomers with one N(H) in the five-membered ring and one or two in the six-membered ring are best described as having delocalized π-sextet in the five-membered ring, and localized bonds or lone pairs in the six-membered ring.

In tautomers without N(H) in the six-membered ring (C2/C6–OH/SH) the π-electronic structure is dominated by a delocalized π-sextet in the six-membered ring, but there is also some delocalization within the five-membered ring. EDDB_P(π) populations of cyclically delocalized electrons are between 1.9 and 3.0e for six-membered rings, and 0.9 and 1.6e for five-membered rings; for comparison, benzene has 5.4e and pyrimidine has 5.2e. This indicates that both structures, with the π-sextet in the six- and five-membered rings, contribute to the wavefunction, but the six-membered one dominates. Importantly, when electrons organize to form a cyclically delocalized π-sextet in one ring, then, from the resonance structures, the electrons in the second ring must localize at bonds and lone pairs.

Interestingly, the electron-withdrawing substituents in the C8 position increase the effectiveness of the five-membered ring delocalization in all cases. This is most likely associated with the X → N(H) resonance interaction (the lone pair at N(H) is donated onto π-withdrawing X, Fig. 2), which has a conjugation path around the five-membered ring circumference. The interaction increases the delocalization of the bonds and the N(H) lone pair in the five-membered ring, significantly improving the cyclic delocalization. In contrast, the π-electron-donating groups disrupt it; e.g. in G 3H,7H, EDDB_P(π) populations for the five-membered ring are 2.39e for X = CN, 1.89e for X = H and 1.30e for X = OH (Fig. S3).

This interaction can also be noticed in ΔBDF_π maps in Fig. 4c and e—in all maps for X = CN the delocalization in the five-membered ring uniformly increases due to substitution. However, strengthening delocalization in one ring weakens it in the other, due to the aforementioned competition between π-sextets in five- and six-membered rings. Consequently, EDDB_P(π) values for the two rings correlate linearly with C8-substituent properties (cSAR(X)), but with opposite slopes (Fig. S9)—an increase in the electron-withdrawing strength of the substituent increases the five-membered ring cyclic delocalization but decreases the six-membered one. So, the electron-withdrawing substituents in the five-membered ring of purine bases shift the relative contribution of resonance forms towards the one with delocalized π-sextet in the five-membered ring (Fig. S9). In 5-aza-7-deazaguanine, all tautomers are characterized by very weak cyclic delocalization in both rings and highly olefinic bonds (Fig. S7).

Cytosine (and isocytosine) tautomers with one or two N(H) in the ring have highly localized olefinic bonds—just one resonance form drawn in, e.g. Fig. 1b represents the π-electronic structure well. The tautomers without endocyclic N(H) groups, i.e. C2/C4–OH, have moderately delocalized π-sextet in the ring, and the EDDB_P(π) is between 1.9 and 2.9e (for comparison, in pyrimidine EDDB_P(π) = 5.2e).²³ This illustrates the disruptive effect of multi-substitution on the ring aromaticity.

Here, we would like to turn attention to the BDF_π maps for unsubstituted systems (X = H). (Iso)cytosines in Fig. 4a and b are nonaromatic—the π-electrons are mostly localized at bonds and lone pairs. However, for the guanine tautomers with N(H) in the five-membered ring and the OH group, i.e. the G 7H and 9H,C6–OH tautomers, the BDF_π map shows a very efficient delocalization within the six-membered ring—it is aromatic. This means that Kekulé resonance structures responsible for the aromaticity of this ring dominate the ground state electronic structure; these Kekulé resonance structures also force the electrons in the five-membered ring to localize. The low-energy aromatic resonance structures are resistant to perturbation by the high-energy zwitterionic structures associated with the SE—the resonance SE in the aromatic ring should be weaker.

Comparing the substituent effect in cytosine and guanine analogues

Fig. 5 and 6 show the substituent effect strength and the properties of CN, Cl and NH₂ groups in all studied heterocycles. In summary, the strongest SE of X on the NH₂ group occurs in C, iC and 6AP C5–X ortho or para conjugated systems (Fig. 5, Table S2). In guanine derivatives (Fig. 6) the SE is weaker, but two isoguanine derivatives with effective C8–X → C6–NH₂ conjugation between stand out: iG 1H,3H (a = −0.444) and iG 1H,9H (a = −0.311). This is because the C6–NH₂ group in iG is closer to C8–X than the C2–NH₂ group in G, tG, adG (Fig. 6).

Fig. 5 Properties (cSAR(X) values in boxes) of X = CN (red), Cl (green), and NH₂ (blue) substituents in all studied tautomers of C5 and C6-substituted cytosine (Ct), isocytosine (iCt) and 6-aminopyridin-2-one (6AP). Intensity of box shading is proportional to the strength of withdrawal or donation of electrons by X. Green and blue curved arrows represent the strength of the substituent effect of X substituent series on exocyclic groups (Y): =O/OH(=S/SH) and NH₂/NH, evaluated from the slopes of corresponding cSAR(Y) vs. cSAR(X) linear correlations, as illustrated in Fig. 1a.

Fig. 6 Properties (cSAR(X) values) of CN (red), Cl (green), and NH₂ (blue) substituents in all studied tautomers of C8-substituted hypoxanthine (Hx), guanine (G), thioguanine (tG), isoguanine (iG) and 5-aza-7-deazaguanine (adG). Darker cell shading indicates stronger withdrawal or donation of electrons. Green and blue curved arrows represent the strength of the substituent effect of X substituent series on exocyclic groups (Y): =O/OH(=S/SH) and NH₂/NH, evaluated from the slopes of the corresponding cSAR(Y) vs. cSAR(X) linear correlations, as illustrated in Fig. 1a.

In contrast, in iG 9H,C2–OH, there is a noticeably weaker SE (a = −0.246) than in iG 1H,3H and iG 1H,9H, despite all having C8–X → C6–NH₂ conjugation. This exemplifies how competition from six-membered ring aromatic resonance structures (present in iG 9H,C2–OH, not present in iG 1H,3H and iG 1H,9H) lowers the contribution of C8–X → C6–NH₂ charge transfer resonance structures. The resonance SE is better transmitted across non-aromatic rings (iG 1H,3H, iG 1H,9H) than the aromatic ones (iG 9H,C2–OH), which has been reported numerous times in the literature, e.g. by comparisons between the SE strength in benzene vs. cyclohexadiene²⁶ or polyene vs. acene,²⁷ among others.²⁸ This trend can also be noticed by comparing 6AP C5–X 1H vs. 6AP C5–X C2–OH and Ct C5–X 3H vs. Ct C5–X C2–OH (Fig. 5).

The substituent effect of X on the C=O group is the strongest in Ct, iCt and 6AP derivatives, moderate in Hx and G derivatives and very weak in adG and iG (Fig. 6, Table S3). The SE on the OH group is the strongest in 6AP C5–X and Ct C5–X (both have π-conjugated X and OH in the para relation). In Hx 9H,C6–OH and 7H,C6–OH, the substituent effect is stronger than in analogous G tautomers, which indicates that adding another functional group (C2–NH₂ in G) weakens the overall effect of X on C6–OH. The weakest SE is present in iG and adG derivatives, due to the distance (iG) or altered ring topology (adG). Generally, in adG derivatives, the SE is weak on both –NH₂/=NH and =O/–OH, owing to the altered ring topology which forbids π-conjugation between substituents—in these systems π-withdrawing groups are confined to a short-range interaction with the N5(H) group, while π-donating groups are not conjugated with any heteroatoms.

In general, =S groups are more sensitive to substituent effects than =O; comparing analogous tG and G systems reveals that SE on =S is stronger than on =O (compare the numbers near the green arrows for tG and G in Fig. 6). The reason behind it is that =S is a stronger electron-withdrawing group than =O, as evidenced by cSAR(=S) and cSAR(=O) (Table S4), as well as the cSAR(X) (for X = NH₂) values (Fig. 6) for tG and G systems. The latter indicates that more electrons are donated by X = NH₂ onto =S than =O. The electron-withdrawing strength of =S is associated with strong resonance interactions, as evidenced by BDF_π maps, which show much better delocalization of the C=S bond compared to C=O (Fig. S3 and S5).

On the other hand, OH groups are always more electron-donating than SH in their C6–OH/SH forms (Table S4), but SH groups are sensitive to substituent effects. This can be noticed when comparing SE strengths in tG C6–SH and G C6–OH (Fig. 6).

Substituent effect on tautomer stability

Finally, the substituent effect is able to shift tautomeric preferences (calculated Gibbs energies, ΔG, are collected in Table S1): (i) in Ct C5–X, the C2–OH form becomes slightly more stable than 1H for X = Br, Cl, and 1H, 3H, C4–NH for X = OH, NH₂; (ii) in Ct C6–X, C2–OH becomes more stable than 1H for all groups except X = H; (iii) in iCt C5–X, 3H becomes more stable than C4–OH for X = F, SH, OH, NH₂; (iv) in G C8–X, 1H,7H becomes slightly more stable than 1H,9H for X = F, SH, OH, NH₂; (v) in Hx C8–X, 1H,9H becomes slightly more stable than 1H,7H for X = NO₂, CN; (vi) in tG C8–X, 1H,7H becomes slightly more stable than 9H,C6–SH for X = NH₂; (vii) in iG C8–X, 1H,3H becomes slightly more stable than 9H,C2–OH for X = NH₂, F; and (viii) in adG C8–X, the 1H form becomes more stable than 9H for X = NO₂, CN, Br, Cl, F, OH, NH₂. In most cases, two tautomers have similar Gibbs energies, with differences <1 kcal mol⁻¹, which indicates that they may coexist in equilibrium. However, it should be remembered that solvation, which was not included in our study, may influence relative tautomeric stabilities.^29–31

Computational details

Geometries were optimized at the DSD-PBEP86-D3BJ (2013 parametrization)³² level with the def2-TZVPP basis set³³ in the Orca 6.0 program.^34,35 Following the optimizations, vibrational frequency calculations were performed to confirm the geometries correspond to the minima on the potential energy surface. Atomic charges were calculated using the Hirshfeld method.³⁶ The properties of exocyclic groups were evaluated using the cSAR method,^20,37 defined as a sum of atomic charges at substituent atoms and the ipso carbon atom to which the substituent is connected: cSAR(X) = q(X) + q(C_ipso). The cSAR values were corrected for the influence of endocyclic N atoms on the substituted carbon atom, the effect on electron density, which is not associated with the SE. The correction was performed in a way that from every cSAR(X) and cSAR(NH₂/NH) value, the cSAR(H) value for the reference system was subtracted (the same molecule, but the H atom replaces the X or NH₂/NH group). Substituent effect strength of X on NH₂/=NH and =O/OH/=S/SH groups was evaluated from linear correlations between the properties of X, cSAR(X), and the properties of the other group, e.g. cSAR(NH₂); the slope (a) of correlation (cSAR(NH₂) = a·cSAR(X) + b) represents the substituent effect strength of X on NH₂ in the given system, as illustrated in Fig. 1a. Natural bond orbital (NBO) calculations were performed using the NBO 7.0 program.³⁸ Electron density of delocalized bonds (EDDB) calculations were performed in the RunEDDB program,²⁴ calculations on .cube files with EDDB grid data to generate bond delocalization function (BDF_π)²⁵ maps were done using the Multiwfn program.³⁹ BDF_π is defined as BDF_π = EDDB_π–EDLB_π, where EDDB_π and EDLB_π are the densities of delocalized and localized electron densities (Fig. 1b) from EDDB density partitioning.^23,40–42

Conclusions

In conclusion, the strength of the substituent effect (SE) depends on: (i) the distance between the groups, which modulates both inductive and resonance contributions, (ii) the presence or absence of a π-conjugated path, required for resonance transmission, (iii) the effectiveness of resonance charge transfer along this path, which can be lower if other, highly-contributing resonance structures (e.g. aromatic ones) interfere with charge-transfer structures, and (iv) interference from other substituents. These principles, summarized in Fig. 7 using adenine as an example, allow qualitative prediction of whether the SE at a given site will be strong or weak in any N-heterocycle. The conjugation (or its lack) and its effectiveness depend on the positions of the double bonds and lone pairs, determined by the tautomeric form and heterocycle topology. The substituents interact with exocyclic groups, but also with endocyclic N atoms, which influences the amount of charge transferred onto (withdrawing) or from (donating) substituents. The π-electronic structure is critical for the SE strength and transmission when comparing heterocycles with the same distance between interacting groups. The bond delocalization function (BDF_π) allows explaining the differences in the SE strength between heterocycles by providing a deep insight into changes in (de)localization within π-electron cloud upon substitution (especially the differential BDF_π, ΔBDF_π). In many cases, ΔBDF_π allows visualizing the resonance interaction pathways, which makes it a powerful complement to traditional analysis of resonance structures. The extensive data on the substituent effect strength obtained in this study could be valuable for the training of Machine Learning models of the SE, while the four basic substituent effect rules should assist in rationalizing and predicting the SE strength in diverse heterocycles.

Fig. 7 Factors influencing the substituent effect transmission, illustrated using adenine (Ad) and guanine (G) tautomers. The SE strengths for Ad were calculated as slopes of cSAR(NH₂) vs. cSAR(X) correlations (X = NO₂, CN, CHO, Cl, F, H, Me, OMe, OH, and NH₂), data taken from ref. 32. The yellow numbers in the rings are the EDDB_P(π) populations of cyclically delocalized electrons, representing aromaticity of rings. Isosurfaces of BDF_π (+0.015, green, delocalization of electrons dominates, −0.015, brown, localization dominates) are superimposed onto chemical structures.

Author contributions

Conceptualization: PAW, HS, and TMK; formal analysis: PAW and HS; funding acquisition: PAW and HS; investigation: PAW and HS; methodology: PAW and HS; supervision: HS; validation: HS; visualization: PAW; writing – original draft: PAW and HS; and writing – review and editing: HS and TMK.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information: additional figures and tables, raw numerical data, and archive (.zip) containing optimized coordinates (.xyz) of all studied molecules. See DOI: https://doi.org/10.1039/d5cp03761a.

Acknowledgements

This research was funded by the Warsaw University of Technology within the Excellence Initiative: Research University (IDUB) programme (grant no. 1820/107/Z01/2023). The authors would like to thank the Wroclaw Center for Networking and Supercomputing for providing computer time and facilities.

References

1. O. Ebenezer, M. A. Jordaan, G. Carena, T. Bono, M. Shapi and J. A. Tuszynski, Int. J. Mol. Sci., 2022, 23, 8117.
2. D. Fan, D. Wang, T. Han and B. Z. Tang, ACS Appl. Polym. Mater., 2022, 4, 3120–3130.
3. A. F. Pozharskii, A. T. Soldatenkov and A. R. Katritzky, Heterocycles in Life and Society: An Introduction to Heterocyclic Chemistry, Biochemistry and Applications, Wiley, 1st edn, 2011.
4. B. Mohite, B. K. Chabhadiya, K. M. Kapadiya, V. M. Khedkar and S. Jauhari, ChemistrySelect, 2024, 9, e202404601.
5. P. Kumar, A. Soni, V. Tomar, T. Singh and M. Nemiwal, Sustainable Chem. Pharm., 2024, 41, 101714.
6. C. M. Marshall, J. G. Federice, C. N. Bell, P. B. Cox and J. T. Njardarson, J. Med. Chem., 2024, 67, 11622–11655.
7. E. Vitaku, D. T. Smith and J. T. Njardarson, J. Med. Chem., 2014, 57, 10257–10274.
8. S. Imran, Y. Gao, W. Wang and J. Zhu, Org. Biomol. Chem., 2025, 23, 4508–4518.
9. L. Simón-Vidal, O. García-Calvo, U. Oteo, S. Arrasate, E. Lete, N. Sotomayor and H. González-Díaz, J. Chem. Inf. Model., 2018, 58, 1384–1396.
10. G. S. Remya and C. H. Suresh, Phys. Chem. Chem. Phys., 2016, 18, 20615–20626.
11. S. Yokoyama and N. Nishiwaki, J. Org. Chem., 2019, 84, 1192–1200.
12. H. Szatylowicz, P. H. Marek, O. A. Stasyuk, T. M. Krygowski and M. Solà, RSC Adv., 2020, 10, 23350–23358.
13. C. Verma, A. Alfantazi, M. A. Quraishi and K. Y. Rhee, Coord. Chem. Rev., 2023, 495, 215385.
14. L. P. Hammett, J. Am. Chem. Soc., 1937, 59, 96–103.
15. L. P. Hammett, Physical Organic Chemistry, McGraw–Hill, New York, 1940.
16. C. Hansch, A. Leo and R. W. Taft, Chem. Rev., 1991, 91, 165–195.
17. O. Exner and S. Bohm, Curr. Org. Chem., 2006, 10, 763–778.
18. H. Szatylowicz, T. Siodla, O. A. Stasyuk and T. M. Krygowski, Phys. Chem. Chem. Phys., 2016, 18, 11711–11721.
19. O. A. Stasyuk, H. Szatylowicz, C. Fonseca Guerra and T. M. Krygowski, Struct. Chem., 2015, 26, 905–913.
20. N. Sadlej-Sosnowska, Pol. J. Chem., 2007, 81, 1123–1134.
21. N. Sadlej-Sosnowska, Chem. Phys. Lett., 2007, 447, 192–196.
22. D. W. Szczepanik, M. Andrzejak, K. Dyduch, E. Żak, M. Makowski, G. Mazur and J. Mrozek, Phys. Chem. Chem. Phys., 2014, 16, 20514–20523.
23. D. W. Szczepanik, M. Andrzejak, J. Dominikowska, B. Pawełek, T. M. Krygowski, H. Szatylowicz and M. Solà, Phys. Chem. Chem. Phys., 2017, 19, 28970–28981.
24. D. W. Szczepanik and M. Solà, Aromaticity, Elsevier, 2021, pp. 259–284.
25. D. Szczepanik and P. Wieczorkiewicz, ChemRxiv, 2025, preprint DOI:10.26434/chemrxiv-2025-m80qm-v4.
26. H. Szatylowicz, A. Jezuita, T. Siodla, K. S. Varaksin, K. Ejsmont, M. Shahamirian and T. M. Krygowski, Struct. Chem., 2018, 29, 1201–1212.
27. P. A. Wieczorkiewicz, M. Shahamirian, T. Kupka, N. Makieieva, T. M. Krygowski and H. Szatylowicz, Chem. Eur. J., 2024, 30, e202303207.
28. T. M. Krygowski, M. Palusiak, A. Płonka and J. E. Zachara-Horeglad, J. Phys. Org. Chem., 2007, 20, 297–306.
29. T. Bartl, Z. Zacharová, P. Sečkářová, E. Kolehmainen and R. Marek, Eur. J. Org. Chem., 2009, 1377–1383.
30. E. D. Raczyńska and B. Kamińska, J. Mol. Model., 2013, 19, 3947–3960.
31. A. Laxer, D. T. Major, H. E. Gottlieb and B. Fischer, J. Org. Chem., 2001, 66, 5463–5481.
32. S. Kozuch and J. M. L. Martin, J. Comput. Chem., 2013, 34, 2327–2344.
33. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297.
34. F. Neese, Wiley Interdiscip. Rev.:Comput. Mol. Sci., 2022, 12, e1606.
35. F. Neese, J. Comput. Chem., 2023, 44, 381–396.
36. F. L. Hirshfeld, Theor. Chim. Acta, 1977, 44, 129–138.
37. N. Sadlej-Sosnowska, Chem. Phys. Lett., 2007, 447, 192–196.
38. E. D. Glendening, J. K. Badenhoop, A. E. Reed, J. E. Carpenter, J. A. Bohmann, C. M. Morales, P. Karafiloglou, C. R. Landis and F. Weinhold, NBO 7.0 Theoretical Chemistry Institute, University of Wisconsin, Madison, 2018.
39. T. Lu and F. Chen, J. Comput. Chem., 2012, 33, 580–592.
40. D. W. Szczepanik, Comput. Theor. Chem., 2016, 1080, 33–37.
41. D. W. Szczepanik, M. Andrzejak, K. Dyduch, E. Żak, M. Makowski, G. Mazur and J. Mrozek, Phys. Chem. Chem. Phys., 2014, 16, 20514–20523.
42. D. W. Szczepanik, E. Żak, K. Dyduch and J. Mrozek, Chem. Phys. Lett., 2014, 593, 154–159.

---