A molecular electron density theory investigation of the mechanism of intramolecular [3+2] cycloaddition (32CA) with the participation of nitrile N-oxide and ethene molecular segments

Abstract

In this study, we have investigated the reaction mechanism of the intramolecular [3+2] cycloaddition (32CA) process involving a nitrile oxide-heterocycle system through molecular electron density theory (MEDT). Density functional theory (DFT) calculations at the B3LYP/6-311++G(d,p) level of approximation have been performed to explore both endo and exo stereoisomeric pathways in benzene as a solvent. The results clearly indicate that the exo path is both thermodynamically and kinetically preferred, consistent with experimental findings. This preference is shown by reduced activation energies and greater negative Gibbs free energies for the exo product relative to the endo product. The bonding evolution theory (BET) analysis, in combination with ELF topological analysis, unraveled a stepwise mechanism that involves the formation of nitrogen lone pairs followed by the formation of C–C and C–O bonds. This mechanistic interpretation points out the asynchronicity of the bond-forming process, and the exo pathway is found to be less asynchronous compared to the endo one. Moreover, molecular docking analyses revealed that the exo product has considerable binding affinity to the 1CIN protease, indicating its potential as a therapeutic inhibitor. Moreover, drug-likeness evaluations verified that the compounds adhere to Lipinski's Rule of Five, signifying advantageous pharmacokinetic characteristics. This extensive work combines theoretical and computational approaches to clarify the intricate processes of 32CA reactions, offering significant insights into their synthetic applications and possible medicinal benefits.

---

1. Introduction

Intramolecular [3+2] cycloaddition (32CA) reactions of nitrile oxides with heterocycles play a critical role in organic synthesis,¹ providing a direct and efficient route for assembling intricate molecular architectures.^2–4 These reactions exhibit exceptional regio- and stereoselectivity, making them invaluable for constructing pharmacologically significant compounds and natural product derivatives.^5–7 However, despite their synthetic importance, the mechanistic complexities of these reactions, particularly the interplay between concerted and stepwise pathways, remain an area of active research.^8–10

One of the most pursued aspects in the reaction mechanism of 32CA has been the question of concertedness and asynchronicity in the formation of bonds, and the factors that control the electronics of the reactants and solvent effects, among other factors.^11,12 Related comparative studies of cycloaddition systems have underlined how advanced computational techniques, such as density functional theory, will be useful in accounting for the intrinsic reaction pathways and energy profiles. For instance, Domingo et al. have illustrated that molecular electron density theory (MEDT)¹³ is a powerful tool to understand the polar character and the electron density flow in most of the cycloaddition processes.^14,15

The bonding evolution theory (BET),¹⁶ based on electron localization function (ELF) analysis,¹⁷ is a powerful theoretical framework that describes bond forming and cleavage events in chemical reactions.^18–20 It allows identification of the elementary chemical events due to dynamical evolution of bonding interactions along the intrinsic reaction coordinate (IRC).^21–23 This analysis has turned out to be quite effective in distinguishing between synchronous and asynchronous mechanisms and for understanding subtle reactivity differences in stereoisomeric pathways.^24,25

In the context of 32CA reactions, the mechanistic pathways are categorized into four distinct types according to the MEDT classification, pseudodiradical (pdr), pseudo(mono)radical (pmr), carbenoid (cb), and zwitterionic (zw). These types are defined based on the nature of the three-atom components (TACs) involved. Specifically, a monosynaptic basin incorporating one electron is classified as a pseudoradical center, while a basin incorporating two electrons is classified as a carbenoid center. TACs containing two pseudoradical centers are termed pseudodiradicals, those containing one pseudoradical center are called pseudo(mono)radicals, and TACs lacking carbenoid or pseudoradical centers are considered zwitterionic. The reactivity order among these types is established as pdr > pmr > cb > zw, with pdr-type reactions advancing more rapidly through early transition states compared to zw-type reactions, which require additional nucleophilic or electrophilic activation to achieve satisfactory reaction rates. This classification system is essential for understanding the diverse mechanistic pathways in 32CA reactions and provides valuable insights into the factors influencing reactivity and selectivity.²⁶

The comparison of the endo and exo pathways in intramolecular 32CA reactions establishes striking differences in their thermodynamic and kinetic inclinations. For example, it has been evidenced that exo pathways are generally preferred when considering lower activation barriers and higher thermodynamic stabilities, a trend which has also been experimentally and theoretically proved in ref. 27 and 28. These analyses are further refined upon inclusion of solvent effects, since solvent polarity can drastically alter the reaction profile by offering stability either to the transition states or to the intermediates.^29,30

The 32CA reaction, involving the interaction of a nitrile N-oxide and ethene molecular fragments (Scheme 1), was performed in benzene at ambient temperature, resulting in a 77% yield of the final bicyclic isoxazolidine product, 5,6-dihydro-[1,2,4]oxadiazolo[3,4-b]indole.³¹ This reaction follows two distinct pathways, leading to the formation of endo 2 and exo 3 cycloadducts (Scheme 2). Experimental data reveal that the exo cycloadduct 3 is the predominant product, highlighting a clear preference for exo selectivity. The experimental conditions were selected based on previous studies suggesting that benzene, a non-polar solvent, promotes the formation of the exo product due to minimized steric hindrance during the cycloaddition step. Temperature control was also crucial, as reactions performed at higher temperatures led to the formation of undesired side products. The origin of this selectivity is the primary focus of this study. These experimental observations are further explored and explained through computational studies presented herein. The motivation for this research arises from the need to understand the mechanisms driving regio- and stereoselectivity in [3+2] cycloaddition reactions involving nitrile N-oxides, which are key intermediates in the synthesis of bioactive heterocyclic compounds.

Scheme 1 Experimental conditions of the intramolecular [3+2] nitrile oxide-heterocycle.

Scheme 2 Endo and exo reaction pathways of the intramolecular 32CA with the participation of nitrile N-oxide and ethene molecular segments.

In this work, density functional theory (DFT) is used to carry out a thorough mechanistic study of the 32CA reaction between nitrile oxide-heterocycle 1 (Scheme 2). The thermodynamic and kinetic parameters of the reaction are analyzed in solution (benzene), providing an in-depth examination of the structural, electronic, and thermochemical properties of the reactants, transition states, and products. A topological analysis based on the ELF within the BET framework is employed to describe the sequence of O–C and C–C bond formations and structural changes along the intrinsic reaction coordinate (IRC). By elucidating these factors, our study offers valuable insights into the fundamental aspects of heterocyclic cycloadditions. Furthermore, molecular docking studies are conducted to assess the potential therapeutic applications of the cycloadducts as protease inhibitors, complementing the mechanistic findings.^32–35 This work bridges the gap between theory and experiment, offering a comparative analysis of the endo and exo pathways, and providing a theoretical framework for optimizing cycloaddition reactions.

2. Theoretical chemistry and computational aspects

The equilibrium structures of the reactants, products, and transition states were fully optimized at the density functional theory (DFT) level using the B3LYP exchange–correlation functional³⁶ and the 6-311++G(d,p) basis set. The B3LYP functional has been widely recognized as an appropriate computational method for studying 32CA reactions in numerous recent studies.^37–39 Transition states (TSs) were identified by the presence of a single imaginary frequency, while minima were confirmed by the absence of any imaginary frequencies. The global electron density transfer (GEDT)⁴⁰ was calculated from the natural bonding orbitals (NBO)⁴¹ at the TSs using the formula:where q denotes the natural atomic charges derived from NBO analysis. To verify the reaction pathway and ensure the minimum energy connection between reagents, transition states, and products, IRC calculations were performed.^42,43 Solvent effects were examined for benzene (ε = 6.0) using the IEFPCM method (integral equation formalism of the polarizable continuum model).⁴⁴ All calculations were carried out at 298.15 K and an initial pressure of 1.0 atm.

Following the Bader idea, Becke and Edgecombe have proposed the ELF, which measures the maximal probability for finding an electron pair in a molecular species (molecule or crystal).¹⁷ The electron localization function (h) is defined by the following equation:

The D_σ(r) and D⁰_σ(r) quantities are measures of the electronic localization for ordinary and homogeneous gases, respectively.where |∇_ρσ|² is defined as the kinetic energy and ρ_σ(r) is the electron density for a given σ spin.

The ELF analysis was carried out along the IRC and led to BET pictures, which describe a chemical reaction as a sequence of elementary chemical processes separated by catastrophes. In this definition of a chemical reaction, the catastrophe corresponds to bond formation or bond breaking. A topological analysis of ELF along the reaction pathway gives attractor basins (domains), in which the probability is maximal for finding an electron pair. The basins can be classified into two types: core and valence basins. The valence basins can be monosynaptic, disynaptic, trisynaptic, and so on, depending on the number of atomic valence shells.⁴⁵ The ELF topological analysis of a structure allows characterization of three types of valence basins: protonated basins or V(A,H), monosynaptic basins V(A) corresponding to a nonbonding region, and V(A,B) corresponding to a bonding region.⁴⁶

For the topological analysis within the framework of the BET theory, the wave function was extracted at each point along the IRC, and the ELF calculation was performed using the TopMod package,⁴⁷ with a cubic grid resolution finer than 0.2 Bohr. This calculation yields the ELF basins, which may be classified as bonding, hydrogenated, or non-bonding, along with the electron population of each basin (i.e., each bond or lone pair). The evolution of these basin populations along the IRC was analyzed, and arrows representing the flow of electron density were drawn accordingly. The changes in basin population along the IRC and the ELF isosurfaces of selected basins were visualized using Drawprofile and Drawmol.^48,49 All calculations were carried out with the Gaussian16 software package.⁵⁰

3. Results and discussion

3.1. Analysis of the global and local CDFT indices

The analysis of conceptual DFT (CDFT) reactivity indices provides an initial understanding of the direction of electronic flow between reactants.^51–54 The global reactivity indices are defined by the following equations:

η = ε_LUMO − ε_HOMO (5)

N = ε_HOMO (nitrile oxide-heterocycle) − ε_HOMO (tetracyano-ethylene (TCE)) (8)

In this context, global reactivity indices show that the electrophilicity (ω) and nucleophilicity (N) indices of nitrile oxide-heterocycle are 1.1 eV and 3.5 eV (Table 1), respectively, being therefore classified as a marginal electrophile and a strong nucleophile within their respective scales.⁵⁵

Table 1 Electronic parameters and global reactivity indices [eV]

ε HOMO	ε LUMO	μ	η	ω	N
−5.8	−0.7	−3.3	5.1	1.1	3.5

---

In a polar 32CA reaction involving non-symmetric species, the preferred reaction pathway typically involves a two-center interaction between the most electrophilic center of one reagent and the most nucleophilic center of the other.⁵⁶ Numerous studies have highlighted the importance of electrophilic (P⁺_k) and nucleophilic (P⁻_k) Parr functions.⁵⁷ These functions, derived from the excess spin electron density accumulated through the GEDT, are among the most reliable tools for analyzing local reactivity in polar and ionic processes.^58,59

From the electron densities of the nitrile oxide-heterocycle 1, as well as its radical anion and cation (Fig. 1), local reactivity descriptors were calculated (Table 2). The results indicate that the O3 oxygen atom of the nitrile N-oxide segment exhibits nucleophilic activation (P⁻_k = −0.17), whereas the N2 nitrogen atom is nucleophilically deactivated (P⁻_k = −0.54). Additionally, the C4 carbon atom of the ethene moiety also displays nucleophilic characteristics (P⁻_k = −0.72), making it the most nucleophilic center of the nitrile oxide-heterocycle 1. Furthermore, the C1 carbon atom of the nitrile N-oxide moiety is the most electrophilically activated site (P⁺_k = 0.98), while the C5 position of the ethene moiety also shows electrophilic activation (P⁺_k = 0.50).

Fig. 1 3D representations of the spin densities (η = 0.02) of the radical anion and radical cation of nitrile oxide-heterocycle 1 together with its nucleophilic (P⁻_k) and electrophilic (P⁺_k) Parr functions.

Table 2 Local electrophile (ω_k) and nucleophile (N_k) indices, local reactivity difference index values (R_k), and nucleophilic (P⁺_k) and electrophilic (P⁻_k) Parr functions of nitrile oxide-heterocycle 1

	Center	P ^+ k	P ^− k	ω k	N k	R k
1	O3	0.11	−0.17	0.20	−0.46	0.66
	N2	−0.47	−0.54	−0.87	−1.49	0.62
	C1	0.98	0.40	1.83	1.09	0.74
	C4	−0.72	−0.72	−1.21	−2.01	0.80
	C5	0.50	−0.11	0.84	−0.31	1.15

---

Notably, C4 is more nucleophilic than O3, and C1 is more electrophilic than C5. This suggests a preferred two-center interaction between C1 and C4 during the intramolecular 32CA reaction.

Note that C4 is a more nucleophilic center than O3, while C1 is a more electrophilic center than C5, thus predicting a feasible two-center interaction between C1 and C4 along the intramolecular 32CA reaction. Consequently, C4 initiates a nucleophilic attack on C1, leading to the formation of the C1–C4 bond, followed by the formation of the O3–C5 bond.

3.2. ELF topological analysis of nitrile oxide-heterocycle 1

Topological analysis of the ELF establishes a quantitative correlation between the electronic structure and the reactivity of TACs participating in 32CA reactions.⁶⁰ The electronic structure of nitrile oxide-heterocycle 1 was then characterized by calculating the natural population analysis (NPA)^61,62 charges and by performing a topological analysis of their most representative ELFs, together with their basin populations, leading to the ELF-based Lewis structures given in Fig. 2. At the nitrile N-oxide moiety, the ELF topology shows the presence of one monosynaptic basin, V(O3), integrating a total population of 5.67e, and two V(N2,C1), and V(O3,N2) disynaptic basins integrating a total population of 6.19e, and 1.54e, respectively. The V(O3) monosynaptic basin can be associated with the non-bonding electron density on the O3 oxygen atom. V(N2,C1) and V(O3,N2) disynaptic basins can be associated, respectively, with the underpopulated N2–C1 triple bond and O3–N2 single bond. In addition, the ELF of 1 shows the presence of another disynaptic basin, V(C4,C5,) integrating a total of 3.32e, associated with the underpopulated C4–C5 double bond at the ethene moiety. Within the nitrile N-oxide framework, NPA charges indicate that the O3 oxygen atom is negatively charged by −0.42e, and the N2 nitrogen and C1 carbon atoms are positively charged by 0.17e and 0.21e, respectively. In addition, at the ethene segment, the C4 carbon exhibits a negligible positive charge of +0.01e, whereas the C5 carbon is negatively charged at −0.27e.

Fig. 2 ELF basin isosurfaces (η = 0.80) and their populations (in e) of key basins of C nitrile oxide-heterocycle 1 and its Lewis-like structures, together with the NPA charges (in e, positive charges in blue, negative in red). Protonated basins are shown in blue, monosynaptic basins are shown in red, disynaptic basins are shown in green, and core basins are shown in magenta color.

3.3. Thermodynamical and geometrical aspects

Table 3 displays the reaction and activation energies for the stationary points involved in the studied intramolecular 32CA reactions (see Scheme 2).

Table 3 Relative energies (ΔE, in kcal mol⁻¹), enthalpies (ΔH, in kcal mol⁻¹), entropies (ΔS, in cal mol⁻¹ K⁻¹), and Gibbs free energies (ΔG, in kcal mol⁻¹) for all stationary points involved in the intramolecular reaction with respect to the reactants as a function of the method of calculation at 25 °C^a

	ΔE		ΔH		ΔS		ΔG
a In parentheses are given endo-to-exo differences of each thermodynamic parameter. [ΔE = X(endo) − X(exo), X = E, H, S, and G].
2	−11	(3.2)	−10.1	(13.1)	−16.7	(1.7)	−5.2	(2.5)
3	−14.2		−13.2		−18.4		−7.7
TSN	23.5	(2.6)	22.3	(2.7)	−11.3	(1.4)	25.7	(2.1)
TSX	20.9		19.6		−12.7		23.4

---

DFT calculations indicate that this 32CA reaction is both exothermic (ΔH(2/3) < 0) and exergonic (ΔG(2/3) < 0), with negative reaction free energies (ΔE), suggesting kinetic control and irreversibility. Moreover, the formation of the exo cycloadduct (3) is slightly more exothermic (ΔH(2) = −10.1 kcal mol⁻¹ > ΔH(3) = −13.2 kcal mol⁻¹, Fig. 3) and more exergonic (ΔG(2) = −5.2 kcal mol⁻¹ > ΔG(3) = −7.7 kcal mol⁻¹, Fig. 4) than the endo product (2). These results indicate that the exo product is thermodynamically slightly more stable than the endo counterpart [ΔX(3) < ΔX(2), with X = E, H°, and G°]. Furthermore, the exo cycloadduct is kinetically preferred over the endo one [ΔX(TSX) < ΔX(TSN), with X = E, H°, and G°]. The above statement is in full agreement with the experimental results, which demonstrate the formation of exo cycloadduct 3.³¹

Fig. 3 Enthalpy profile for the 32CA reaction of the nitrile oxide-heterocycle 1.

Fig. 4 Gibbs free energy profile for the 32CA reaction of the nitrile oxide-heterocycle 1.

The optimized geometries of the TSs and cycloadducts are depicted in Fig. 5 together with the length of new forming C1–C4 and O3–C5 bonds (see Scheme 2 for the nomenclature). For TSN, the C1–C4 and O3–C5 bond lengths forming are 1.964 Å and 2.769 Å, respectively. In TSX, these distances are 1.990 Å and 2.721 Å, respectively. This indicates that the C1–C4 distance is shorter than the distance between O3 and C5 atoms.

Fig. 5 Optimized geometries of the TSs and cycloadducts (bond lengths in Å) involved in these intramolecular 32CA reactions.

To assess the asynchronicity of the studied 32CA reaction, an analysis was conducted on all stationary points. Key parameters, such as the interatomic distance between reaction centers (r), the distance progress index (l) and the asymmetry index Δl, as well as the dipole moments are presented in Table 4.^40,63

Table 4 Key parameters of the critical structure parameters of the studied 32CA reaction

	r C1–C4 (Å)	l C1–C4	r O3–C5 (Å)	l O3–C5	Δl	GEDT (e)	Dipole moment (D)
1	3.395		4.299			0.29	5.91
TSN	1.964	0.75	2.769	0.12	0.63	0.39	5.45
2	1.568		1.473			0.28	3.89
TSX	1.99	0.69	2.721	0.14	0.55	0.39	5.27
3	1.517		1.465			0.28	3.97

---

The distance progress indexes of the forming l_C1–C4 and l_O3–C5 are 0.75 and 0.12 for TSN, respectively. For TSX, the distance progress indexes of the forming l_C1–C4 and l_O3–C5 are 0.69 and 0.14, respectively. These values indicate that for both the TSN and TSX, the formation of the C1–C4 bond is more advanced than the O3–C5 bond. This information shows that both TSN and TSX present a significant asynchronicity, with the TSN slightly more asynchronous than the TSX. In addition, it can be concluded that the more favorable TSX is slightly earlier than TSN. On the other hand, the dipole moment of the endo TS (5.45 D) is greater than that of the exo one (5.27 D), which explains the higher activation Gibbs free energy of TSN compared to TSX (Table 3).

To assess the polar nature of this intramolecular reaction, the GEDT at the TSs was evaluated (Table 4). The GEDT values are 0.39e for both the TSN and TSX pathways, indicating a polar character for this 32CA reaction, following Domingo's classification⁴⁰ where reactions with GEDT values of 0.0e are considered non-polar, while those with GEDT values greater than 0.2e are classified as polar processes. The positive GEDT value shows the electron density flowing nitrile N-oxide segment to the ethene one, supporting a forward electron density flux (FEDF) reaction.⁶⁴

3.4. Topological analysis of the ELF at the TSs

The ELF of the TSs was studied to comprehend their electronic structures. Fig. 6 illustrates the basin attractor locations and ELF localisation domains for the transition states in the intramolecular 32CA reactions. This analysis highlights the O3–C5 and C1–C4 bonding interactions.

Fig. 6 ELF localization domains and the positions of the TS attractor basins (η = 0.80). Monosynaptic basins are shown in red, protonated basins in blue, disynaptic basins in green, and attractor positions in purple.

The ELF analysis indicates that both TSN and TSX are characterized by the presence of a monosynaptic basin V(O3) associated with the non-bonding electron density at the O3 oxygen, integrating a total population of 5.68e and 5.67e, respectively. Additionally, the ELFs show the presence of the V(N2,C1) disynaptic basin associated with the N2–C1 bonding region integrating total population of 3.41e and 3.40e in TSN and TSX, respectively. Furthermore, the two TSs show the presence of a V(N2) monosynaptic basin, integrating 2.20e (TSN) and 2.19e (TSX), associated with the non-bonding electron density created at the N2 nitrogen. The V(O3,N2) disynaptic basin integrating a total population of 1.45e in both TSs, which can be associated with the underpopulated O3–N2 single bond. ELF of TSs shows the presence of a first monosynaptic basin, V(C1), integrating 0.44 and 0.45e, for TSN and TSX, respectively. At the same time, we can observe the presence of a second monosynaptic basin, V(C4), integrating 0.47, and 0.49, for the TSN and TSX, respectively. These basins correspond to the two pseudoradical centers at the C1 and C4 carbon atoms needed for the subsequent formation of the C1–C4 single bonds.

3.5. AIM and NCI topological analysis of the TSs

Bader and colleagues developed the topological analysis of the atoms in molecules (AIM) method,^65,66 which is employed to investigate interactions between atoms during transition states. The Laplacian of the electron density Δ²ρ(r), and the total electron density (ρ) of the bond critical points (BCPs) are presented in Table 5. These values are associated with the formation of the O3–C5 bond. Specifically, Δ²ρ(r) for CP1(O3-C5) in TSN is 0.0496 and in TSX is 0.0464. Additionally, empirical data suggests that an interatomic distance greater than 2.0 Å is essential for bond formation. Further insights can be gained from the reduced density gradient (RDG) plot in Fig. 7 and the non-covalent interaction (NCI) gradient isosurfaces depicted in Fig. 8.

Table 5 The nitrile oxide-heterocycle undergoes a 32CA reaction, with the total electron density (ρ, in a.u.) and the Laplacian of the electron density Δ²ρ(rc) at the critical points in the transition states corresponding to this reaction

Structure	CP(O3–C5)
	ρ	∇^2ρ
TSN	0.01896	0.0496
TSX	0.01762	0.0464

---

Fig. 7 NCI gradient isosurfaces (η = 0.08) of the optimized TSs (TSN and TSX) associated with the 32CA reaction of nitrile oxide-heterocycle 1.

Fig. 8 Reduced density gradient for the optimized TSs (TSN and TSX) associated with the 32CA reaction of nitrile oxide-heterocycle 1.

3.6. BET analysis along the endo and exo-approaches

To better understand the mechanism of the 32CA intramolecular reaction of the nitrile oxide-heterocycle, BET analyses were performed for all the studied pathways. For each BET analysis, the reported basins are those with an evolving electron population along the IRC.

3.6.1. BET analysis of the endo path (TSN). The BET study of the TSN stereoisomeric pathway associated with the 32CA intramolecular reaction indicates that the process takes place along seven structural stability domains (SSDs), as we can see in Fig. 9, which represents the populations of the basins involved in the cycloaddition. In the SSD-I [d(C1–C4) = 3.01 Å and d(O3–C5) = 3.43 Å], the system shows the topologies of the reactant as reported in Table S1 (ESI†). Four basins are observed: the disynaptic basin V(C4,C5) found between C4 and C5 atoms with a total population of 3.34e, the disynaptic basin V(N2,C1) found between N2 and C1 atoms with a total population of 6.13e, the disynaptic basin V(O3,N2) with a total population of 1.67e and the V(O3) basin associated with the three O3 oxygen lone pairs integrating a total of 5.72e (Table S1, ESI†). At SSD-II [d(C1–C4) = 2.27 Å and d(O3–C5) = 2.84 Å], the population of the V(N2,C3) basin drops (to attain 4.73e) due to the formation of a new monosynaptic basin V(N2) (with a population of 1.68e) via a fold-F type catastrophe. In the following domain (SSD-III) [d(C1–C4) = 2.03 Å and d(O3–C5) = 2.74 Å], we observe the creation of two monosynaptic basins, V(C1) and V(C4) with 0.36e and 0.39e, that arise from the simultaneous drop of V(N2,C1) and V(C4,C5) disynaptic basin populations, respectively, by means of two fold-F type catastrophes (Fig. 10). These two basins join together in domain IV [d(C1–C4) = 1.77 Å and d(O3–C5) = 2.65 Å], to form the new disynaptic V(C1,C4) basin, which corresponds to the C1–C4 σ-bond by the means of a cusp-C type catastrophe. Then, the V(C1,C4) basin population increases while the V(C4,C5) decreases by 0.20e (Table S1, ESI†). In the SSD-V [d(C1–C4) = 1.55 Å and d(O3–C5) = 2.06 Å], the population of the V(C4,C5) basin drops due to the formation of the new monosynaptic V(C5) basin (of a population of 0.14e) by means of a fold-F catastrophe. Then, at the end of this domain, the populations of V(C1,C4) and V(N2) increase to attain 1.96e and 2.74e, respectively (Fig. 9).

Fig. 9 Population (in electrons) evolution of selected basins along the IRC associated with the TSN stereoisomeric channel together with the potential energy surface along the reaction coordinate.

Fig. 10 ELF basin isosurfaces (η = 0.75) of each of the SSDs found along the IRC associated with TSN. Color labeling of the basins is adopted according to Fig. 9, and, for each point, the corresponding intrinsic reaction coordinate is provided.

In the SSD-VI [d(C1–C4) = 1.55 Å and d(O3–C5) = 1.90 Å], a second monosynaptic basin V′(O3) (with a population of 0.27e) appears on the oxygen atom (Fig. 9) due to the drop of the V(O3) basin population. The V(O3,C5) basin is formed in the VII domain [d(C1–C4) = 1.53 Å and d(O3–C5) = 1.75 Å], by the merger of V′(O3) (which loses 0.48e) and V(C5) (which loses 0.27e) basins (Fig. 9) by means of another cusp-C type catastrophe. Then, the population of the V(O3–C5) basin rises from 0.85e at the beginning of the VII domain to reach 1.17e at the end of the domain (Fig. 9).

3.6.2. BET analysis of the exo path (TSX). Like the TSN stereoisomeric channel, we have also conducted BET study along the TSX pathway, and seven SSDs have also been found for this analysis (see Fig. 11), and the same is observed for the endo pathway, with only slight differences in the population of the basins (Fig. 11 and 12). At SSD-I [d(C1–C4) = 3.05 Å and d(O3–C5) = 3.43 Å] the system shows the topologies of the reactant with three basins: the V(C4,C5) basin (3.34e) along the C–C double bond, the V(O3) basin on the oxygen atom lone pairs (5.75e), a V(N2,C1) basin with 6.11e for the triple C–N bond, and the V(O3,N2) basin with 1.65e for the single O–N bond (Table S2, ESI†).

Fig. 11 Population (in electrons) evolution of selected basins along the IRC associated with the TSX stereoisomeric channel together with the potential energy surface along the reaction coordinate.

Fig. 12 ELF basin isosurfaces (η = 0.75) of each of the SSDs found along the IRC associated with TSX. Color labeling of the basins is adopted according to Fig. 9, and, for each point, the corresponding intrinsic reaction coordinate is provided.

In SSD-II [d(C1–C4) = 2.34 Å and d(O3–C5) = 2.83 Å] as in TSN, a pseudoradical center on the N2 atom (1.59e) appears first coming from the drop of the V(C1,N2) (from 4.73 to 4.49e) populations (Fig. 10). The beginning of domain SSD-III [d(C1–C4) = 2.03 Å and d(O3–C5) = 2.74 Å] deals with the creation of a two new other monosynaptic basins, V(C1) and V(C4) with 0.33e and 0.37e, corresponding to a fold-type F catastrophe. Its population comes from the reduction in the population of disynaptic V(N2,C1) and V(C4,C5). Furthermore, the electron population of the monosynaptic basins V(C1), V(C4) and V(N2) keeps growing steadily and hits the pot at 0.60, 0.65, and 2.33e, respectively, at the end of the domain (Table S2, ESI†). The passage from SSD-III to SSD-IV [d(C1–C4) = 1.81 Å and d(O3–C5) = 2.61 Å] describes the formation of the C1–C4 single bond with the transformation of the C4–C5 double bond into a single one and the N2–C1 triple bond into a double one, which corresponds to the cusp-type catastrophe, Fig. 12. The population of 1.36e, for the new V(C1,C4) disynaptic basin, comes from the disappearance of the V(C1) and V(C4) monosynaptic basins. In the following domain (SSD-V) [d(C1–C4) = 1.56 Å and d(O3–C5) = 2.08 Å], another C5 pseudoradical center integrating a population of 0.16e, is created (Table S2, ESI†). The electron density of this pseudoradical center is a consequence of the continuous depopulation of the C4–C5 bonding region by about 0.24e. At SSD-VI [d(C1–C4) = 1.54 Å and d(O3–C5) = 1.86 Å], a second monosynaptic V′(O3) (with 0.37e), the oxygen atom O3 is created simultaneously with the drop of the population of the V(O3) basin. Finally, at the beginning of the last domain, SSD-VII [d(C1–C4) = 1.53 Å and d(O3–C5) = 1.71 Å], the two new pseudoradical center V′(O3) and V(C5) basins disappear by combining together in order to form the disynaptic V(O3,C5) basin, which illustrates the formation of an O3–C5 bond (Fig. 11). The new V(O3–C5) basin holds a population of 0.83e and reaches 1.22e at the end of the domain (Table S2, ESI†).

Drawing on the BET analysis for the two pathways, the electron flows leading from TS to the cycloadduct can be divided into four stages, as depicted in Scheme 3: (i) depopulation of V(C1,N2), giving rise to monosynaptic basins on N2 and on C3; (ii) depopulation of V(C4,C5), with the creation of a monosynaptic basin on C5; (iii) formation of the C1–C4 single bond from the C1 and C4 monosynaptic basins and creation of a C4 monosynaptic basin by depopulating the C4–C5 double bond; and (iv) formation of the C5–O3 single bond from the C5 and O3 monosynaptic basins.

Scheme 3 Flow of the electrons along the two reaction paths TSN and TSX. (The half arrows symbolize the transfer of a portion of electrons).

In addition to the global topology analysis, we compare the reaction coordinates of the appearance of the basins involved in the formation of the C1–C4 and O3–C5 bonds (Table S3, ESI†). For both TSN and TSX channels, the reaction coordinate of the appearance of the second topological change (V(C1) and V(C4)) is almost the same (0.32 Bohr ). For the V(C1,C4) one, they appear at 0.96, and 0.95 Bohr , for the endo and exo pathways, respectively. For the V(C5) basin, the reaction coordinates of appearance are all different and range from 4.79 to 4.43 Bohr . Therefore, the difference of reaction coordinate between the appearances of the V(C1) and V(C5) basins decreases when going from TSN to TSX, in agreement with the asynchronicity criterion (Δd) discussed in the Thermodynamical and Geometrical Aspects section.

3.7. Drug-likeness assessment and ADMET predictions

To assess the drug-likeness of the compounds, their physicochemical properties are often evaluated using specific filter modifications. Table 6 highlights the physicochemical parameters generated via the ADMETlab 3.0 Web server. Similarly, the radar charts in Fig. 13 show the physicochemical properties of the two compounds positioned between the lowest (red) and highest (brown) bounds. The two compounds that complied with Lipinski's Rule of Five underwent further analysis using the SwissADME Web server (Table 6). This evaluation included additional drug-likeness filters, such as Lipinski's filter rule.

Table 6 SwissADME-computed drug-likeness predictions for compounds 2 and 3

Ligands	MW	NROT	NHA	NHD	TPSA	Log P	Lipinski's violations
2	202.11	0	3	1	24.5	1.579	Accepted
3	202.11	0	1	1	24.5	1.607	Accepted

---

Fig. 13 Physico-chemical radar chart of the selected lead compounds in the dataset.

For the theoretical prediction of drug-like properties of compounds 2 and 3, SwissADME was employed, incorporating Lipinski's rule to identify the best-hit compounds. It is important to note that small molecules failing to meet multiple conditions of these rules may face bioavailability challenges. The pharmacokinetic analysis of the novel compounds is detailed in Table 6. Notably, each compound meets all five criteria of Lipinski's Rule of Five,^67,68 including the molecular weight requirement, which limits the molecular weight to a maximum of 500 g mol⁻¹. Furthermore, each of the best-hit compounds contained one hydrogen bond acceptor, which falls within the acceptable range. The number of hydrogen bond donors was also within the permissible limit, remaining below 10. Additionally, the calculated Log P values stayed within the acceptable ranges. Based on these filtering criteria, the best-hit compounds were classified as drug-like, as they did not exceed the established threshold values or violate more than one of the filtering conditions. Similarly, both compounds, 2 and 3, meet all the filtering requirements, indicating that they do not present any bioavailability issues. Fig. 11 demonstrates that both compounds share characteristics with the boiled egg model, reflecting properties typical of a drug.

However, the BBB (blood–brain barrier) permeability values for compounds 2 and 3 are 0.575 and 0.918, respectively, suggesting that compound 3 is more likely to penetrate the blood–brain barrier (BBB+) (Table S4, ESI†). Compound 2 has a higher PPB (Plasma Protein Binding) value (90.34) compared to compound 3 (52.59), indicating that 2 may exhibit a stronger binding affinity to plasma proteins, which could influence its therapeutic index (Table S4, ESI†). Additionally, compound 2 shows a higher clearance rate (9.731 mL min⁻¹ kg⁻¹) than compound 3 (8.268 mL min⁻¹ kg⁻¹), suggesting that 2 is cleared from the plasma more rapidly. Both compounds fall into the moderate clearance category (Table S4, ESI†). In summary, compound 3 demonstrates better BBB penetration, while compound 2 has higher plasma protein binding and a slightly higher plasma clearance rate. Table S5 (ESI†) provides SwissADME-computed drug-likeness predictions for the metabolism of compounds 2 and 3, focusing on their interactions with various cytochrome P450 (CYP) enzymes. The table evaluates the likelihood of these compounds acting as inhibitors or substrates for different CYP enzymes, which are crucial for drug metabolism. In CYP1A2, both compounds exhibit a high likelihood of acting as inhibitors and substrates, with compound 2 showing slightly higher values. In contrast, for CYP2C19, compound 2 demonstrates a moderate probability of inhibition, while compound 3 shows a low probability of inhibition, but both compounds are likely to act as substrates. Regarding CYP2D6, compound 3 has a high probability of being both an inhibitor and a substrate, while compound 2 shows moderate inhibition with a higher substrate potential. For CYP3A4, both compounds display a low to moderate probability of being inhibitors and substrates. In CYP2B6, both compounds are likely inhibitors, but their substrate potential is minimal. Finally, in CYP2C8, both compounds exhibit a low probability of inhibition.

On the other hand, Table S6 (ESI†) outlines the ADMET (absorption, distribution, metabolism, excretion, and toxicity) and pharmacokinetic features of the best-hit compounds, 2 and 3. The table evaluates key properties related to their drug-likeness and potential for human use. Regarding Caco-2 permeability, both compounds exhibit moderate permeability, with compound 2 (−4.999) showing slightly higher values than compound 3 (−4.688), suggesting their potential to cross intestinal cells. In the MDCK assay, both compounds demonstrate low permeability, indicating limited passive diffusion across cellular barriers. For human intestinal absorption (HIA), both compounds show a 0.0 probability of being highly absorbed in the human intestine, pointing to poor absorption. In terms of bioavailability, for F20%, compound 2 (0.051) has a slightly higher probability of having less than 20% bioavailability compared to compound 3 (0.004), suggesting that neither compound is likely to achieve high bioavailability. For F30%, compound 2 (0.887) exhibits a higher probability of having less than 30% bioavailability compared to compound 3 (0.573), indicating that both compounds are likely to display low bioavailability, with compound 2 showing a more pronounced tendency. In summary, the data suggest that both compounds, 2 and 3, possess limited permeability, poor absorption, and low bioavailability, which may significantly affect their efficacy as oral drugs.

Table S7 (ESI†) presents the predicted toxicity profiles of compounds 2 and 3, calculated using OSIRIS Property Explorer and PreADMET. The table evaluates various toxicity endpoints, providing insight into the potential hazards of these compounds. For LC50DM (48-hour Daphnia magna LC50), both compounds show similar toxicity to aquatic organisms, with compound 2 (4.524) being slightly more toxic than compound 3 (4.482). In terms of hepatotoxicity, both compounds have a moderate probability of causing liver toxicity, with compound 2 (0.618) slightly higher than compound 3 (0.606). Regarding carcinogenicity, both compounds also show a moderate probability of being carcinogenic, with compound 3 (0.594) exhibiting a slightly higher likelihood than compound 2 (0.561). For Ames toxicity, both compounds are likely to be mutagenic (Ames positive), with compound 2 (0.807) showing a slightly higher probability than compound 3 (0.789). In hematotoxicity, compound 2 (0.542) has a higher probability of causing blood-related toxicity compared to compound 3 (0.348). In respiratory toxicity, compound 3 (0.682) has a higher probability of being a respiratory toxicant compared to compound 2 (0.628). In summary, both compounds exhibit moderate to high toxicity across various endpoints, with compound 2 showing slightly higher hepatotoxicity and hematotoxicity, while compound 3 has a higher likelihood of being a respiratory toxicant and carcinogen. These toxicity predictions are crucial for assessing the safety profiles of compounds 2 and 3.⁶⁹

3.8. Molecular docking against 1CIN

A molecular docking study was conducted to investigate the binding modes of the tested cycloadducts (2 and 3) with the primary protease of 1CIN, aiming to elucidate their mechanism of action (Fig. 14). These proteins, critical for viral protein production, were selected as targets due to their crucial role in the viral life cycle. Targeting these proteins through therapeutic interventions offers potential benefits for virus elimination. As shown in Table 7, the co-crystallized ligand interacted with the protein's active site by forming five hydrogen bonds, resulting in a free binding energy of −11.00 kcal mol⁻¹. In comparison, compound 2 exhibited a lower binding affinity, with a free binding energy of −9.64 kcal mol⁻¹. Compound 3 demonstrated an even lower binding affinity, with a free binding energy of −7.98 kcal mol⁻¹. These findings are visually represented in Fig. 14.

Fig. 14 Three-dimensional (3D) and two-dimensional (2D) interactions of compounds 2 and 3 with the 1CIN protein.

Table 7 Docking affinity of compounds 2 and 3 as well as the co-crystal ligands for 1CIN, given in kcal mol⁻¹

Compound	Affinity (kcal mol^−1)
2	−9.64
3	−7.98
Co-crystal ligand	−11

---

Based on the molecular docking results, compound 2 showed a stronger binding energy to 1CIN compared to the control molecules, suggesting its potential as a more effective inhibitor. While 3 showed a comparatively lower binding affinity, it remains a promising candidate for inhibiting the target protease.

4. Conclusions

The mechanisms of the [3+2] intramolecular cycloaddition reactions through two stereoisomeric pathways (endo and exo) have been studied. Only one TS was found for each reaction path, implying a one-step asynchronous mechanism. Analysis of the relative Gibbs free energies indicated that the reaction path associated with the formation of the exo-isomer is thermodynamically and kinetically preferred, in agreement with experimental outcomes.

For a deeper understanding of the bond forming/breaking during the chemical reaction mechanism, we performed BET analysis along this intramolecular cycloaddition. The reaction mechanism is globally the same for the two paths and can be summarized as follows: (i) the V(N2,C1) basin population drops with the appearance of the V(N2) basin associated with the N lone pair, (ii) then, two monosynaptic basins (V(C1) and V(C4)) arise simultaneously together with a fall of the V(C1,N2) and V(C4,C5) basin populations, (iii) these two basins then merge to form the C1–C4 bond, and (iv) finally, the V(O3,C5) basin is created from the merger of V(C5) and V′(O3) basins formed in the meantime from the drop of V(O3) and V(C4,C5) basin populations. The molecular docking study revealed that compound 2 demonstrated a stronger binding affinity to the 1CIN protease, suggesting its potential as a more effective inhibitor compared to the control molecules. While compound 3 showed a lower binding affinity, it remains a promising candidate for targeting the protease in viral therapeutic interventions. The compounds 2 and 3 meet all the required filtering criteria, including Lipinski's Rule of Five, suggesting their potential as drug-like candidates with no bioavailability issues. Their physicochemical properties, as illustrated by the radar charts and pharmacokinetic analysis, further support their suitability for drug development.

Author contributions

Mohamed Chellegui, Sofiane Benmetir, Raad Nasrullah Salih: writing, investigation, validation, methodology; Haydar A. Mohammad-Salim, Jesus Vicente De Julian Ortiz: validation, supervision; Ali Ben Ahmed: editing, reviewing, and supervision. All authors have read and approved the published version of the manuscript.

Data availability

Publicly available datasets were analyzed in this study.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

M. C. is grateful to UNamur for granting him the status of scientific collaborator. J. V. deJ.-O. acknowledges the financial support provided by the Universitat de València's special actions program, reference ComPPETE UV-INV-AE-3677056. M. C. also acknowledges the computational resources provided by the “Consortium des Équipements de Calcul Intensif (CÉCI)” (https://www.ceci-hpc.be), including those of the “UNamur Technological Platform of High-Performance Computing (PTCI)” (https://www.ptci.unamur.be). We gratefully acknowledge the financial support from the FNRS-FRFC, the Walloon Region, and the University of Namur (Conventions No. U.G006.15, U.G018.19, U.G011.22, RW1610468, RW/GEQ2016, RW1117545, and RW2110213).

References

1. K. V. Gothelf and K. A. Jørgensen, Chem. Rev., 1998, 98, 863.
2. P. DeShong, J. M. Leginus and S. W. Lander Jr, J. Org. Chem., 1986, 51, 574.
3. T. Kametani, S. D. Chu and T. Honda, J. Chem. Soc., Perkin Trans. 1, 1988, 1593.
4. A. Dondoni and A. Marra, Chem. Rev., 2000, 100, 4395.
5. M. E. Filkina, D. N. Baray, E. K. Beloglazkina, Y. K. Grishin, V. A. Roznyatovsky and M. E. Kukushkin, Int. J. Mol. Sci., 2023, 24, 1289.
6. M. O. Zouaghi, M. Y. Doggui and Y. Arfaoui, J. Mol. Graphics, 2021, 107, 107960.
7. R. Arhin, I. Ofori, A. Fosu, R. Tia, E. Adei and A. Aniagyei, J. Mol. Graphics, 2023, 121, 108451.
8. A. I. Adjieufack, V. Liegeois, I. Ndassa Mboumbouo, J. Ketcha Mbadcam and B. Champagne, J. Phys. Chem. A, 2018, 122, 7472.
9. E. Dresler, A. Wróblewska and R. Jasiński, Molecules, 2024, 29, 3042.
10. L. R. Domingo, B. Mondal, H. A. Mohammad-Salim and N. Acharjee, ChemPhysChem, 2025, e202401106.
11. R. Boutiddar, K. Abbiche, M. D. Mellaoui, A. Imjjad, M. Alahiane, Y. Ait Albrimi and M. Hochlaf, J. Comput. Chem., 2024, 45, 284.
12. R. Boutiddar, K. Abbiche, S. Izzaouihda, A. A. Taleb, N. El Baraka, M. Errai and M. Hochlaf, Comput. Theor. Chem., 2024, 114665.
13. L. R. Domingo, Molecules, 2016, 21, 1319.
14. A. Mondal, H. A. Mohammad-Salim and N. Acharjee, SciRad, 2023, 2, 75.
15. A. Mondal and N. Acharjee, Chem. Heterocycl. Compd., 2023, 59, 145.
16. X. Krokidis, S. Noury and B. Silvi, J. Phys. Chem. A, 1997, 101, 7277.
17. A. D. Becke and K. E. Edgecombe, J. Chem. Phys., 1990, 92, 5397.
18. L. R. Domingo, M. Ríos-Gutierrez, E. Chamorro and P. Perez, ChemistrySelect, 2016, 1, 6039.
19. M. Chellegui, M. Trabelsi, B. Champagne and V. Liegeois, ACS Omega, 2025, 10, 833.
20. M. Chellegui, A. I. Adjieufack, M. Trabelsi, V. Liegeois and B. R. Champagne, ChemPhysChem, 2025, e202400896.
21. K. Fukui, Acc. Chem. Res., 1981, 14, 363.
22. A. I. Adjieufack, V. Liegeois, I. M. Ndassa and B. Champagne, RSC Adv., 2021, 11, 10083.
23. A. I. Adjieufack, B. Champagne and V. Liegeois, Synthesis, 2023, 2070.
24. J. Andres, V. S. Safont, M. Oliva, K. L. Caster and F. Goulay, Int. J. Quantum Chem., 2022, 122, e26892.
25. M. Ríos-Gutierrez, P. Perez and L. R. Domingo, RSC Adv., 2015, 5, 58464.
26. R. Jasiński and E. Dresler, Organics, 2020, 1, 49.
27. H. Mohammad-Salim, J. V. de Julián-Ortiz, K. A. Dahlous, M. S. Islam, T. M. Almutairi and S. Benmetir, Struct. Chem., 2025, 36, 339.
28. H. A. Mohammad-Salim, N. Acharjee and H. H. Abdallah, Theor. Chem. Acc., 2021, 140, 1.
29. S. A. M. Salih, H. A. Basheer, J. V. de Julián-Ortiz and H. A. Mohammad-Salim, Int. J. Mol. Sci., 2023, 24, 9102.
30. A. Tawiah, G. B. Pipim, R. Tia and E. Adei, J. Mol. Model., 2021, 27, 1.
31. W. Dehaen and A. Hassner, J. Org. Chem., 1991, 56, 896.
32. U. Atmaca, R. Saglamtas, Y. Sert, M. Çelik and İ. Gülçin, ChemistrySelect, 2023, 8, e202204578.
33. G. Lotfy, H. El Sayed, M. M. Said, Y. M. A. Aziz, A. Al-Dhfyan, A. M. Al-Majid and A. Barakat, J. Photochem. Photobiol., B, 2018, 180, 98.
34. J. Panda, B. P. Raiguru, S. Nayak, S. Mohapatra, S. Mohapatra, S. Rout and C. R. Sahoo, New J. Chem., 2024, 48, 13609.
35. A. O. Aitouna, H. Mohammad-Salim, A. Zeroual, A. Syed, A. H. Bahkali and J. V. de Julián-Ortiz, Comput. Theor. Chem., 2023, 1228, 114283.
36. A. D. Becke, J. Chem. Phys., 1993, 98, 5648.
37. L. R. Domingo, M. J. Aurell and P. Perez, Tetrahedron, 2015, 71, 1050.
38. L. R. Domingo, K. Kula and M. Ríos-Gutierrez, Eur. J. Org. Chem., 2020, 5938.
39. M. Ríos-Gutierrez, A. Barakat and L. R. Domingo, Organics, 2022, 3, 122.
40. L. R. Domingo, RSC Adv., 2014, 4, 32415.
41. A. E. Reed, L. A. Curtiss and F. Weinhold, Chem. Rev., 1988, 88, 899.
42. C. Gonzalez and H. B. Schlegel, J. Phys. Chem., 1990, 94, 5523.
43. K. Fukui, J. Phys. Chem., 1970, 74, 4161.
44. J. Tomasi, B. Mennucci and R. Cammi, Chem. Rev., 2005, 105, 2999.
45. B. Silvi and A. Savin, Nature, 1994, 371, 683.
46. B. Silvi, J. Mol. Struct., 2002, 614, 3.
47. S. Noury, X. Krokidis, F. Fuster and B. Silvi, Comput. Chem., 1999, 23, 597.
48. V. Liegeois, DrawMol; UNamur, https://www.unamur.be/drawmol, 2023.
49. V. Liegeois, DrawProfile; UNamur, https://www.unamur.be/drawprofile, 2023.
50. M. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci and G. A. Petersson, et al., Gaussian 16, Revision A. 03, Gaussian Inc., Wallingford, CT, 2016.
51. R. G. Parr and W. Yang, Density Functional Theory of Atoms and Molecules, Oxford University Press, New York, 1989.
52. L. R. Domingo, M. Ríos-Gutierrez and P. Perez, Molecules, 2016, 21, 748.
53. L. R. Domingo, M. J. Aurell, P. Perez and R. Contreras, Tetrahedron, 2002, 58, 4417.
54. L. R. Domingo and P. Perez, Org. Biomol. Chem., 2011, 9, 7168.
55. M. Ríos-Gutierrez, A. Saz Sousa and L. R. Domingo, J. Phys. Org. Chem., 2023, 36, e4503.
56. M. J. Aurell, L. R. Domingo, P. Perez and R. Contreras, Tetrahedron, 2004, 60, 11503.
57. L. R. Domingo, P. Perez and J. A. Sáez, RSC Adv., 2013, 3, 1486.
58. L. R. Domingo and J. A. Sáez, Org. Biomol. Chem., 2009, 7, 3576.
59. P. K. Chattaraj, S. Duley and L. R. Domingo, Org. Biomol. Chem., 2012, 10, 2855.
60. L. R. Domingo, M. J. Aurell and P. Perez, Tetrahedron, 2014, 70, 4519.
61. E. Reed, R. B. Weinstock and F. Weinhold, J. Chem. Phys., 1985, 83, 735.
62. E. Reed, L. A. Curtiss and F. Weinhold, Chem. Rev., 1988, 88, 899.
63. R. Jasiński, Tetrahedron Lett., 2015, 56, 532.
64. L. R. Domingo and M. Ríos-Gutierrez, SciRad, 2023, 2, 1.
65. R. F. Bader, Acc. Chem. Res., 1985, 18, 9.
66. R. F. Bader and H. Essen, J. Chem. Phys., 1984, 80, 1943.
67. M. Sadowski and K. Kula, Molecules, 2024, 29, 5066.
68. M. Sadowski, B. Synkiewicz-Musialska and K. Kula, Molecules, 2024, 29, 542.
69. D. A. Filimonov, A. A. Lagunin, T. A. Gloriozova, A. V. Rudik, D. S. Druzhilovskii, P. V. Pogodin and V. V. Poroikov, Chem. Heterocycl. Compd., 2014, 50, 444.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5nj01058f

---