Chemical reactivity of quinclorac employing the HSAB local principle - Fukui function

Abstract In the present work we have calculated several DFT reactivity descriptors for quinclorac at the B3LYP/6- 311++G(2d,2p) and MP2/6-311++G(2d,2p) levels of theory in order to analyze its reactivity. Reactivity descriptors such as ionization energy, molecular hardness, electrophilicity, condensed Fukui function and total energies were determined to predict the reactivity of quinclorac. The influence of the solvent was taken into account employing the PCM model. The results indicate that the solvation modifies the values of quinclorac reactivity descriptors. The Fukui function values predict that an electrophilic attack on quinclorac might cause a dechlorination, while a nucleophilic attack might lead to a decarboxylation and a free radical attack would cause a hydrogen substitution on the quinoline ring. Quinclorac in deprotonated form would be susceptible to decarboxylation through an electrophilic attack while nucleophilic and free radical attacks would cause an attack on the hydrogens of the ring. Graphical Abstract


Introduction
Auxinic herbicides are widely used to control weeds in crops because they were the first selective herbicides developed [1,2]. Quinclorac (3,7-Dichloro-8quinolinecarboxylic acid, see Fig. 1) is probably one of the auxinic herbicides most used to increase rice production because it can be applied throughout the growing season [3]. Quinclorac has been used effectively in controlling dicotyledon and monocotyledon weeds such as grass species of Echinochloa, Digitaria, and Setaria, with excellent crop safety [4][5][6]. In recent years, there have been major concerns about quinclorac herbicide residues and associated food safety issues, their negative impacts on the environment and the increasing occurrence of herbicide resistance in weed populations [7][8][9][10][11]. Marchezan et al. found quinclorac in the water of rivers in considerable concentrations, sufficient to harm the local benthic community [8]. Additionally, quinclorac may be expected to have greater effects on algae and aquatic plants, and may be toxic to birds, bees or earthworms [9]. Although quinclorac belongs to an important class of herbicides, its precise mode of action is still controversial [2,6]. Originally, its mode of action was proposed to be auxin-like, largely based on its morphological response [10]. However, Koo et al reported that quinclorac did not exhibit auxin-like activity in barnyard grass based on a mesocotyl elongation assay [11][12][13]. In other studies, it was found that quinclorac induces oxidative damage in grass, suggesting a relationship between reactive oxygen species and herbicidal activity [14,15]. Also, an enhancement of quinclorac activity under illumination in maize seedlings has been reported [16][17][18]. On the other hand, some studies have revealed that quinclorac is stable to aqueous hydrolysis [19]. It degrades to 3-Chloro-8quinoleincarboxylic acid resulting from a dechlorination reaction [19]. Also, it has been reported that quinclorac is stable to hydrolytic degradation between pH 3-9, but it may be degraded in natural water (pH 8) when it is exposed to UV light (254 nm) or natural sunlight. Multiple photodegradation products are generated, since the main product is identified as 3,7 dichloroquinoline, which results from the decarboxylation of the parent molecule [19]. Recently, two possible quinclorac degradation pathways were reported. In the first pathway the degradation process begins with an oxidative dechlorination at carbon 3 (see Fig. 1 [20]). The second pathway begins with an electrophilic aromatic substitution at carbon 2 prior to the opening of the pyridine ring [20]. This indicates that more progress in the elucidation of the molecular reactivity of quinclorac is required to identify and understand its specific reactivity. To our best knowledge, a quantum chemical study of quinclorac to evaluate its global and local reactivity descriptors is still missing. Thus, in the present work we have analyzed the molecular reactivity of quinclorac through global reactivity descriptors and the Fukui Function derived from Density Functional Theory (DFT) [21]. We consider that this kind of study will contribute to a better understanding of the chemical behavior of this important herbicide both in the gas and solution phases.

Theoretical details
From DFT it is possible to define universal concepts and descriptors of molecular stability and reactivity such as the electronic chemical potential (μ), absolute hardness (η), and global electrophilicity index (ω) [21][22][23][24][25][26][27]. The electronic chemical potential μ (which is equal to the negative of the electronegativity of atoms and molecules) was defined by Parr and Pearson [28] as (1) where I is the vertical ionization energy and A stands for the vertical electron affinity. Absolute hardness can be shown to be [27-31]: (2) Furthermore, the global electrophilicity index ω was introduced by Parr [25,32,33] and may be calculated using the electronic chemical potential μ and the absolute hardness η: According to this definition ω measures the propensity of a species to accept electrons. Thus, a good nucleophile has low values of ω while a good electrophile is characterized by high values. Also, the hard and soft acids and bases principle (HSAB) has been useful to predict the reactivity of chemical systems [30,31,33,[34][35][36][37][38]. The HSAB principle has been used in a local sense in terms of DFT local properties such as the Fukui function [39,40]. Gázquez where is the electronic density, N is the number of electrons and ν is the external potential exerted by the nuclei. This parameter is related to the frontier orbitals within the frozen core approximation [34]. This approximation considers that when there is a variation in the number of electrons, only the respective frontier orbital is affected, and when N diminishes to N -dN the Fukui function for an electrophilic attack turns out to be: is the electronic density of the highest occupied molecular orbital (HOMO). When N increases to N + dN the Fukui function for a nucleophilic attack becomes: , where is the electronic density of the lowest unoccupied molecular orbital (LUMO). The Fukui function is a local reactivity descriptor that indicates the preferred regions where a chemical species will change its density when the number of electrons is modified [45][46][47][48][49][50]. Different procedures have been reported in the literature to obtain the change in the electronic density in each atom [45][46][47][48][49][50]. Bultinck et al. proposed that the density on an atom, A, may be evaluated as:  where ( ) is a weight function dependent on the number of electrons contained in the molecule. Thus, this weight function will be different whether the molecule has N, N -dN or N + dN electrons even if the geometry is the same. Also, it is important to mention that the weight functions for all atoms in a molecule, M, always sum to unity, and that usually a positive definite weight function is used [45][46][47][48][49][50].
(8) Also, Bultinck et al. proposed that the condensed Fukui function can be computed using either the fragment of molecular response (FMR) approach or the response of molecular fragment (RMF) approach [50]. In the case of the Mulliken approach, the weight factor derived from Bultinck et al is not dependent on the number of electrons and the two approaches lead to the same expressions to evaluate the condensed Fukui function. Here it is convenient to remember that Mulliken charges and natural charges are both based on orbital occupancies and similar expressions may be obtained to evaluate the condensed Fukui function. Thus, under these conditions, it is possible to define the corresponding condensed or atomic Fukui function on the atom A for an electrophilic attack as [50]: is the population on the atom A in the molecule with N electrons, and it can be evaluated through , and ( ) is the atomic charge on the atom A for the molecule M with N electrons. Note that Eq. 9 is the atom condensed Fukui function introduced by Yang and Mortier [45]. The nucleophilic attack is given by: For a free radical attack the expression is given by: Here it is important to mention that independently of the approximations used to calculate the Fukui function, all of them comply with the exact closure equation [51]: which is important for the use of the Fukui function as an intramolecular reactivity index.

Methodological procedure
A starting geometry was generated using the PM6 method [52] implemented in Mopac2012 [53]. The optimal conformation was subjected to full geometry reoptimization in the gas phase employing the hybrid B3LYP functional [54-56] and the 6-311++G(2d,2p) basis set [57,58]. The optimized herbicide molecules in the gas phase were further reoptimized at the B3LYP/6-311++G(2d,2p) level employing the PCM solvation model [59,60]. The vibrational frequencies were computed to make sure that the stationary points were minima in the potential energy surface. Finally, the optimized B3LYP/6-311++G(2d,2p) geometries were reoptimized employing the second order Moller Plesset theory (MP2) [61] with the 6-311++G(2d,2p) basis set. The energy values and atomic charges of the ionic states of quinclorac, cation and anion, were calculated at the B3LYP/6-311++G(2d,2p) and MP2/6-311++G(2d,2p) levels, in gas and aqueous phases, by using the geometry of the neutral system. The final atomic charges for quinclorac in gas and aqueous phases were obtained in the framework of MP2 theory. In all calculations of these atomic charges, the DENSITY=MP2 option was used.

Results and discussion
The quinclorac molecule was first optimized at the B3LYP/6-311++G(2d,2p) level in the gas phase. Note that this molecule can exhibit two possible conformers as shown in Figs. 2a and 2b. The total energy calculated for the structure depicted in Fig. 2a, conformer I, was -1509.92033023 hartrees, while the energy for conformer II (Fig. 2b) was -1509.92033024 hartrees, suggesting that the two structures are equivalent.
In order to analyze the effect of solvation on the electronic properties of quinclorac, the optimized structures in the gas phase were used as starting points to reoptimize at the B3LYP/6-311++G(2d,2p) level without any symmetry constraints employing the PCM solvation model [59,60]. In this model, the solvent is treated as an unstructured polarizable continuum, characterized only by its dielectric constant which is 78.5 for water at 25°C. For both conformers I and II, the total energy was -1509.94479921. The energy difference between the structures in the gas phase and the aqueous phase is 15.35 kcal mol -1 which corresponds to the calculated solvation energy of quinclorac. It is important to mention that no significant differences in either distances or angles were obtained when the solvent effect was considered for conformers I and II. Moreover, the bond distances were very similar to those reported for the experimental structure [65]. All frecuency values calculated at the B3LYP/6-311++G(2d,2p) level in the gas and aqueous phases were positive. All vibrations are in good agreement with values reported in the literature [66]. A summary of the main bands is depicted in Fig. 3. The optimized structures at the B3LYP/6-311++G(2d,2p) level were further reoptimized at the MP2/6-311++G(2d,2p) level in gas and aqueous phases. For both conformers I and II, the total energy in the gas phase is the same. Also, an identical result was obtained in the aqueous phase. In Table 1, we report the values of the electronic energies calculated for quinclorac with charges +1, 0 and -1 at the MP2/6-311++G(2d,2p) level of theory. From the values obtained at this level, it can be observed that there is not any significant energy difference between conformers I and II. Also, note that the B3LYP/6-311++G(2d,2p) calculations reproduce the experimental geometry well and the same energy trend is obtained at the MP2 level. This result compares favorably with the values reported by Guo et al. [65], where quinclorac appears as a twin structure of both conformers in a 1:1 ratio in its X-ray crystal structure, which could be expected for two isoenergetic conformers. However, here it is important to mention that one conformer may be converted into another via an inversion symmetry operation. Thus, under the Born-Oppenheimer approximation, these conformers have an identical electronic structure and energy and the observed differences may be due to numerical imprecision.
The reactivity descriptors for quinclorarc were calculated from the values reported in Table 1. It is possible to calculate the value of the vertical electron affinity energy as are the total ground-state energies in the neutral N and singly charged ( ) In a similar way, the ionization potential was calculated as . The values of m, h, w were calculated employing equations 1, 2 and 3 respectively (see Table 2). From the results reported in Table 2, it may be observed that the global hardness decreases when solvent is taken into account. Also, note that the value of w is larger in the gas than in the solution, suggesting that the electrophilic behavior of quinclorac decreases in the aqueous phase. However, it is important to mention that the solvation may change the value of w relative to the gas phase. One of the reasons for this is that a polar solvent like water might stabilize charged species relative to the gas phase, and this would cause an increase in the value of w without having a direct implication on quinclorac reactivity. Therefore, it is possible to conclude that an analysis of quinclorac reactivity based only on the value of w is not definitive.

Frontier molecular orbital analysis (FMO)
The distribution of the electrophilic sites in a system can be derived from the theory of frontier molecular orbitals [67] within the frozen core approximation [34], Eq. 5. To find out the distribution of electrophilic sites, we analyzed the sites where the HOMO frontier orbital attains its larger absolute value on quinclorac (Conformer I). Here it is important to mention that the HOMO and LUMO maps obtained at the B3LYP/6-311++G(2d,2p) and HF/6-311++G(2d,2p) levels were qualitatively similar in the gas and aqueous phases. In Fig. 4a, it may be observed that the greatest extension value of HOMO, calculated at the HF/6-311++G(2d,2p) level, is on the 3C, 10N, 8C, 5C, 16Cl and 11Cl atoms. These atoms are the more susceptible sites to an electrophilic attack. Note that the oxygen atoms are not reactive to this kind of attack according to the HOMO space distribution. In the case of LUMO distributions (see Fig. 4b), the sites where the LUMO attains its larger values are the 7C, 10N, 2C, 3C, 9C, 5C atoms. These sites correspond to the more nucleophilic active sites on quinclorac according to Eq. 6.

Condensed Fukui function
From Fig. 4 it may be observed that the HOMO and LUMO orbitals do not allow a quantitative analysis of    Table 3 we report the values of the Fukui function obtained in the gas phase, as evaluated through NBO natural charges which were calculated at the B3LYP/6-311++G(2d,2p) and MP2/6-311++G(2d,2p) levels of theory. We carried out these calculations in order to compare the results provided by the B3LYP hybrid functional with those given by more sophisticated calculations such as the MP2 theory [61]. From Table 3 it appears that at the DFT level the most susceptible site to an electrophilic attack is located on 16Cl. In the case of a nucleophilic attack, the most reactive site is predicted on 7C, and for a free radical attack the most reactive site is 16Cl. If we analyze the distribution of reactive sites through the condensed Fukui function evaluated based on MP2/6-311++G(2d,2p) theory, we note some differences. For an electrophilic attack the most reactive atom is 7C, for a nucleophilic attack the most reactive site is 2C, and for a free radical attack the most reactive site is 7C.
In Table 4 we report the values of the condensed Fukui Function of quinclorac in the aqueous phase. Note that in the case of the B3LYP method the more reactive sites for the electrophilic and nucleophilic cases are on 11Cl and 7C respectively, while a free radical attack should be aimed at 7C. In the case of MP2 calculations, the electrophilic reactivity is on 8C while the nucleophilic and free radical activity are focused on 3C and 2C, respectively. These results indicate a strong influence of solvent on the reactivity exhibited by quinclorac. We analyzed the effect of the chosen charge calculation scheme on the condensed Fukui function and we have compared the MEP and Mulliken charges (data not shown). At the MP2 level, the values of the Fukui Function evaluated through MEP charges were similar to those obtained from NBO charges: in the aqueous phase, the more reactive sites are located on 8C, 3C and 2C for an electrophilic, nucleophilic and free radical attack respectively. The Fukui function values (FF) calculated by employing Mulliken charges show the larger electrophilic reactivity on 11Cl, while the nucleophilic and free radical activity are focused on 7C and 11Cl. The FF evaluated at the DFT level did not show any consistency and the more reactive sites are different when NBO, MEP or Mulliken charges are employed. From these results, it is possible to conclude that the FF obtained from NBO Natural charges at the MP2 level is a better option to evaluate the Fukui function. Moreover, the FF values obtained with MP2 calculations predict reactivity on 8C, which has been reported as one of the reactive sites in some degradation pathways of quinclorac [20]. Thus, from these results it is possible to suggest that MP2 calculations improve the prediction of Fukui function values in comparison to those obtained from the B3LYP method. From our calculations, it is possible to predict that the most electrophilic site is located on atom 8C. Under these conditions an electrophile would attack at 8C, via an electrophilic addition to the double bond of the benzene ring, followed by the 8C-11Cl scission [74][75][76]. A nucleophilic aromatic attack to 3C might induce an addition-elimination mechanism, where the carboxylic group would be substituted by the nucleophile, while a free radical attack might cause a hydrogen substitution on 2C [74]. However, it must be remembered that the full pathway depends on the type of electrophile, nucleophile or free radical employed during the attacks.
In an aqueous environment at pH ≈ 7 carboxyl groups are usually deprotonated, meaning they lose a H + , and become negatively charged. Therefore the reactivity analysis of the deprotonated (anionic) form of quinclorac (DQ) was studied. The DQ structure was optimized at the MP2/6-311++G(2d,2p) level of theory, see Fig. 5. In Table 5 we report the values of the condensed Fukui Function of DQ in the aqueous phase. Note that the more reactive sites are located on 3C, 7C and 7C for an electrophilic, nucleophilic and free radical attack, respectively. These results suggest that a decarboxylation is expected during an electrophilic
The results indicate that the solvation modifes the values of the reactivity descriptors of quinclorac in the aqueous phase in comparison to the gas phase. The values of the Fukui function at the MP2 level indicated that the more reactive sites of quinclorac in the aqueous phase are located on 8C, 3C and 2C for an electrophilic, nucleophilic and free radical attack, respectively. For the case of deprotonated quinclorac, the more reactive sites are located on 3C, 7C and 7C for electrophilic, nucleophilic and free radical attacks, respectively. Acknowledgement: C.H.R.R. is grateful for a postdoctoral fellowship from CONACYT. We gratefully acknowledge financial support from CONACyT project APOY-COMPL2008 No. 91261 and from the Universidad Autónoma del Estado de Hidalgo in projects PIFI 200813M8U0017T-04-01 and PIFI-2009-13MSU0017T-0401. We are also grateful to the reviewers of the manuscript for valuable suggestions.