Exploring the Relationship between Reactivity and Electronic Structure in Isorhodanine Derivatives Using Computer Simulations

The electronic structure and reactivity of 22 isorhodanine (IsRd) derivatives in the Diels–Alder reaction with dimethyl maleate (DMm) were investigated under two different environments (gas phase and continuous solvent CH3COOH), using free Gibbs activation energy, free Gibbs reaction energy, and frontier molecular orbitals to analyze their reactivity. The results revealed both inverse electronic demand (IED) and normal electronic demand (NED) characteristics in the Diels–Alder reaction and also provided insights into the aromaticity of the IsRd ring by employing HOMA values. Additionally, the electronic structure of the IsRd core was analyzed through topological examination of the electron density and electron localization function (ELF). Specifically, the study demonstrated that ELF was able to successfully capture chemical reactivity, highlighting the potential of this method to provide valuable insights into the electronic structure and reactivity of molecules.

The Diels-Alder (DA) reaction is classified as a [4+2] cycloaddition due to the interaction of the 4π electron system (the diene structure) with the 2π electron system (the dienophile structure). According to molecular orbital theory, the highest occupied molecular orbital (HOMO) of the diene overlaps with the lowest unoccupied molecular orbital (LUMO) of the dienophile. A key requirement for this reaction is the presence of a strong dienophile with electron-acceptor properties, which decreases the energy difference between the diene's HOMO and the dienophile's LUMO. These considerations define the normal electron-demand DA reaction. A conceptually different type of this reaction is the inverse electron-demand DA reaction, which occurs between an electron-rich dienophile and an electron-poor diene. The energy gap between the diene's LUMO and the dienophile's HOMO is smaller than between the diene's HOMO and the dienophile's LUMO, making the interaction between these two orbitals the most energetically significant stabilizing orbital interaction. According to Conceptual Density Theory, the electronic properties and reactivity of a molecule can be attributed to the energy difference between HOMO and LUMO orbitals, which determine the molecule's ability to donate or accept electrons. The energy gap between HOMO and LUMO orbitals affects ionization energies and the corresponding changes in charge densities, which are essential in deriving various reactivity parameters, such as electronegativity, hardness, and softness [19,20]. Lately, Domingo and others have applied the ideas on DA and various other reactions, while Shaik has developed a detailed model for organic reactivity that utilizes a valence bond approach and VB diagrams to demonstrate the evolution of different valence bond structures throughout reaction pathways [21][22][23]. In recent years, the Distortion/Interaction Activation Strain model has emerged as a novel approach to understanding reactivity. This model explains the origin of activation barriers in chemical reactions based on the distortion of the reactants and the interaction between the distorted species [24].
The mechanism of DA reaction was investigated using chemical catastrophe theory, which provides a complementary approach to understanding the reaction mechanism based on quantum chemical calculations. By mapping out the reaction pathway and analyzing the electronic and geometric changes that occur during the reaction, catastrophe theory has provided valuable insights into the complex reaction mechanism [25]. Interestingly, to date, no correlation was found between reactivity and the populations of chemical bonds or lone pair electrons calculated through the electron localization function (ELF). Additionally, it should be noted that the reactivity of dienes in DA reaction strongly depends on the electron density on the multiple bonds. The main aim of this work is to present successfully captured chemical reactivity using ELF map in DA reaction, where charge-transfer is the leading term [26]. This demonstrates that ELF method can provide valuable insight into the electronic structure and reactivity of molecules.
The chemical reactivity of the 5-arylidene derivative of 4-sulfanylidene-1,3-thiazolidine-2-one, commonly known as isorhodanine (IsRd), and its isomer, 2-sulfanylidene-1,3thiazolidine-4-one, commonly named rhodanine, has been previously investigated using quantum chemical methods, in both gas-phase and solvent (acetic acid) models, revealing low energy barrier for IsRd in DA reaction [27]. The current study focuses on analyzing the reactivity between IsRd derivatives with dimethyl maleate (DMm). In this paper, the analysis of geometrical structures and energetic parameters, such as the free Gibbs activation energy and the free Gibbs reaction energy, is followed by a topological analysis of the electron density using AIM analysis proposed by Bader [28] and ELF proposed by Silvi and Savin [29,30]. Finally, the results obtained from these analyses are discussed and summarized.

Materials and Methods
The geometry optimizations were performed using the density functional theory (DFT) within Gaussian 16 (G16) programme (version C.01) [31] and Minnesota group functional, M06-2X [32], included in G16. The Pople basis set, 6-311+G(d,p) [33,34], with diffuse and polarisation functions were used in these calculations. The study considered two distinct environments: the gas phase and acetic acid used as a solvent. To simulate the acetic acid effects, Polarizable Continuum Model (PCM, ε = 6.2528) was used, which was previously validated for organic molecules [35][36][37][38][39][40]. The free Gibbs activation energy (∆G a ) was calculated as a difference between sum of electronic and thermal free energies between TS and optimised structures of the reagents. The free Gibbs reaction energy, ∆G r , has been calculated as difference of sum of electronic and thermal free energies between products and reagents. The energy calculations were performed using gas-phase and PCM models at a temperature of 298.15 K and a pressure of 1 atm.
The electron density, ρ(r), and ELF, η(r), were analyzed to determine chemical bonds using topological analysis. A grid points with a step size of 0.05 bohr was used and calculations were carried out with Multiwfn [41] and TopMod09 [42] packages. The Multiwfn program was also used to calculate the harmonic oscillator model of aromaticity (HOMA) values [43,44] on DFT optimized molecules. Finally, graphical representations of optimized structures and isosurfaces were generated using the Chimera program [45].

Results and Discussion
The DA reaction is a powerful and well-known method for forming cyclic compounds through the [4+2] cycloaddition of a conjugated diene and a dienophile. The mechanism of DA reaction is characterized by the breaking of three π-bonds and the formation of two single σ-bonds and one π-bond (see Figure 1). The reaction process proceeds through a single transition state, which starts with the formation of the pre-reaction complex. The optimized −H, −NO 2 , and −NH 2 derivatives of DA reaction are presented in Figure 2. These derivatives, also known as substituents, can have a significant impact on the reactivity and selectivity of DA reaction by modifying the electronic and steric properties of the reacting molecules.

Geometrical Structures and Energetics
The lengths of the C3-C4 and C1-S1 bonds for the optimized structures of the TS are shown in Table 1. The C1-S1 and C3-C4 bond lengths, r opt (C, S) and r opt (C, C), in the gas-phase study varies between 1.951 (−O − ) to 2.545 (−NO 2 ) Å and 2.205 (−NH 2 ) to 2.293 (−SO 3 H) Å, respectively. The values of r opt (C, S) are significantly larger than r opt (C, C). However, for the −O − , −S − , −N(CH 3 ) 2 , −N(CH 2 CH 3 ) 2 , and −N(Pr) 2 calculated r opt (C, S) is smaller than r opt (C, C). The application of the PCM reveals similarity of r opt (C, S) and r opt (C, C) values to those obtained in the gas-phase and range 1.974-2.543 Å and 2.181-2.346 Å, respectively. Although obtained C1-S1 and C3-C4 bond lengths differ for −O − , −S − , −N(CH 3 ) 2 , −N(CH 2 CH 3 ) 2 , and −N(Pr) 2 , the r opt (C, S) and r opt (C, C) values are in the typical range of carbon-sulphur and carbon-carbon interactions [46] and the selected computational method, DFT(M06-2X)/6-311+G(d,p), can be used to study the reactivity of IsRd derivatives with DMm. Table 1. The values of the selected bond lengths, r opt (Å), between the interacting atoms in IsRd · · · DMm TS derivatives, the free Gibbs activation energy, ∆G a (kcal/mol), and the free Gibbs reaction energy, ∆G r (kcal/mol), calculated for DA reaction with DMm.

Gas Phase
Solvent (Acetic Acid) The values of the interaction energy (E int ) are presented in Table 2, corrected using the counterpoise correction method [47]. The counterpoise correction method is a computational technique that is used to correct the binding energy of a complex for the basis set superposition error (BSSE), which is an artifact that arises when the monomers of a complex are treated with different basis sets. This method is based on the concept of ghost atoms, which are fictitious atoms that are introduced to cancel out the BSSE by making the monomers of the complex equivalent to the complex. Depending on the derivative, the value of E CP int + ∆ZPVE ranges between −14.82 (−SO 3 H) and −6.41 (−OH) kcal/mol, indicating that the interaction energy does not show significant deviation across the studied compounds. The values of the free Gibbs activation energy, ∆G a , and the free Gibbs reaction energy, ∆G r , in the gas phase are presented in Table 1. The value of ∆G a ranges from 18.47 (−SO 2 CF 3 ) to 28.67 (−NH 2 ) kcal/mol, depending on the considered derivative. The value of ∆G r for meta-directing deactivators (electron withdrawing groups, EWGs) ranges from −22.08 (−CN) to −34.09 (−NO 2 ) kcal/mol and is smaller than that obtained for −H, −18.67 kcal/mol. The free Gibbs reaction energy for −Br and −Cl groups is in the range of EWGs. The value of ∆G r for ortho-and para-directing activators, EDGs, ranges between −25.07 (−Ph) and 23.39 (−O − ) kcal/mol and is notably higher than for EWGs. In all cases, except for −O − and −S − , the product of the reaction is stabilized due to a larger free activation energy on the side of the product, and the probability of the reaction occurring in both directions is small. To conclude this part, the cycloaddition of IsRd with DMm is an exothermic process, and the formation of the product is favored. In the −O − and −S − cases, the cycloaddition of IsRd with DMm is an endothermic process, and the forwards reaction is not favored.
In the presence of the solvent (acetic acid, CH 3 COOH), the value of ∆G a ranges be- The reactivity and selectivity of a molecule can be well understood by the examination of frontier orbitals [48]. Specifically, HOMO and LUMO can provide insight into the reactivity of IsRd derivatives by determining the energy differences between them for optimized geometries, which is vital in understanding the mechanism of DA reaction. The 3D plots of HOMO and LUMO orbitals for −H are shown in Figure 3, and the energy values of the MOs are compared and stored in   The character of DA reaction can be understood by looking at the energy differences between HOMO and LUMO orbitals of diene and dienophile. When the energy gap between LUMO of the dienophile and HOMO of the diene is greater than the energy gap between LUMO of the diene and HOMO of the dienophile, it suggests that the reaction has an inverse electronic demand (IED) character. On the other hand, if the energy gap between LUMO of the dienophile and HOMO of the diene is smaller, it suggests that the reaction has a normal electronic demand (NED) character. The energy gaps for EWGs between DMm LUMO -IsRd HOMO are higher than IsRd LUMO -DMm HOMO , and suggest the inverse electronic demand (IED) character of DA reaction. Opposite conclusions were obtained for compounds with EDGs due to smaller DMm LUMO -IsRd HOMO to IsRd LUMO -DMm HOMO values. This result indicates the normal electronic demand (NED) character.

The Electronic Properties of Isorhodanine Derivatives
One of the computational methods to investigate the electronic structure of chemical bonds is the topological analysis of electron density [28] and ELF [30], grouped under the umbrella of the quantum chemical topology (QCT) approach [49]. The topological analysis of ELF divides the molecular space into regions that correspond to atomic cores, lone pairs, and chemical bonds, thus having a clear chemical meaning. These spaces are described by basins and critical points and have a one-to-one representation of the expected chemical objects. The critical points are characterized by the number of positive eigenvalues of the Hessian. The local maxima called attractors are critical points where three negative eigenvalues of the Hessian occur. The basins have at least one attractor and may be of two kinds: core basins, C(A), which correspond to the electron density around nuclei; and valence basins, V(A, B, . . . ), which describe valence electrons. The valence basins are characterized by their synaptic order, the number of core basins with which a common boundary is shared. The concept of synapticity has been described by Silvi et al., in reference [50]. Monosynaptic basins are associated with lone pairs, and disynaptic ones with two-center covalent bonds. Polysynaptic basins are characteristic of multicenter bonds.

Electron Density
Electron Localization Function ρ ( The delocalization index is a quantitative measure of the number of electron pairs that are delocalized between two atomic basins. It is often considered to be equivalent to the topological bond order [51]. The median number of electron pairs exchanged between the atomic basin ranges between 0.728 (N-H1) and 1.723 (C1-S1). Thus, for N-H1, C4-H2, N-C2, N-C1, C1-C3, C2-S2, and C3-S2, the DI suggests the presence of a single bond. The highest DI median values were obtained for C1-S1 (DI = 1.723) and C3-C4 (DI = 1.615) bonds. According to the DI interpretation, they may be considered as a double type. The classification of C2-O (DI = 1.340) interaction is ambiguous due to the significantly higher DI than in a single type bonding, 0.728-1.116, but also lower than for a double type, 1.615-1.723. In the second step, topological analysis has been carried out for ELF. The electronic structure of the heterocyclic ring in all studied derivatives was found to be characterized by five core and five valence attractors. Figure [30], sharedelectron interactions, such as covalent, dative, and metallic bonds occur when there is at least one bond attractor between the core attractors of the atoms involved in the bond. Thus, the localization of those basins proof that the N-C2, N-C1, C1-C3, C2-S2, and C3-S2 interactions are covalent, with shared electron density. The bonding attractors are localized in the regions with high electron localization. The median population equals 2.08e, 2.11e, 2.28e, 1.93e, and 1.77e.
The formal Lewis structure of IsRd presents double character of C1-S1, C2-O, and C3-C4 bonds, which are described by V(C1,S1), V(C2,O), and V i=1,2 (C3, C4) basins. In the EDGs case, such as −N(CH 3 ) 2 , −N(CH 2 CH 3 ) 2 , −N(Pr) 2 , −NHCH 3 , −NH 2 , −S − , and −O − the electron density between C3 and C4 atoms is characterized by a single V(C3,C4) basin. The basin population equals 2.55e, 2.46e, and 3.48e, respectively. The topological bond order for C1-S1, C2-O, and C3-C4 equals 1.28, 1.23, and 1.74. For C1-S1 and C2-O bond orders, they are closer to 1, which indicates the single nature of those interactions. The C-O and C-S bonds in H 2 C = X (X = O, S, Se, Te) have been already studied by Berski et al. [52]. The delocalization of the electron density in H 2 C = X is dominated by an exchange of electrons between lone pairs of chalcogens for X = (O, S) and delocalization between the lone pairs and core basin for X = (Se, Te). For −H the standard deviations of theN[V(C1,S1)] and N[V(C2,O)] are the same and yields 1.16. A high value of the σ confirms the mechanism of fluctuating electron density in the C1-S1 and C2-O bonds. Such results support that the nature of carbon-oxygen and carbon-sulfur interaction should be described by two resonance structures, C + O − , C − O + and C + S − , C − S + , rather than double bonding pairs, C=O and S=O.
The hetero-DA reaction is a variant of DA reaction that involves the use of heteroatoms, such as carbonyls and imines, in π-systems. Heteroatoms, defined as atoms other than carbon and hydrogen, play a crucial role in this reaction by introducing new chemical and electronic properties to the reacting molecules. In the hetero-DA reaction, these atoms are utilized to modify the reactivity and selectivity of DA process. In IsRd π-systems are represented by bonding, disynaptic V(C1,S1) and V i=1,2 (C3, C4) basins. The populations of these basins range from 2.26 to 2.68e and 2.59 to 3.63e, respectively, depending on the studied derivative. It is noteworthy that the population of the V i=1,2 (C3, C4) basin is affected by the character of nearby chemical substituents, with this effect being more pronounced than for the V(C1,S1) basin. The reactivity of dienes is found to significantly increase when EDG is attached. To investigate this phenomenon, linear regression was utilized to show the relationship between ∆G r andN∑[V(C3,C4),V(C1,S1)] using data obtained at the DFT(M06-2X) level. The results of this analysis are presented in Figure 5, which shows that diene reactivity correlates with electron density on both double bonds, with a coefficient of determination, R 2 , equal to 0.93. Unfortunately, a similar correlation based on the free Gibbs activation barriers andN∑[V(C3,C4),V(C1,S1)] cannot be presented. The highest value ofN∑[V(C3,C4),V(C1,S1)] was observed for −NO 2 (6.26e: −34.09 kcal/mol) and the lowest for −O − (4.85e: 23.39 kcal/mol). To conclude, topological analysis of ELF reveals a strong impact of electron density on both multiple bonds on diene reactivity in DA reaction. Finally, the aromaticity of the IsRd derivatives is analyzed using HOMA values [43,44]. This method compares bond lengths between chemical bonds in the studied rings and the idealized bond length in an aromatic ring in benzene. To gain insight into the influence of chemical groups on electron density in the heterocycle ring, we calculated the HOMA index for all studied derivatives at the DFT(M06-2X)/6-311+G(d,p) computational level. All calculations were performed for isolated molecules in the gas phase. The calculated HOMA values are presented in Table 4. The aromaticity of IsRd with −H is intermediate between meta-directing deactivators and ortho-and para-directing activators, and it equals −0.009. For molecules with strong aromatic rings (benzene), the HOMA index should be close to 1.00. Thus, the derivatives with withdrawing groups can be considered slightly antiaromatic, as all HOMA values are clearly smaller than 1.00. When donating groups are attached, HOMA values increase, and the highest values were observed for −O − , −S − , −NH 2 , −NHCH 3 , −N(Pr) 2 , −N(CH 2 CH 3 ) 2 , −N(CH 3 ) 2 , −NHCOCH 3 . For these substituents, the HOMA index ranges between 0.306 and 0.467, indicating a much higher aromaticity character of the heterocycle ring. Interestingly, −Br and −Cl do not behave as typical withdrawing groups, where the HOMA index is negative, but as donating groups with positive values, 0.049 and 0.066, respectively.

Conclusions
The reactivity and electronic structure of 22 IsRd derivatives were examined using thermodynamics, frontier molecular orbitals, topological analysis of ELF and electron density, and HOMA. This research enhances our understanding of IsRd reactivity. The main findings can be summarized as follows: 1.
The withdrawing and donating substituents produce similar ∆G a values, and one may assume similar reactivity. However, different conclusions were reached by analyzing ∆G r values. The free Gibbs reaction energy depends heavily on the chemical group next to the C=C bond. Smaller ∆G r values were obtained for reactions with EWGs, indicating that the cycloaddition of IsRd with DMm is more favored for derivatives with withdrawing groups. Only for the −O − derivative in the gas phase and PCM model, ∆G r is positive, making the cycloaddition an endothermic process. Comparison of ∆G a and ∆G r in acetic acid and in the gas phase shows that acetic acid has a minimal impact on the energetics of the reaction between IsRd and DMm.

2.
The analysis of frontier DFT orbitals performed on molecules with EWGs shows a decrease in the electronic energy of HOMO and LUMO orbitals, indicating that a NED character of DA reaction is expected. For EDGs, the energy of frontier MOs increases, resulting in an IED type of cycloaddition reaction.

3.
Topological analysis of ρ(r) shows positive values of the Laplacian for BCPs only for the C1-S1 and C2-O bonds, indicating a depletion of electron density around BCP. The Laplacian values for other bonds are negative, indicating a concentration of electron density around BCP. Thus, C1-S1 and C2-O bonds exhibit a different nature than observed for other covalent bonds where both atoms share electron density.

4.
Topological analysis of ELF partly confirms the results based on electron density. The topological bond orders for C1-S1 and C2-O are 1.28 and 1.23, respectively, thus one may consider these interactions to be of a single nature. Analysis of oe suggests that the C + O − , C − O + and C + S − , C − S + representations for the carbon-oxygen and carbon-sulphur interactions, instead of the classical C=O and S=O formulas, provide a better description of the bond nature. 5.
The most intriguing discovery was obtained by analyzing the relationship between ∆G r and the basin populations for the C=S and C=C bonds in the IsRd molecule. The analysis of ELF was performed on isolated molecules. According to the proposed DA reaction mechanism, the electron density from these bonds is distributed to the regions of the two new C-S and C-C bonds that are formed between the IsRd and DMm molecules. The regression analysis applied to this relationship shows that the values of the sum of the populations for the V(C1,S1) and V i=1, 2 (C3, C4) basins,N∑[V(C3, C4), V(C1, S1)], correlates with the value of the free Gibbs reaction energy. Large negative values of ∆G r correspond to large values of the sum of basin populations. The favorability of the studied DA reaction, which is associated with a rearrangement of chemical bonds and electron density, is high when the C=S and C=C bonds are better saturated with electron density. 6.
The HOMA values indicate low aromaticity of the IsRd ring.