How the mechanism of a [3 + 2] cycloaddition reaction involving a stabilized N-lithiated azomethine ylide toward a π-deficient alkene is changed to stepwise by solvent polarity? What is the origin of its regio- and endo stereospecificity? A DFT study using NBO, QTAIM, and NCI analyses

Abstract

A theoretical study at the MPWB1K/6-31G(d) level was performed on the [3 + 2] cycloaddition (32CA) reaction of a stabilized N-lithiated azomethine ylide (AZY 1) toward a π-deficient methyl acrylate (MA 11) in tetrahydrofuran (THF) to elucidate the regio- and stereoselectivity as well as the mechanistic aspects involved in this reaction. The energetics results, in excellent agreement with the experimental findings, indicate that between the two competitive C3–C4 and C3–C5 regioisomeric channels, the former is both kinetically and thermodynamically preferred over the latter, in which the electrophilic attack of the C4 carbon atom of MA 11 on the C3 carbon atom of AZY 1 through an endo approach mode passes the endo transition state TS1n, whose Gibbs free energy is 9.8 kcal mol⁻¹ higher than those of the separate reagents. The analysis of the intrinsic reaction coordinate (IRC) curves of TS1n in the gas phase and in THF demonstrates that unlike the gas phase, the noticeable charge separation caused by the electron density delocalization from the carbonyl oxygen atom of the MA 11 moiety toward the Li cation of the AZY 1 moiety is sufficiently stabilized by the solvent polarity, leading to the formation of the zwitterionic intermediate IN1n according to a stepwise mechanism. Upon the formation of IN1n, during which the first C3–C4 single bond forms, the attractive electrostatic force between the oppositely charged C1 (+0.08e) and C5 (−0.51e) carbon atoms promotes the ring closure step via TS11n with a non-appreciable barrier of 1.3 kcal mol⁻¹, leading to the formation of a second C1–C5 single bond at the corresponding lithiated cycloadduct. The calculated electrophilic and nucleophilic Parr functions at the reactive sites of the reagents together with the presence of large steric hindrances along the unfavourable C3–C5 regioisomeric channel enables us to explain the C3–C4 regiospecificity observed experimentally. In addition, the origin of the endo-stereospecificity experimentally displayed by the 32CA reaction under consideration can be nicely portrayed using non-covalent interactions (NCIs) analysis at the exo TS1x and endo TS1n transition states involved in the preferred C3–C4 regioisomeric channel.

---

1. Introduction

The [3 + 2] cycloaddition (32CA) reaction, one of the most well-known and useful reactions in the chemist's arsenal, is being increasingly used in the vast majority of areas of chemistry, such as materials chemistry,¹ drug discovery,² and biological chemistry.³ The 32CA reaction occurs when a three-atom-component (TAC) is attacked by an unsaturated bond to generate five-membered heterocyclic compounds in a highly regio- and stereoselective fashion, provided that adequate functional groups are substituted on the skeletons of the reagents.⁴

The electronic nature of the TACs involved in the 32CA reaction allows its classification into two different types: (i) pseudodiradical-type (pr-type) and (ii) zwitterionic-type (zw-type).⁵ While in pr-type 32CA reactions, the very high pseudodiradical (singlet diradical) character of the involved TAC acts as a driving force to promote the reaction via a non-polar early transition state (TS), the highly zwitterionic character of the given TAC acts as a driving force in zw-type 32CA reactions, promoting the reaction through a polar and relatively late TS.^5,6

Ess and Houk⁷ have explained the reactivity of TACs toward 32CA reactions in terms of a distortion/interaction model. In Ess and Houk's proposed model, whose applicability was examined in the 32CA reaction of well-known TACs, diazonium, nitrilium, and azomethine betaines, toward ethylene and acetylene, the activation energy (ΔE^≠) is divided into two additive terms: distortion energy (ΔE^≠_d) and interaction energy (ΔE^≠_i). While the distortion energy is associated with the energy required to distort reagents from their ground state geometries to the corresponding geometries in the TS, the interaction energy corresponds to the interaction between two fragments in the TS. The relationship between the activation, distortion, and interaction energies is graphically depicted in Scheme 1. Ess and Houk found a reasonable linear correlation coefficient of R² = 0.97 between the B3LYP/6-31G(d) activation enthalpies and distortion energies in the investigated 32CA reactions and concluded that the different reactivity of the studied TACs mainly arises from the distortion contribution of TACs in the distortion energy, i.e. ΔE^≠_d = ΔE^≠_d (TAC) + ΔE^≠_d (ethylene/acetylene) ≈ ΔE^≠_d (TAC).

Scheme 1 A schematic relationship between the activation, ΔE^≠, distortion, ΔE^≠_d, and interaction, ΔE^≠_i, energies in the 32CA reaction of different TACs toward ethylene and acetylene, investigated by Ess and Houk.⁷

Despite the good correlation coefficient found between the activation enthalpies and distortion energies, it appears that the Ess and Houk model should be revisited because this model fails when is employed to explain why the activation energy depends on the geometry, or vice versa.^6b Indeed, Domingo has recently indicated “while the distribution of the electron density is responsible for the molecular shape and physical properties, the capability for changes in electron density is responsible for the reactivity”.⁸ In other words, dividing the TS geometry into two separate fragments, as proposed by Ess and Houk, does not make any physical sense because it is well known that through density functional theory (DFT), the energy of a given system, such as a TS, is not only the functional of the electron-density distribution over the whole of the system but highly depends on the spatial position of the nuclei. In consequence, the energy of the two separate fragments cannot be correlated with the energy of the corresponding TS, in consideration of the very important fact that “the external potential created by one fragment over the other one is lost when one of them is removed”.⁹

Recently, Braida and co-workers established a reactivity model for the TACs studied by Ess and Houk in terms of their singlet-diradical character at their ground state geometries.¹⁰ Using the “breathing-orbital valence bond (BOVB)” approach, they calculated the singlet-diradical weight of nine TACs participating in the 32CA reaction with ethylene and/or acetylene that were considered in the Ess and Houk study. They found that the reactivity of TACs is well correlated with their singlet-diradical character when the barriers of the corresponding 32CA reactions toward ethylene and/or acetylene, calculated by Ess and Houk, are plotted against singlet-diradical weights. Meanwhile, in the Braida study, the uncommon and very difficult-to-use Xiamen VB (XMVB) package is employed; very recently, we have calculated the singlet-diradical character for azomethine betaines (AZBs) using the very common, user friendly, and familiar GAUSSIAN and Gen NBO.05W software packages.^6b Some high linear correlation coefficients were found between the singlet-diradical character of the AZBs and the corresponding MP4(SDTQ)/6-311++G(d,p) activation energies toward the 32CA reaction with ethylene, indicating that this singlet-diradical character is responsible for the high reactivity of AZBs in 32CA reactions. The singlet-diradical weights of the AZBs (CH₂ = NH⁺–X⁻, X = CH₂, NH, and O) among five feasible Lewis resonance structures, (I) through (V), calculated by Braida using the BOVB method¹⁰ and via natural resonance theory (NRT) employing Gen NBO.05W,^6b as well as the computed activation barriers of the corresponding 32CA reactions toward ethylene, are given in Scheme 2. When the CBS-QB3 and MP4(SDTQ)/6-311++G(d,p) calculated activation energies associated with the 32CA reaction of AZBs toward ethylene are plotted against the corresponding singlet-diradical weights calculated by the BOVB methodology employed by Braida, as well as the NRT approach (see values presented in Scheme 2), the resulting diagram, as shown in Scheme 3, contains some appealing points.

Scheme 2 The singlet-diradical weights of AZBs among feasible five Lewis resonance structures, (I) through (V), computed by Braida using the BOVB method and via the NRT approach employing Gen NBO.05W, as well as the calculated activation barriers of the corresponding 32CA reactions toward ethylene.

Scheme 3 Linear correlation between the calculated activation energies corresponding to the 32CA reaction of AZBs toward ethylene and the computed singlet-diradical character (weight) of AZBs at their ground state geometries (see the text for details).

As illustrated in Scheme 3, plotting CBS-QB3 activation energies (taken from the Ess and Houk study)⁷ against the singlet-diradical weights of AZBs (calculated by Braida via the BOVB method)¹⁰ results in a linear correlation coefficient of R² = 0.9646, while a similar plot of the MP4(SDTQ)/6-311++G(d,p) activation energies against the singlet-diradical weights of AZBs calculated via the NRT approach (taken from our recent study)^6b gives a more reliable and better linear correlation coefficient of R² = 0.9858. It is worth noting that this simple comparison between two R² values given in Scheme 3 is not intended to call Braida's highly valuable study into question, as the methodologies employed to obtain the activation energies are quite different in both studies. On the other hand, if the singlet-diradical weights of AZBs calculated through the NRT approach are plotted against those provided by Braida employing the BOVB method (in order to avoid overcrowding and deviation from the main goal of the current study, the corresponding plot is not given here), an excellent linear correlation coefficient of 0.9955 is obtained, emphasizing not only that our recent methodology is in complete harmony with Braida's but also that the NRT approach can be used as a much simpler approach to obtain singlet-diradical weights, at least for AZBs.

Among the studied AZBs, the formoazomethine ylide (X = CH₂ in Scheme 2) is the simplest azomethine ylide, with the largest singlet-diradical character (41.3% calculated by the BOVB method or 23.9% calculated via the NRT approach); thus, it displays the highest reactivity in the 32CA reaction toward ethylene, as verified by the lowest activation energies of 1.2 kcal mol⁻¹ calculated at the CBS-QB3 level and 1.0 kcal mol⁻¹ at the MP4(SDTQ)/6-311++G(d,p) level. Azomethine ylides (AZYs), in which four electrons are delocalized over a continuous C–N–C framework with a planar and allylic structure, are one of the most well-known TACs; their 32CA reactions toward C–C double bonds furnish pyrrolidine derivatives as important building blocks in natural products synthesis, pharmaceuticals and chiral ligands.¹¹

AZYs can be prepared through several common procedures, such as¹² (i) proton abstraction from α-iminoesters, whose remarkably acidic α-carbon hydrogens, arising from the coordination of a Lewis acid catalyst to the heteroatoms, can easily be attacked by an appropriate base; (ii) photolysis or thermolysis of aziridines; and (iii) acid catalyzed decomposition of N-alkyl-N-methoxymethyl-N-(trimethylsilyl) methylamines (see Scheme 4). As mentioned, the participation of AZYs in the 32CA reaction toward a C–C double bond can be used as a valuable tool to synthesize polyfunctional pyrrolidine ring systems which are found to be the main structural elements of many alkaloids and pharmacologically active compounds.¹³ Most of the various experimental procedures available to synthesize pyrrolidines via the 32CA reaction of AZYs toward π-deficient alkenes are associated with the use of stabilized N-metalated AZYs, usually formed in situ by coordination of a metal transition salt as a Lewis acid catalyst with the corresponding deprotonated α-iminoester.¹⁴ This method promotes the target 32CA reaction under mild conditions with a high degree of stereo- and regioselectivity. Moreover, among the employed metal cations, it has well been established that silver(I) and lithium(I) metal cations are the most appropriate choices to facilitate the investigated 32CA reaction along with an excess of base, such as a tertiary amine.¹⁵

Scheme 4 In situ generation of AZYs through (i) proton abstraction from an α-iminoester (LA indicates a Lewis acid catalyst); (ii) photolysis or thermolysis of an aziridine; and (iii) acid catalyzed decomposition of an N-alkyl-N-methoxymethyl-N-(trimethylsilyl)methylamine.

On the basis of a theoretical study supported by experimental evidence, the mechanistic aspects of the 32CA reaction of stabilized N-metalated AZYs toward π-deficient alkenes have been clarified by Vivanco and co-workers (see Scheme 5).¹⁶ As depicted in Scheme 5, α-iminoester 1 in the presence of tertiary amine R₃N as a base and a metal salt MX is firstly converted into stabilized N-metalated AZY 2, whose 32CA reaction toward the π-deficient alkene 3 bearing an electron-withdrawing group (EWG) generates a zwitterionic intermediate, 4.

Scheme 5 Mechanism for the 32CA reaction of stabilized N-metalated AZYs toward π-deficient alkenes proposed by Vivanco and co-workers on the basis of a theoretical study supported by experimental evidence.¹⁶

The next step involves a ring closure process at intermediate 4, affording N-metalated formal [3 + 2] cycloadduct 5 which, in turn, is converted into the final product 6 in the presence of an aluminium salt, followed by the recycling of the tertiary amine together with the metal salt. It is worth noting that the presence of a strong EWD group on the skeleton of alkene 3 leads to extra stabilization of the negative charge appearing in intermediate 4; this favours the stepwise mechanism, particularly when a polar solvent is employed.

Note that when a bulky coordinating chiral ligand is employed together with a metal salt, generation of the final product 6 (see Scheme 5) can also occur in a highly enantioselective manner. This enantiomeric achievement can be explained by considering the fact that one of two attacking re- and/or si-faces is significantly shielded by the bulky chiral ligand coordinated to the metal (M) of N-metalated AZY; thus, the π-deficient alkene is allowed to attack the given AZY only along the less hindered re- and/or si-face, leading to the formation of a single enantiomer. For instance, the 32CA reaction of imine 7 with N-methylmaleimide 8 to yield cycloadduct 9 has very recently been studied by Ponce and co-workers both experimentally and theoretically (see Scheme 6).¹⁷

Scheme 6 The 32CA reaction between imine 7 and N-methylmaleimide 8 furnishing cycloadduct 9, studied both experimentally and theoretically by Ponce and co-workers.¹⁷

They experimentally found that the 32CA reaction between imine 7 and N-methylmaleimide 8 in the presence of silver(I) acetate (AgOAc), as the metal salt, and Taniphos, a P,P axially chiral ligand, in methyl tert-butyl ether (MTBE) as the solvent without employment of base affords exo and endo cycloadduct 9 in a good yield of 83% via a highly endo-stereoselective (endo 9/exo 9 > 20/1) and enantioselective (ee = 91%) manner. Using a theoretical approach, they rationalized the high enantioselectivity observed experimentally. The optimized structure of the most stable AZY complex derived from imine 7 and N-methylmaleimide 8 is given in Scheme 7. As shown, when Ag ion is additionally coordinated with Taniphos, the top attacking face is strongly shielded by some parts of the bulky Taniphos ligand so that N-methylmaleimide 8 can more easily attack the AZY via the bottom face, leading to the formation of the corresponding cycloadduct in a highly enantioselective fashion (ee = 91%).

Scheme 7 Optimized structure of the most stable AZY complex derived from imine 7 and N-methylmaleimide 8 in the presence of the metal salt AgOAc and the P,P axially chiral ligand Taniphos provided by Ponce and co-workers, in which hydrogens are omitted for clarity.¹⁷

A recent experimental study on the regio- and endo stereospecific 32CA reaction between AZY 1, in situ generated from α-imino phosphonate 10, and the π-deficient methyl acrylate MA 11 to yield α-amino phosphonate derivative 12 has been investigated by Dondas et al.,¹⁸ in which lithium bromide (LiBr), 1,8-diazabicyclo[5.4.0]undec-7-ene (DBU), and tetrahydrofuran (THF) were employed as the metal salt, base, and solvent, respectively (see Scheme 8).

Scheme 8 The regio- and endo-stereospecific 32CA reaction between AZY 1, in situ generated from α-imino phosphonate 10, and π-deficient methyl acrylate 11 affording α-amino phosphonate derivative 12 experimentally studied by Dondas et al.¹⁸

Herein, the 32CA reaction of AZY 1 toward π-deficient methyl acrylate MA 11 in THF affording α-amino phosphonate derivative 12 (Scheme 8) is theoretically studied at the MPWB1K/6-31G(d) level in order to characterize the energetics, the regio- and endo stereospecificity observed experimentally, and the reaction mechanism of this process.

2. Computational details

While B3LYP¹⁹ is the most popular density functional and is widely used in the computational chemistry community, several works have shown that the reaction exothermicity is underestimated by this functional.²⁰ Truhlar's group has proposed some hybrid meta-density functionals, such as M06-2X and MPWB1K.²¹ While overestimation of the activation energies by the M06-2X functional has very recently been proved by Jasiński in a [4 + 2] cycloaddition reaction,²² the MPWB1K functional gives good results for thermochemistry and thermochemical kinetics as well as excellent saddle point geometries.²¹ Consequently, in the present study, DFT computations were carried out by means of the MPWB1K exchange–correlation functional, together with the standard 6-31G(d) basis set.²³ Berny analytical gradient optimization method using GEDIIS²⁴ was employed in the geometry optimization steps. The stationary points were characterized by frequency calculations in order to verify that the TSs have one and only one true imaginary frequency along the reaction coordinates. The intrinsic reaction coordinate (IRC) paths²⁵ were traced in order to check the energy profiles connecting each TS to the two associated minima of the proposed mechanism using the second order González–Schlegel integration method.²⁶

The solvent effects of THF (ε = 7.42) were taken into account through re-optimization of the gas phase stationary points using the polarizable continuum model (PCM) developed by Tomasi's group²⁷ in the framework of the self-consistent reaction field (SCRF).²⁸ Then, in order to explore whether the kinetics and thermodynamics as well as the selectivities of the studied 32CA reaction can be affected by diffuse functions and to verify the validity of the employed 6-31G(d) basis set, a single point-energy calculation was performed over the optimized geometries using the 6-31+G(d) basis set.²³ In this way, considering the statistical formulation of thermodynamics²⁹ and applying some appropriate manipulations, the corrected Gibbs free energy, G(corrected), of the stationary points in THF can simply be obtained using the following equation:

G(corrected) = {(U_T − U₀) + ZPVE}_{MPWB1K/6-31G(d)} + E^ele_{MPWB1K/6-31+G(d)} + T{R − (S)_{MPWB1K/6-31G(d)}} (1)

The first term in the right-hand side of eqn (1) is the thermal correction to energy, in which U_T, U₀, and ZPVE indicate the internal energy at T kelvin, the internal energy at zero kelvin, and the zero-point vibrational energy (zero-point correction), respectively, attained through MPWB1K/6-31G(d) frequency calculations. The second term, E^ele_{MPWB1K/6-31+G(d)}, is the total electronic energy taken from the MPWB1K/6-31+G(d) single-point energy calculation over the MPWB1K/6-31G(d) optimized geometry, and the last term includes the temperature T (298 K is taken into account according to the experimental conditions of the studied reaction), the universal gas constant R, and the entropy S obtained via MPWB1K/6-31G(d) frequency calculations.

The electronic structures of the stationary points were analyzed by the natural bond orbital (NBO) method.³⁰ Non-covalent interactions (NCIs) were computed by evaluating the pro-molecular density using a methodology previously described.^31,32 All computations were carried out with the Gaussian 09 suite of programs.³³

The global electrophilicity index ω³⁴ is given by the following expression, ω = μ²/2η, based on the electronic chemical potential, μ, and the chemical hardness, η. Both quantities may be approached in terms of the one-electron energies of the frontier molecular orbitals, HOMO and LUMO, ε_H and ε_L, as μ ≈ (ε_H + ε_L)/2 and η ≈ (ε_L − ε_H), respectively.³⁵ The global nucleophilicity index N,³⁶ based on the HOMO energies obtained within the Kohn–Sham scheme,³⁷ is defined as N = ε_HOMO(Nu) − ε_HOMO(TCE), in which Nu denotes the given nucleophile. This relative nucleophilicity index refers to tetracyanoethylene (TCE), as it presents the lowest HOMO energy in a long series of molecules already investigated in the context of polar organic reactions. This choice enables convenient handling of a nucleophilicity scale of positive values.

Nucleophilic P⁻_k and electrophilic P⁺_k Parr functions³⁸ were obtained through analysis of the Mulliken atomic spin density (MASD) of the radical cation of AZY 1 and the radical anion of MA 11, respectively.

3. Results and discussion

The present theoretical study is divided into four parts: (i) first, an analysis of the DFT reactivity indices at the ground state of the reagents involved in the 32CA reaction of AZY 1 toward MA 11 is performed in order to predict the reactivity and regioselectivity of the studied cycloaddition reaction; (ii) in the second part, the coordination of Li cation, an electron-density acceptor center, with nitrogen atoms and the oxygen atom of the P=O functional group in AZY 1, as electron-density donating centers, is discussed in terms of electron-density delocalization to obtain a deeper insight into the interactions between the abovementioned centers using NBO together with Bader's quantum theory of atoms-in-molecules (QTAIM) analyses; (iii) in the third part, the potential energy surfaces (PESs) associated with the reaction paths involved in the 32CA reaction between AZY 1 and MA 11 are studied to elucidate the energetics and mechanistic aspects as well as the regio- and stereospecificity of this reaction; and finally, (iv) in the fourth part, the origin of the endo-stereospecificity in the 32CA reaction of AZY 1 toward MA 11 is explored using NCI analysis.

3.1. Analysis of the global and local DFT reactivity indices at the ground state of the reagents involved in the 32CA reaction between AZY and MA 11

Global reactivity indices defined within the conceptual DFT³⁹ are powerful tools to explain the reactivity of cycloaddition reactions. Since the global electrophilicity and nucleophilicity scales are given at the B3LYP/6-31G(d) level, reagents were optimized at this DFT level. The global reactivity indices, namely, electronic chemical potential (μ), chemical hardness (η), global electrophilicity (ω) and global nucleophilicity (N), for AZY 1 and MA 11 are presented in Table 1.

Table 1 B3LYP/6-31G(d) electronic chemical potential, μ, chemical hardness, η, global electrophilicity index, ω, and global nucleophilicity index, N, in eV, for AZY 1 and MA 11

	μ	η	ω	N
AZY 1	−2.64	3.32	1.05	4.81
MA 11	−4.31	6.22	1.49	1.70

---

As shown in Table 1, the higher electronic chemical potential of AZY 1, −2.64 eV, than that of MA 11, −4.31 eV, predicts that the corresponding 32CA reaction global electron density transfer (GEDT),⁴⁰ as a measure of reaction polarity, will occur from AZY 1 toward MA 11, in good agreement with the GEDT analysis performed at the exo TS1x (0.29e) and endo TS1n (0.27e) transition states along the favourable C3–C4 regioisomeric channel (see later). These GEDT values indicate that the 32CA reaction of AZY 1 toward MA 11 does not display a highly polar character. AZY 1 has a low global electrophilicity ω index, 1.05 eV, and a high nucleophilicity N index, 4.81 eV, being classified as a weak electrophile and a strong nucleophile within the electrophilicity⁴¹ and nucleophilicity⁴² scales, respectively. On the other hand, MA 11 exhibits a low global electrophilicity ω index, 1.49 eV, and a low nucleophilicity N index, 1.70 eV, being classified as a weak electrophile and a weak nucleophile. Analysis of these global indices indicates that in a low-polarity 32CA reaction, AZY 1 and MA 11 will act as the nucleophile and electrophile, respectively.

When an electrophile–nucleophile pair approach each other, the most favourable reactive channel is that associated with the initial two-center interaction between the most electrophilic center of the electrophile and the most nucleophilic center of the nucleophile. Recently, Domingo has proposed the nucleophilic P⁻_k and electrophilic P⁺_k Parr functions³⁸ derived from the excess of spin electron-density reached via the GEDT process from the nucleophile to the electrophile as a powerful tool in the study of the local reactivity of polar processes. Since AZY 1 and MA 11, respectively, act as the nucleophile and electrophile in the studied 32CA reaction, the nucleophilic P⁻_k Parr functions on the reactive sites of AZY 1 and the electrophilic P⁺_k Parr functions on the reactive site of MA 11 were calculated and are depicted in Fig. 1.

Fig. 1 MPWB1K/6-31G(d) map of the Mulliken atomic spin density (MASD) of the radical cation AZY 1˙⁺ and the nucleophilic P⁻_k Parr function on the reactive sites of AZY 1; map of the MASD of the radical anion MA 11˙⁻ and the electrophilic P⁺_k Parr function on the reactive sites of MA 11.

Considering the calculated Parr functions given in Fig. 1, while the C3 carbon atom in AZY 1 is the most nucleophilic center, possessing a nucleophilic Parr function of 0.51, the most electrophilic center in MA 11 is the C4 carbon atom, possessing an electrophilic Parr function of 0.58 (see Scheme 8 for atom numbering). Therefore, it is predicted that the most favourable electrophile–nucleophile interaction along the electrophilic attack of MA 11 on AZY 1 should occur between the most electrophilic center of MA 11, the C4 carbon atom, and the most nucleophilic center of AZY 1, the C3 carbon atom. While the electronic effects characterized by the calculated Parr functions predict a C3–C4 regioselectivity in the considered 32CA reaction, the C3–C4 regiospecificity observed experimentally (see Scheme 8)¹⁸ cannot be decisively related to electronic effects alone because the nucleophilic Parr functions calculated at the reactive sites of AZY 1 (see the left illustration in Fig. 1) present an extremely low (even negligible) difference of 0.04. Indeed, geometrical assessments of TS1n (involved in the C3–C4 regioisomeric channel) and of TS2n (involved in the C3–C5 regioisomeric channel) indicates that the experimentally observed C3–C4 regiospecificity is mainly caused by some noticeable steric hindrances between MA 11 and AZY 1 fragments of TS2n along the unfavourable C3–C5 regioisomeric approach mode. Consequently, it seems reasonable to indicate that electronic effects and, more importantly, steric hindrances are responsible for the C3–C4 regiospecificity in the studied 32CA reaction.

3.2. NBO and QTAIM analyses of electron-density delocalization between heteroatoms and lithium cation in AZY 1

As depicted in Scheme 8, the Li cation in AZY 1 interacts with two nitrogen atoms and an oxygen atom. Which orbitals are involved in these interactions? How strong are these interactions? Which of these interactions is the strongest? Answering such questions can be attractive for every chemist, particularly when the abovementioned interactions are graphically portrayed. In this way, one of the most important concepts in chemistry, namely electron-density delocalization, is visualized, particularly to lead experimentalists to a deeper insight into the chemical interactions between a metal cation and a heteroatom. In order to attain this purpose, in this section, the electron-density delocalization between Li cation and the heteroatoms in AZY 1 is studied using NBO and QTAIM analyses.

Fig. 2 displays the shape of the donor pre-orthogonal natural bond orbitals (PNBOs)⁴³ associated with heteroatoms and of the acceptor PNBO corresponding to the Li cation in AZY 1, obtained via NBO computation. Analysis of the NBO results indicates that LP(N2) is located in an sp^4.57 hybridized NBO mainly built from the valence 2S, 2P_x, 2P_y, and 2P_z natural atomic orbitals (NAOs) of the N2 nitrogen atom, in which the largest coefficient of 2P_y NAO, 0.70, causes LP(N2) to display an orientation mainly along the Y-axis. On the other hand, LP1 (O7) is located in an sp^0.65 hybridized NBO in which the largest contribution of 0.77 corresponds to the valence 2S NAO of the O7 oxygen atom, which explains why this hybridized NBO highly resembles an S-shaped orbital which is mainly extended along the Y-axis due to the largest coefficient of 0.42, displayed by the Y-component of the 2P NAOs of O7.

Fig. 2 MPWB1K/6-31G(d) shapes of the donor PNBOs LP(N2), LPs(O7), and LP(N9) and of the acceptor LP*1(Li8), in which LP and LP* denote filled Lewis type NBOs associated with the lone pair of the corresponding heteroatom and the vacant non-Lewis type NBO of the Li cation, respectively.

Similarly, LP2(O7) is an sp^99.99 hybridized NBO (almost a pure P atomic orbital) with contributions of the 2P_x, 2P_y, and 2P_z NAOs of the O7 oxygen atom of 0.77, 0.48, and 0.41, respectively. The LP(N9) NBO has an sp^2.91 hybridization mainly built from the valence 2S, 2P_x, 2P_y, and 2P_z NAOs of the N9 nitrogen atom, with contributions of 0.50, 0.63, 0.57, and 0.16, respectively. LP*1(Li8) NBO, the only detectable vacant natural bond orbital by the NBO computation performed over the optimized structure of AZY 1, which participates in the interaction with heteroatoms, is the valence 2S NAO of Li cation, with a contribution of 0.99 to the corresponding s^99.16 hybridized NBO.

According to Weinhold's NBO calculation, hyperconjugation has a stabilizing effect arising from the delocalization of electron density from filled (bonding or lone pair) Lewis type NBOs to neighboring vacant orbitals (non-Lewis type NBOs, such as anti-bonding or Rydberg) when these orbitals are properly oriented. For each donor NBO (i) and acceptor NBO (j), the stabilization energy can be described by means of the second-order perturbation interaction energy, E⁽²⁾, and evaluated by the following expression:

where q_i is the donor orbital occupancy, ε_i and ε_j are diagonal elements (orbital energies) and F_i,j is the off-diagonal NBO Fock matrix elements.⁴³ Table 2 summarizes the shape of the donor NBO (i)–acceptor NBO (j) interactions together with the corresponding E⁽²⁾ values in AZY 1. As presented in the second column of Table 2, while some parts of the constructive overlap between the LP1(O7) and LP*1(Li8) NBOs with a positive sign, green–green overlap, are canceled by the destructive overlap of the NBOs with an opposite sign, i.e. positive (green)–negative (pink) overlap, there is no destructive overlap in the case of other NBO interactions because the negative LOB of the heteroatom NBO is completely out of reach of the LP*1(Li8) NBO, with a highly spherical positive sign. From the third column of Table 2, it is evident that the greatest donor–acceptor interaction takes place between LP(N9) and LP*1(Li8), with an E⁽²⁾ value of 16.7 kcal mol⁻¹, while the LPs(O7)–LP*1(Li8) interactions, whose stabilizing E⁽²⁾ values have a sum of 20.3 kcal mol⁻¹, permits us to consider the trend of O7–Li8 > N9–Li8 > N2–Li8 for the strength of the investigated donor–acceptor interactions.

Table 2 MPWB1K/6-31G(d) shapes of donor NBO (i)–acceptor NBO (j) interactions together with the corresponding stabilizing E⁽²⁾ values in AZY 1

Donor-NBO (i) → acceptor-NBO (j) interaction	Shape of donor-NBO (i) → acceptor-NBO (j) interaction	E^(2) value (kcal mol^−1)
LP(N2) → LP*1(Li8)		13.5
LP1(O7) → LP*1(Li8)		13.0
LP2(O7) → LP*1(Li8)		7.3
LP1(N9) → LP*1(Li8)		16.7

---

A complementary view to that presented by the NBO computations can be achieved via topology analysis of the electron density in the framework of Bader's QTAIM, in which the sign of the Laplacian of the electron density, ∇²ρ(r), and the electronic energy density, H(r), at the bond critical point (BCP) of a given interaction provide valuable information about the strength of an interaction. On the other hand, the magnitude of the −G(r)/V(r) ratio, in which the numerator and denominator indicate the kinetic and potential electron energy densities at the BCP of a given interaction, respectively, is used to determine the magnitude to which an interaction is covalent in nature. In this way, weak interactions exhibit positive values for both ∇²ρ(r) and H(r), while medium interactions display positive ∇²ρ(r) but negative H(r).

Similarly, in strong interactions, both ∇²ρ(r) and H(r) are negative. When the −G(r)/V(r) ratio is greater than one, the given interaction should be classified as non-covalent, while for 0.5 < −G(r)/V(r) < 1, the interaction becomes partly covalent.⁴⁴ The MPWB1K/6-31G(d) generated wave function of optimized AZY 1 was analyzed using AIM2000,⁴⁵ and the corresponding results are sketched in Fig. 3.

Fig. 3 MPWB1K/6-31G(d) molecular graph of AZY 1, in which the BCPs of N9–Li8, N2–Li8, and O7–Li8 are denoted A, B, and C, respectively; some topological descriptors are given in the corresponding boxes (see the text for details).

As shown in Fig. 3, not only are the ∇²ρ(r) and H(r) topological descriptors both positive, but the value of the −G(r)/V(r) descriptor is also greater than one at BCPs A, B, and C, implying that each Li–heteroatom pair interaction in AZY 1 should be classified as weak and non-covalent.

3.3. Study of the PESs associated with the reaction paths involved in the 32CA reaction of AZY 1 toward MA 11 in THF

Due to the asymmetry of both AZY 1 and MA 11, four reactive channels can compete in the corresponding 32CA reaction. These reactive channels correspond to the two stereoisomeric approach modes of the carbonyl moiety of MA 11 relative to the Li(I) cation of AZY 1, named exo, in which the carbonyl group and the Li(I) cation are located on two opposite sides, and endo, in which the carbonyl group is located over the Li(I) cation; and the two regioisomeric approach modes of the C3 carbon atom of AZY 1 toward the C4 (via C3–C4 attack) or C5 (via C3–C5 attack) carbon atoms of MA 11. Exploration of the IRC curves of the corresponding TSs evidently demonstrates that two different mechanistic aspects are portrayed by the 32CA reaction of AZY 1 toward MA 11 in THF. Interestingly, while the exo stereoisomeric approach follows a one-step mechanism along both the C3–C4 and C3–C5 regioisomeric pathways, the endo stereoisomeric approach follows a stepwise mechanism, accompanied with the formation of corresponding zwitterionic intermediates. It is worth mentioning that under thermal conditions, some ylides react with π-deficient ethene derivatives via a stepwise mechanism involving zwitterionic species.⁴⁶ In addition, a similar mechanism has also been observed in the case of 32CA reactions involving other allylic-type TACs.⁴⁷

As depicted in Scheme 9, six transition state structures, TS1x, TS1n, TS11n, TS2x, TS2n, and TS22n, two stable intermediates, IN1n and IN2n, and four lithiated formal [3 + 2] cycloadducts, CA1x, CA1n, CA2x, and CA2n, in which the “x” and “n” notations indicate exo and endo, respectively, were located and characterized on the PES of the investigated 32CA reaction in THF. The MPWB1K/6-31G(d) calculated relative Gibbs free energies are given in Scheme 9; the values of G(corrected) calculated according to eqn (1), followed by the corresponding energy profiles and a detailed explanation about how the inclusion of diffuse functions leads to results that cannot properly justify the experimental outcomes, are presented in Table S1 and Scheme S1,† respectively, of the ESI.† A similar justification of the failure of diffuse functions to modify selectivities has been very recently provided by Nacereddine and co-workers through an interesting theoretical investigation on the 32CA reaction between (Z)-C-phenyl-N-methylnitrone and dimethyl 2-benzylidenecyclopropane-1,1-dicarboxylate.⁴⁸

Scheme 9 Reaction paths involved in the 32CA reaction between AZY 1 and MA 11, including the MPWB1K/6-31G(d) calculated relative Gibbs free energy values in THF.

Fig. 4 illustrates the IRC plots of TS1x, TS1n, TS2x, and TS2n, by which the involvement of the stepwise and one-step mechanisms, respectively, in the endo and exo approach modes of MA 11 toward AZY 1 along both the C3–C4 and C3–C5 regioisomeric attack routes can be clearly understood. These IRC profiles were obtained using 150 points along the forward and reverse directions (total of 301 points), taking a very small step size of 0.020 bohr. The resulting structures at the end points of the IRC profiles of TS1x, TS1n, TS2x, and TS2n convert into CA1x, IN1n, CA2x, and IN2n, respectively, when they are fully optimized.

Fig. 4 MPWB1K/6-31G(d) IRC profiles of TS1x, TS1n, TS2x, and TS2n involved in the 32CA reaction of AZY 1 toward MA 11 in THF, including the resulting structures E, F, G, and H associated with the end points of the corresponding IRC profiles.

The relative Gibbs free energy profile of the 32CA reaction between AZY 1 and MA 11 in THF is depicted in Fig. 5. This profile evidently demonstrates that when the reactants come sufficiently close, the PES of the reaction increases by 9.8 kcal mol⁻¹, reaching the endo TS1n by electrophilic attack of the C4 carbon atom of MA 11 on the C3 carbon atom of AZY 1 to yield the endo intermediate IN1n, whose Gibbs free energy value is 5.2 kcal mol⁻¹ higher than those of the separate reagents.

Fig. 5 MPWB1K/6-31G(d) relative Gibbs free energy profile of the studied 32CA reaction between AZY 1 and MA 11 in THF, in which the differences of 5.2 and 10.1 evidently explain why this reaction experimentally occurs in a C3–C4 regiospecific and an endo-stereospecific manner, respectively, at ambient temperature.

Upon the formation of IN1, an exergonic ring closure step, ΔG = −10.3 kcal mol⁻¹, passing through TS11n with a non-appreciable barrier of 1.3 kcal mol⁻¹, leads to the formation of the endo lithiated cycloadduct CA1n, whose Gibbs free energy is 5.1 kcal mol⁻¹ lower than those of the separate reagents. In excellent agreement with the experimental findings,¹⁸ the energetic results indicate that the studied 32CA reaction in THF takes place via a C3–C4 regiospecific and an endo-stereospecific pathway along which CA1n is generated as the only reachable cycloadduct.

It is worth mentioning that the energetic preference of the C3–C4 regioisomeric pathway is in agreement with the justifications provided in Section 3.1. On the other hand, at first glance, the difference of 5.2 kcal mol⁻¹ between the relative Gibbs free energies of TS2n and TS1n (see Fig. 5) may seem much lower than a value that can explain the regiospecificity experimentally displayed by the investigated 32CA reaction.¹⁸ Considering the activation Gibbs free energies associated with the formation of intermediates IN1n and IN2n as well as cycloadducts CA1n and CA2n, it is evident that along both the C3–C4 and C3–C5 regioisomeric pathways, the formation of intermediates should be taken as the rate-determining step (RDS), i.e. the rates of formation of CA1n and CA2n are controlled by the rates of conversion of the reactants into IN1n and IN2n, respectively. According to basic kinetics, the rates of formation of IN1 and IN2 can be written as k₂[AZY 1][MA 11] and k′₂[AZY 1][MA 11], respectively, where k₂ and k′₂ denote the second-order rate constants. Consequently, the concentration ratio of IN1 to IN2 in the reaction mixture becomes:

Using transition state theory (T.S.T) in solution,⁴⁹ k₂ and k′₂ are related to the relative (activation) Gibbs free energies of TS1n and TS2n, respectively, as follows:

where K_B, h, and C° are the Boltzmann constant, Planck constant, and standard molar concentration (1 mol L⁻¹), respectively. Substitution of eqn (4) into (3) gives

Given ΔG′^≠ = 15.0 and ΔG^≠ = 9.8 kcal mol⁻¹ at 298 K, the ratio becomes 6515, i.e. the percentage becomes almost 100%, evidently indicating that IN1n is the only intermediate present in the reaction medium; hence, the 32CA reaction of AZY 1 toward MA 11 in THF actually takes place via a regiospecific pathway to generate CA1n. Consequently, due to the difference of 5.0 kcal mol⁻¹ in the relative activation Gibbs free energies of the two competitive paths, only one of these may be considered to be formally allowed from the kinetic point of view. However, the relative Gibbs free energy of the generated cycloadduct CA1n, −5.1 kcal mol⁻¹, obviously demonstrates that this species is thermodynamically more stable than the separate reagents. Inclusion of this value in eqn (6) leads to an equilibrium constant (K_eq) larger than 5500 at 298 K, i.e. the generation of CA1n (compared with other CAs) is not only kinetically preferred but is also highly supported thermodynamically. Moreover, the relatively large energy content of 11.6 kcal mol⁻¹ released in the formation of CA1n from TS11n is the thermodynamic driving force of the C3–C4 channel.

The two different mechanistic aspects involved in the studied 32CA reaction can be rationalized by an assessment of the natural atomic charges distributed over the C1, C3, C4, C5, O6, and Li atoms (see Scheme 9 for atom numbering). In other words, we are going to answer this question: “why does the studied 32CA reaction take place via a one-step mechanism along the exo approach but stepwise along the endo one?” In order to address this question, a natural population analysis (NPA) was performed on the optimized structures of AZY 1, MA 11, TS1x, TS1n and IN1n in THF; the results are given in Table 3.

Table 3 MPWB1K/6-31G(d) NPA atomic charges over the considered atoms in AZY 1, MA 11, TS1x, TS1n and IN1n

	NPA atomic charge (e)
	C1	C3	C4	C5	O6	Li
AZY 1	−0.14	−0.39				0.89
MA 11			−0.33	−0.37	−0.64
TS1x	−0.02	−0.34	−0.41	−0.49	−0.73	0.83
TS1n	−0.04	−0.34	−0.39	−0.47	−0.75	0.71
IN1n	0.08	−0.27	−0.48	−0.51	−0.81	0.71

---

From Table 3, one can see that at TS1x, while the C1 and C3 atomic charges shift toward more positive values, a shift toward a more negative atomic charge is experienced by C4 and C5 compared with the corresponding values at AZY 1 and MA 11; this, in accord with the discussion provided in Section 3.1, indicates that the GEDT takes place from the AZY 1 moiety toward the MA 11 moiety. This GEDT flux not only leads to increased polarization of the C4–C5 double bond toward C5, characterized with a more negative charge of −0.49e, but also strongly increases the negative charge on O6 from −0.64e at MA 11 to −0.73e at TS1x, implying a high polarization of the carbonyl double bond toward O6. Despite the mentioned changes, as presented in the last column of Table 3, a slight reduction of 0.06e is observed in the positive charge of Li cation due to the lack of any interaction (electron-density delocalization) between O6 and Li cation along TS1x. The higher difference in the atomic charges of C4 and C5, together with a considerable reduction in the positive charge of Li cation at TS1n compared to TS1x, demonstrate that the electron-density delocalization from LP(O6) toward LP*(Li) leads to a higher polarization of the C4–C5 double bond at TS1n than at TS1x. The donor-LP(O6) → acceptor-LP*(Li) interaction, leading to an E⁽²⁾ value of 13.0 kcal mol⁻¹, is graphically represented in Fig. 6. The interaction depicted in Fig. 6 is one of the determining factors in the predominance of endo-stereospecificity over exo-stereospecificity along the investigated 32CA reaction (see the NCI analysis provided in Section 3.4). Furthermore, analysis of the NPA atomic charges at IN1 obviously shows a positive charge of +0.08e at C1 and a negative charge of −0.51e at C5 (see the last row of Table 3), portraying a highly zwitterionic species.

Fig. 6 MPWB1K/6-31G(d) shape of donor-LP(O6) → acceptor-LP*(Li) interaction leading to an E⁽²⁾ value of 13.0 kcal mol⁻¹ at TS1n located on the PES of the 32CA reaction of AZY 1 toward MA 11 in THF.

In Fig. 7, the zwitterionic nature of IN1n is portrayed by a molecular electrostatic potential (MEP) map in which the negatively charged C5 atom is surrounded by a yellow-red isodensity surface while the positively charged C1 atom is surrounded by a light-green isodensity surface.

Fig. 7 MPWB1K/6-31G(d) MEP map of IN1n portraying a zwitterionic structure in which the C1 and C5 NPA charges are +0.08 and −0.51e, respectively.

The attractive electrostatic force between the C1 and C5 centers allows the ring closure step to easily take place passing through TS11n, with a non-appreciable barrier of 1.3 kcal mol⁻¹ (see earlier). It is highly important to note that the key role of THF, as a relatively polar solvent, in the stabilizing of IN1n should not be ignored. In fact, very interestingly, the IRC profile of TS1n located in a vacuum (gas phase) does not reach a stable intermediate at its end point and stops at a structure that highly resembles cycloadduct CA1n. As nicely portrayed with the aid of the IRC profile of TS1n in a vacuum (see Fig. 8), the shoulder-like region appearing after TS1n indicates that the strong attractive electrostatic force between the oppositely charged C1 and C5 centers cannot effectively be reduced and, therefore, the ring closure tends to spontaneously proceed along a severe downward slope. The key role of solvent polarity in the stabilization of IN1n can simply be described using Coulomb's law, which states “the electrostatic force between two charges is inversely proportional with the dielectric constant of the medium”. Considering this law, it can easily be concluded that the attractive electrostatic force between the C1 and C5 centers at TS1n is reduced by more than seven times when the studied 32CA reaction is transferred from a vacuum (ε = 1.00) into a medium surrounded by THF (ε = 7.42). Consequently, one can conclude that the O6–Li interaction at TS1n leading to a significant charge separation at the resulting zwitterionic IN1n (+0.08e at C1 and −0.51e at C5), which is sufficiently stabilized by the employed solvent, is responsible for the stepwise mechanistic aspect explored along the endo approach mode of the reagents in the considered 32CA reaction.

Fig. 8 MPWB1K/6-31G(d) IRC profile of TS1n in a vacuum.

The MPWB1K/6-31G(d) optimized structures and corresponding Cartesian coordinates of the stationary points involved in the 32CA reaction between AZY 1 and MA 11 in THF (see Scheme 9 for details) together with the unique imaginary frequency of TSs, as well as some selected bond distances, are given in the ESI.†

3.4. The origin of the complete endo stereoselectivity experimentally observed in the 32CA reaction between AZY 1 and MA 11 in THF

The energetic results presented in Section 3.3 perceptibly demonstrate that the 32CA reaction of AZY 1 toward MA 11 in THF initializes with the C4 electrophilic attack of MA 11 on the C3 carbon atom of AZY 1 in an endo-stereospecific manner passing through the endo TS1n (see Fig. 5).

Although stereoselectivity is generally controlled by a combination of stereoelectronically-induced bond orientations in the TS, this factor can be affected by non-bonded (non-covalent) weak interactions between two interacting fragments in the TS.⁵⁰ Recently, it has well been established that reduced density gradient (RDG) index can be used as a very helpful tool to identify and characterize weak interactions (NCIs) of various strengths as chemically intuitive RDG isosurfaces portraying both stabilizing and destabilizing interactions.^50–54 In order to confirm the presence of various stabilizing weak NCIs appearing between two interacting fragments at TS1n and TS1x, NCI analysis of the pro-molecular electron density of these TSs was performed using the NCIPLOT 3.0 software package,^52b and the corresponding NCI gradient isosurfaces are portrayed in Fig. 9. From Fig. 9, it is obvious that there is a much larger green isosurface extended across the region between the two interacting fragments in TS1n than in TS1x, indicating the existence of more stabilizing weak van der Waals NCIs, e.g. the CH⋯π interaction between the CH₃ of the methoxy group in the MA 11 fragment and the pyridine ring of the AZY 1 fragment. Similarly, in accord with the NBO and QTAIM analyses, the gradient isosurfaces that correspond to the weak stabilizing NCIs between Li(I) cation and the N2, N9, and O7 heteroatoms in both TS1x and TS1n can evidently be identified (see violet arrows in Fig. 9). On the other hand, the presence of a unique and distinguishing NCI isosurface arising from the electron-density flux of the O6 lone pairs toward the unfilled orbital of Li cation at TS1n (characterized with a red arrow in Fig. 9), which is absent in TS1x, leads to greater stabilization in TS1n than in TS1x. As a result, the presence of the unique O6–Li NCI together with the large van der Waals NCIs is responsible for the strong predominance of endo TS1n over the exo TS1x, promoting the studied 32CA reaction in an endo stereospecific fashion, in complete agreement with the experimental outcomes.¹⁸

Fig. 9 NCI gradient isosurfaces with an isovalue of 0.3 for TS1x and TS1n involved in the more favourable C3–C4 regioisomeric channel of the 32CA reaction considered in this study. The type and strength of the NCIs are simply characterized using the given colored scale bar.

4. Conclusion

The present work comprises a theoretical investigation on the 32CA reaction between a stabilized N-lithiated azomethine ylide (AZY 1) and the π-deficient methyl acrylate (MA 11) in THF, experimentally studied by Dondas and co-workers, from the energetics, regio- and stereoselectivity, and mechanistic points of view.

The electron densities shared from the lone pairs of N(2), O(7), and N(9) heteroatoms toward the unfilled orbital of Li(8) cation at the ground state geometry of AZY 1 were analyzed and graphically portrayed in terms of NBO and QTAIM approaches to explore how a very reactive azomethine ylide is sufficiently stabilized in the reaction medium.

In excellent agreement with the experimental outcomes, the MPWB1K/6-31G(d) calculated relative Gibbs free energies indicate that this 32CA reaction in THF takes place via a C3–C4 regiospecific and an endo-stereospecific pathway, yielding the corresponding α-amino phosphonate derivative as the only reachable cycloadduct. The validity of the employed 6-31G(d) basis set to obtain the energetic results was examined using the larger basis set of 6-31G+(d). It should be noted that, as discussed in the ESI,† the inclusion of diffuse functions using the 6-31+G(d) basis set not only cannot change the selectivities or improve the kinetics and thermodynamics of the studied 32CA reaction, but also, taking into account that this reaction experimentally takes place at ambient temperature, leads to results that cannot properly justify the experimental outcomes. Consequently, one may conclude that the inclusion of diffuse functions in cycloaddition reactions (at least in the case of 32CA) independent of the involved mechanistic aspects should be avoided, although this claim needs to be proved more decisively by many other investigations.

Mechanistically, the stepwise 32CA reaction of AZY 1 toward MA 11 is initialized when the C3 carbon atom of AZY 1 is electrophilically attacked by the C4 carbon atom of MA 11 along an endo approach mode passing through the endo transition state TS1n (located 9.8 kcal mol⁻¹ over separate reagents) to afford the zwitterionic intermediate IN1n in which the first C3–C4 single bond forms. Unlike the gas phase, the presence of IN1n on the potential energy surface (PES) of the studied reaction in THF is a direct consequence of the solvent polarity by which the significant charge separation at the resulting zwitterionic IN1n (+0.08e at C1 and −0.51e at C5) arising from the electron-density delocalization between the carbonyl oxygen atom of MA 11 and the Li cation of AZY 1 is sufficiently stabilized. Upon the formation of IN1n, the attractive electrostatic forces between the oppositely charged C1 and C5 carbon atoms allow the ring closure step to proceed by the formation of a second C1–C5 single bond via TS11n with a non-appreciable barrier of 1.3 kcal mol⁻¹.

Although the interaction between the C3 carbon atom at AZY 1, the most nucleophilic center, possessing a nucleophilic Parr function of 0.51, and the C4 carbon atom at MA 11, the most electrophilic center, possessing an electrophilic Parr function of 0.58, is able to predict a C3–C4 regioselectivity, the extremely small (even negligible) difference of 0.04 in the calculated nucleophilic Parr functions of the C1 and C3 atoms at AZY 1 cannot explain the C3–C4 regiospecificity experimentally observed in the 32CA reaction under study. Indeed, the large steric hindrances between the CO(OCH₃) substituent in MA 11 and the PO(OCH₃) substituent in AZY 1 along the C3–C5 regioisomeric channel, together with the less important favourable electronic effects along the C3–C4 regioisomeric channel, is responsible for the complete predominance of the latter channel over the former, justifying the C3–C4 regiospecificity.

Finally, non-covalent interactions (NCIs) analysis of the electron density of TS1x and TS1n, which are involved, respectively, in the competitive exo and endo stereoselective approach modes along the favourable C3–C4 regioisomeric channel of the given 32CA reaction provides an explanation of the endo-stereospecificity observed experimentally. In fact, the presence of weak attractive NCIs between the carbonyl oxygen atom of the MA 11 moiety and the Li cation of the AZY 1 moiety at TS1n, which is absent in TS1x, together with some larger attractive van der Waals interactions extending between two interacting fragments at TS1n but not at TS1x, are responsible for the high preference of the endo approach over the exo approach.

References

1. A. E. Speers, G. C. Adam and B. F. Cravatt, J. Am. Chem. Soc., 2003, 125, 4686.
2. A. Krasinski, Z. Radic, R. Manetsch, J. Raushel, P. Taylor, K. B. Sharpless, H. Kolb and C. Kolb, J. Am. Chem. Soc., 2005, 127, 6686.
3. T. S. Seo, X. Bai, H. Ruparel, Z. Li, N. J. Turro and J. Ju, Proc. Natl. Acad. Sci. U. S. A., 2004, 101, 5488.
4. W. Carruthers, Cycloaddition Reactions in Organic Synthesis, Pergamon, Oxford, 1990.
5. L. R. Domingo and S. R. Emamian, Tetrahedron, 2014, 70, 1267.
6. (a) S. Emamian, RSC Adv., 2015, 5, 72959; (b) S. Emamian, T. Lu and F. Moeinpour, RSC Adv., 2015, 5, 62248.
7. D. H. Ess and K. N. Houk, J. Am. Chem. Soc., 2008, 130, 10187.
8. L. R. Domingo, RSC Adv., 2014, 4, 32415.
9. L. R. Domingo, M. Ríos-Gutiérrez, M. Duque-Noreña, E. Chamorro and P. Pérez, Theor. Chem. Acc., 2016, 135, 160.
10. B. Braida, C. Walter, B. Engels and P. C. Hiberty, J. Am. Chem. Soc., 2010, 132, 7631.
11. J. Adrio and J. C. Carretero, Chem. Commun., 2011, 47, 6784.
12. I. Coldham and R. Hufton, Chem. Rev., 2005, 105, 2765.
13. W. Gao, X. Zhang and M. Raghunath, Org. Lett., 2005, 7, 4241.
14. (a) S. Eröksüz, Ö. Doğan and P. Garner, Tetrahedron: Asymmetry, 2010, 21, 2535; (b) A. S. Gothelf, K. V. Gothelf, R. G. Hazell and K. A. Jørgensen, Angew. Chem., 2002, 114, 4410; (c) R. Grigg, Tetrahedron: Asymmetry, 1995, 6, 2475; (d) J. Shi, M. Zhao, Z. Lei and M. Shi, J. Org. Chem., 2008, 73, 305; (e) S. Carbrera, R. Gomez and J. C. Carretero, J. Am. Chem. Soc., 2005, 127, 16394; (f) H. Wu, B. Wang, H. Liu and L. Wang, Tetrahedron, 2011, 67, 1210.
15. J. M. Longmire, B. Wang and Z. Zhang, J. Am. Chem. Soc., 2002, 124, 13400.
16. S. Vivanco, B. Lecea, A. Arrieta, P. Prieto, I. Morao, A. Linden and F. P. Cossío, J. Am. Chem. Soc., 2000, 122, 6078.
17. A. Ponce, I. Alonso, J. Adrio and J. C. Carretero, Chem.–Eur. J., 2016, 22, 4952.
18. H. A. Dondas, Y. Durust, R. Grigg, M. J. Slatera and M. A. B. Sarker, Tetrahedron, 2005, 61, 10667.
19. (a) C. Lee, W. Yang and R. G. Parr, Phys. Rev. B, 1988, 371, 785; (b) A. D. Becke, J. Chem. Phys., 1993, 98, 5648.
20. (a) C. E. Check and T. M. Gilbert, J. Org. Chem., 2005, 70, 9828; (b) G. O. Jones, V. A. Gune and K. Houk, J. Phys. Chem. A, 2006, 110, 1216; (c) G. A. Griffith, I. H. Hillier, A. C. Moralee, J. M. Percy, R. Roig and M. A. Vincent, J. Am. Chem. Soc., 2006, 128, 13130.
21. Y. Zhaoand and D. G. Truhlar, J. Phys. Chem. A, 2004, 108, 6908.
22. R. Jasiński, Comput. Theor. Chem., 2014, 1046, 93.
23. W. J. Hehre, L. Radom, P. V. R. Schleyer and J. A. Pople, Ab initio molecular orbital theory, Wiley, New York, 1986.
24. X. Li and M. J. Frisch, J. Chem. Theory Comput., 2006, 2, 835.
25. K. Fukui, J. Phys. Chem., 1970, 74, 4161.
26. (a) C. González and H. B. Schlegel, J. Phys. Chem., 1990, 94, 5523; (b) C. González and H. B. Schlegel, J. Chem. Phys., 1991, 95, 5853.
27. (a) J. Tomasi and M. Persico, Chem. Rev., 1994, 94, 2027; (b) B. Y. Simkin and I. Sheikhet, Quantum Chemical and Statistical Theory of Solutions-A Computational Approach, Ellis Horwod, London, 1995.
28. (a) E. Cances, B. Mennucci and J. Tomasi, J. Chem. Phys., 1997, 107, 3032; (b) M. Cossi, V. Barone, R. Cammi and J. Tomasi, Chem. Phys. Lett., 1996, 255, 327; (c) V. Barone, M. Cossi and J. Tomasi, J. Comput. Chem., 1998, 19, 404.
29. J. B. Foresman and A. Frisch, Exploring chemistry with electronic structure methods, Gaussian, Inc., Pittsburgh, PA, 2nd edn, 1998.
30. (a) A. E. Reed, R. B. Weinstock and F. Weinhold, J. Chem. Phys., 1985, 83, 735; (b) A. E. Reed, L. A. Curtiss and F. Weinhold, Chem. Rev., 1988, 88, 899.
31. E. R. Johnson, S. Keinan, P. Mori-Sanchez, J. Contreras-Garcia, J. Cohen and A. W. Yang, J. Am. Chem. Soc., 2010, 132, 6498.
32. (a) J. R. Lane, J. Contreras- Garćıa, J. P. Piquemal, B. J. Miller and H. G. Kjaergaard, J. Chem. Theory Comput., 2013, 9, 3263; (b) J. Contreras-Garcia, E. R. Johnson, S. Keinan, R. Chaudret, J. P. Piquemal, D. N. Beratan and W. Yang, J. Chem. Theory Comput., 2011, 7, 625.
33. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Son-nenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hase-gawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery Jr, J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Starov-erov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Ad-amo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Mar-tin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, O. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski and D. J. Fox, Gaussian 09 (Revision A.02-SMP), Gaussian Inc., Wallingford, CT, 2009.
34. R. G. Parr, L. V. Szentpaly and S. B. Liu, J. Am. Chem. Soc., 1999, 121, 1922.
35. (a) R. G. Parr and R. G. Pearson, J. Am. Chem. Soc., 1983, 105, 7512; (b) R. G. Parr and W. Yang, Density Functional Theory of Atoms and Molecules, Oxford University Press, New York, 1989.
36. L. R. Domingo, E. Chamorro and P. Pérez, J. Org. Chem., 2008, 73, 4615.
37. W. Kohn and L. Sham, J. Phys. Rev., 1965, 140, 1133.
38. L. R. Domingo, P. Pérez and J. A. Sáez, RSC Adv., 2013, 3, 1486.
39. P. Geerlings, F. De Proft and W. Langenaeker, Chem. Rev., 2003, 103, 1793.
40. L. R. Domingo, RSC Adv., 2014, 4, 32415.
41. L. R. Domingo, M. J. Aurell, P. Pérez and R. Contreras, Tetrahedron, 2002, 58, 4417.
42. P. Jarmillo, L. R. Domingo, E. Chamorro and P. Pérez, J. Mol. Struct.: THEOCHEM, 2008, 865, 68.
43. F. Weinhold and C. R. Landis, Valency and Bonding: A Natural Bond Orbital Donor–Acceptor Perspective, Cambridge University Press, New York, 2005.
44. S. Bagchi, D. Mandal, D. Ghosh and A. K. Das, J. Phys. Chem. A, 2013, 117, 1601.
45. F. Biegler-König and J. Schönbohm, AIM2000 Version 2.0, 2002.
46. (a) R. Jasiński, RSC Adv., 2015, 5, 101045; (b) A. F. Khlebnikov, A. S. Koneva, A. A. Virtseva, D. Yufit, G. Mlostoń and H. Heimgartner, Helv. Chim. Acta, 2014, 7, 453.
47. (a) R. Jasiński, Tetrahedron, 2013, 69, 927; (b) R. Jasiński, Tetrahedron Lett., 2015, 56, 532.
48. A. K. Nacereddine, C. Sobhi, A. Djerourou, M. Ríos-Gutiérrez and L. R. Domingo, RSC Adv., 2015, 5, 99299.
49. (a) S. Glasstone, K. J. Laidler and H. Eyring, The theory of rate processes, McGraw-Hill, New York, 1941; (b) K. J. Laidler, Theories of chemical reaction rates, McGraw-Hill, New York, 1941.
50. A. Armstrong, R. A. Boto, P. Dingwall, J. Contreras-Garćıa, M. J. Harvey, N. J. Mason and H. S. Rzepa, Chem. Sci., 2014, 4, 2057.
51. E. R. Johnson, S. Keinan, P. Mori-Sanchez, J. Contreras-Garcia, J. Cohen and A. W. Yang, J. Am. Chem. Soc., 2010, 132, 6498.
52. (a) J. R. Lane, J. Contreras- Garćıa, J. P. Piquemal, B. J. Miller and H. G. Kjaergaard, J. Chem. Theory Comput., 2013, 9, 3263; (b) J. Contreras-Garcia, E. R. Johnson, S. Keinan, R. Chaudret, J. P. Piquemal, D. N. Beratan and W. Yang, J. Chem. Theory Comput., 2011, 7, 625.
53. L. R. Domingo, M. Ríos-Gutiérrez and S. Emamian, RSC Adv., 2016, 6, 17064.
54. A. K. Nacereddine, C. Sobhi, A. Djerourou, M. Ríos-Gutiérrez and L. R. Domingo, RSC Adv., 2015, 5, 99299.

---

Footnote

† Electronic supplementary information (ESI) available: Absolute Gibbs free energy values, G in a.u., obtained via MPWB1K/6-31G(d) frequency calculation, total electronic energy, E^ele_{MPWB1K/6-31G(d)} in a.u., attained through MPWB1K/6-31+G(d) single-point energy calculation over the MPWB1K/6-31G(d) optimized structures, and absolute Gibbs free energy values, G(corrected) in a.u., calculated using eqn (1), of the stationary points involved in the studied 32CA reaction in THF. The relative corrected Gibbs free energy profile, ΔG(corrected), obtained via eqn (1), of the studied 32CA reaction in THF followed by a detailed discussion about how the inclusion of diffuse functions leads to results that cannot properly justify the experimental findings. MPWB1K/6-31G(d) computed total electronic energies (E, in a.u.) and Cartesian coordinates of the optimized structures involved in the 32CA reactions of AZY 1 toward MA 11 in THF (see Scheme 9 in the text for details). While the unique imaginary frequency of the TSs is also given in cm⁻¹, the bond distances in the given optimized structures are in Å. See DOI: 10.1039/c6ra13913b

---