Behavior of the E–E’ Bonds (E, E’ = S and Se) in Glutathione Disulfide and Derivatives Elucidated by Quantum Chemical Calculations with the Quantum Theory of Atoms-In-Molecules Approach

The nature of the E–E’ bonds (E, E’ = S and Se) in glutathione disulfide (1) and derivatives 2–3, respectively, was elucidated by applying quantum theory of atoms-in-molecules (QTAIM) dual functional analysis (QTAIM-DFA), to clarify the basic contribution of E–E’ in the biological redox process, such as the glutathione peroxidase process. Five most stable conformers a–e were obtained, after applying the Monte-Carlo method then structural optimizations. In QTAIM-DFA, total electron energy densities Hb(rc) are plotted versus Hb(rc) − Vb(rc)/2 at bond critical points (BCPs), where Vb(rc) are potential energy densities at BCPs. Data from the fully optimized structures correspond to the static nature. Those containing perturbed structures around the fully optimized one in the plot represent the dynamic nature of interactions. The behavior of E–E’ was examined carefully. Whereas E–E’ in 1a–3e were all predicted to have the weak covalent nature of the shared shell interactions, two different types of S–S were detected in 1, depending on the conformational properties. Contributions from the intramolecular non-covalent interactions to stabilize the conformers were evaluated. An inverse relationship was observed between the stability of a conformer and the strength of E–E’ in the conformer, of which reason was discussed.


Introduction
E-E' bonds (E, E' = S and Se) play a crucial role in biological redox processes [1]. High energy levels of HOMO and low energy levels of LUMO of the E-E' bond must be the driving force for the high reactivity in the redox processes. The HOMO and LUMO of E-E' would correspond to n p (E/E') and σ* (E-E'), respectively, where n p (E/E') denote the p-type lone pair orbitals of E and/or E', while σ* (E-E') corresponds to the σ*-orbital of E-E'. Glutathione disulfide (GSSG: 1) has been widely used as a redox reagent in vitro. To facilitate the protein folding process, a mixture of 1 and glutathione (GSH) is often confirmed as the optimum condition if concentrations similar to those observed in vivo [2] are employed [3][4][5][6]. The reduced form of ribonuclease A will undergo disulfide-coupled folding and gain in structural stability in the presence of 1, for example [7]. The detoxification of hydroperoxides in the glutathione peroxidase (GPx) process must be one of most important biological redox processes [8][9][10][11][12][13]. Scheme 1 summarizes a catalytic mechanism proposed for the antioxidant activity of GPx, which is a typical example of the intervention of E-E' (E, E' = S, Se) in biological reactions. According to this mechanism, two equivalents of GSH are oxidized to the corresponding oxidized disulfide in the overall process, while the hydroperoxide is reduced to water [14,15].
reactions. According to this mechanism, two equivalents of GSH are oxidized to the corresponding oxidized disulfide in the overall process, while the hydroperoxide is reduced to water [14,15]. Scheme 1. Catalytic mechanism, proposed for the antioxidant activity of GPx.
The behavior of the S-S, S-Se and Se-Se bonds should be clarified, bearing in mind the role of these bonds in the antioxidant mechanism. It is highly important to elucidate the behavior of the S-S bond in 1 together with S-Se and Se-Se in two derivatives of 1 (compounds 2 and 3, respectively). Scheme 2 illustrates the structures of 1-3. There are many possibilities for the formation of intramolecular hydrogen bonds (HBs) in 1-3, although the intermolecular HBs of the solute-solute and solute-solvent interactions must also be important in the real system. HBs in 1-3 must be considered in assessing the basic properties of 1-3 based on the calculated results, rather than performing calculations for a single molecule in vacuum. Scheme 2 also shows the structures of Rcystine and its derivatives 4-6 and MeEE'Me (compounds 7-9). Scheme 2. Structures of glutathione disulfide (1) and its derivatives (2 and 3) and R-cystine (4) and its derivatives (5 and 6), together with MeEE'Me (7)(8)(9).
The structures of 1 and 4, determined by X-ray crystallographic analysis, have been reported, although 4 is in the di-protonated form. Figure 1 shows these structures. The structure of 1 was observed as a half-extended form close to C2 symmetry with the formation of zwitterions [16]. The structure of 4 was reported as an extended form [17]. The extended form in the observed structure of 4 may be the result of the electrostatic repulsion of the positive charges developed on 4 2+ . Many conformers must exist in such compounds, primarily due to the intramolecular HBs. Scheme 1. Catalytic mechanism, proposed for the antioxidant activity of GPx.
The behavior of the S-S, S-Se and Se-Se bonds should be clarified, bearing in mind the role of these bonds in the antioxidant mechanism. It is highly important to elucidate the behavior of the S-S bond in 1 together with S-Se and Se-Se in two derivatives of 1 (compounds 2 and 3, respectively). Scheme 2 illustrates the structures of 1-3. There are many possibilities for the formation of intramolecular hydrogen bonds (HBs) in 1-3, although the intermolecular HBs of the solute-solute and solute-solvent interactions must also be important in the real system. HBs in 1-3 must be considered in assessing the basic properties of 1-3 based on the calculated results, rather than performing calculations for a single molecule in vacuum. Scheme 2 also shows the structures of R-cystine and its derivatives 4-6 and MeEE'Me (compounds 7-9).
Molecules 2018, 23, 443 2 of 19 reactions. According to this mechanism, two equivalents of GSH are oxidized to the corresponding oxidized disulfide in the overall process, while the hydroperoxide is reduced to water [14,15]. Scheme 1. Catalytic mechanism, proposed for the antioxidant activity of GPx.
The behavior of the S-S, S-Se and Se-Se bonds should be clarified, bearing in mind the role of these bonds in the antioxidant mechanism. It is highly important to elucidate the behavior of the S-S bond in 1 together with S-Se and Se-Se in two derivatives of 1 (compounds 2 and 3, respectively). Scheme 2 illustrates the structures of 1-3. There are many possibilities for the formation of intramolecular hydrogen bonds (HBs) in 1-3, although the intermolecular HBs of the solute-solute and solute-solvent interactions must also be important in the real system. HBs in 1-3 must be considered in assessing the basic properties of 1-3 based on the calculated results, rather than performing calculations for a single molecule in vacuum. Scheme 2 also shows the structures of Rcystine and its derivatives 4-6 and MeEE'Me (compounds 7-9). Scheme 2. Structures of glutathione disulfide (1) and its derivatives (2 and 3) and R-cystine (4) and its derivatives (5 and 6), together with MeEE'Me (7)(8)(9).
The structures of 1 and 4, determined by X-ray crystallographic analysis, have been reported, although 4 is in the di-protonated form. Figure 1 shows these structures. The structure of 1 was observed as a half-extended form close to C2 symmetry with the formation of zwitterions [16]. The structure of 4 was reported as an extended form [17]. The extended form in the observed structure of 4 may be the result of the electrostatic repulsion of the positive charges developed on 4 2+ . Many conformers must exist in such compounds, primarily due to the intramolecular HBs. Reactions of 1 and/or 4 in vivo proceed under conditions with very large and highly complex species. However, the essence of the elementary processes is expected to be close to that of the typical chemical reactions. Therefore, it would be instructive to start with less complex species to clarify the behavior of the E-E' bonds (E, E' = S and Se). We reported the dynamic and static behavior of S-S in R-cystine (4) and S-Se and Se-Se in the derivatives of 4 (5 and 6, respectively), together with MeEE'Me 7-9 as references [18]. It is challenging to clarify the nature of the E-E' bonds (E, E' = S and Se) in glutathione disulfide and its derivatives (1-3), although the structures of 1-3 are considerably more complex relative to 4-6, respectively. Structures 1-3 will have much more plausible HBs than 4-6.
The QTAIM approach, introduced by Bader [19][20][21][22][23][24][25][26][27][28], enables us to analyze the nature of chemical bonds and interactions [11][12][13][14][15][16]. A bond critical point (BCP, * ) is an important concept in QTAIM. BCP is a point where ρ(r) (charge density) reaches a minimum along the interatomic (bond) path, while it is a maximum on the interatomic surface separating the atomic basins. ρ(r) at BCP is denoted by ρb(rc) and other QTAIM functions are denoted in a similar way. Interactions seem to be defined by the corresponding bond paths (BPs), but we must be careful to use the correct terminology with the concept. Interactions would be easily imaged by means of QTAIM if they can be defined as the corresponding BPs, especially for experimental chemists. However, it is demonstrated that the detection of the BPs between two atoms in a molecule emerging from natural alignment of the gradient vector held of the one-electron density of a molecule is neither necessary nor a sufficient condition for the presence of a chemical bond between those atoms [29][30][31][32][33][34]. In this connection, it is pointed out that the terms line paths (LPs) and line critical points (LCPs) should be used in place of BPs and BCPs, respectively [30]. Consequently, the dynamic and static nature in this work should be regarded as the investigation performed at LCPs on LPs corresponding to the E-E' interactions. Nevertheless, the interactions expected for E-E' are clearly detected by BPs with BCPs, which is another reason to use BPs and BCPs in this work. The structures of species can be described by molecular graphs, which are the sets of attractors (atoms) and BPs, together with BCPs, ring critical points (RCPs) and cage critical points (CCPs). We recently proposed QTAIM-DFA [35][36][37][38][39], as a tool for experimental chemists to analyze their own results concerning chemical bonds and interactions using their own images. QTAIM-DFA provides an excellent possibility to evaluate, classify, characterize and understand weak to strong interactions in a unified manner [40][41][42][43][44][45]. Hb(rc) are plotted versus Hb(rc) − Vb(rc)/2 in QTAIM-DFA, where Hb(rc) and Vb(rc) are the total electron energy densities and potential energy densities, respectively, at BCPs. The QTAIM-DFA treatment can incorporate the classification of interactions based on the signs of Hb(rc) and ∇ 2 ρb(rc) (Laplacian ρ), as shown in Scheme S1 of the Supplementary Materials, since the signs of Hb(rc) − Vb(rc)/2 must be equal to those of ∇ 2 ρb(rc), where (ћ 2 /8m) ∇ 2 ρb(rc) = Hb(rc) − Vb(rc)/2 (= Gb(rc) + Vb(rc)/2 while Hb(rc) = Gb(rc) + Vb(rc), Gb(rc): kinetic energy densities at BCPs) (see Equations (S1), (S2) and (S2') of the Supplementary Materials). In our treatment, data for the perturbed structures around fully optimized structures are employed for the plots in addition to data for the fully optimized structures [35][36][37][38][39]. Data from the fully optimized structures are analyzed by the polar coordinate (R, θ) representation, which The QTAIM approach, introduced by Bader [19][20][21][22][23][24][25][26][27][28], enables us to analyze the nature of chemical bonds and interactions [11][12][13][14][15][16]. A bond critical point (BCP, * ) is an important concept in QTAIM. BCP is a point where ρ(r) (charge density) reaches a minimum along the interatomic (bond) path, while it is a maximum on the interatomic surface separating the atomic basins. ρ(r) at BCP is denoted by ρ b (r c ) and other QTAIM functions are denoted in a similar way. Interactions seem to be defined by the corresponding bond paths (BPs), but we must be careful to use the correct terminology with the concept. Interactions would be easily imaged by means of QTAIM if they can be defined as the corresponding BPs, especially for experimental chemists. However, it is demonstrated that the detection of the BPs between two atoms in a molecule emerging from natural alignment of the gradient vector held of the one-electron density of a molecule is neither necessary nor a sufficient condition for the presence of a chemical bond between those atoms [29][30][31][32][33][34]. In this connection, it is pointed out that the terms line paths (LPs) and line critical points (LCPs) should be used in place of BPs and BCPs, respectively [30]. Consequently, the dynamic and static nature in this work should be regarded as the investigation performed at LCPs on LPs corresponding to the E-E' interactions. Nevertheless, the interactions expected for E-E' are clearly detected by BPs with BCPs, which is another reason to use BPs and BCPs in this work. The structures of species can be described by molecular graphs, which are the sets of attractors (atoms) and BPs, together with BCPs, ring critical points (RCPs) and cage critical points (CCPs). We recently proposed QTAIM-DFA [35][36][37][38][39], as a tool for experimental chemists to analyze their own results concerning chemical bonds and interactions using their own images. QTAIM-DFA provides an excellent possibility to evaluate, classify, characterize and understand weak to strong interactions in a unified manner [40][41][42][43][44][45]. H b (r c ) are plotted versus H b (r c ) − V b (r c )/2 in QTAIM-DFA, where H b (r c ) and V b (r c ) are the total electron energy densities and potential energy densities, respectively, at BCPs. The QTAIM-DFA treatment can incorporate the classification of interactions based on the signs of H b (r c ) and ∇ 2 ρ b (r c ) (Laplacian ρ), as shown in Scheme S1 of the Supplementary Materials, since the signs of H b (r c ) − V b (r c )/2 must be equal to those of kinetic energy densities at BCPs) (see Equations (S1), (S2) and (S2') of the Supplementary Materials). In our treatment, data for the perturbed structures around fully optimized structures are employed for the plots in addition to data for the fully optimized structures [35][36][37][38][39]. Data from the fully optimized structures are analyzed by the polar coordinate (R, θ) representation, which corresponds to the static nature of interactions. Data from the perturbed structures and a fully optimized structure are used to construct a curve. Each curve is analyzed in terms of the (θ p , κ p ) parameters: θ p corresponds to the tangent line of the plot and κ p is the curvature. We proposed the concept of the "dynamic nature of interactions" based on (θ p , κ p ) [35][36][37][38][39]. QTAIM-DFA is applied to typical chemical bonds and interactions. Rough criteria have been established that distinguish the chemical bonds and interactions in question from others. QTAIM-DFA and the criteria are explained in the Supplementary Materials, employing Schemes S1 and S2, Figure S1 and Equations (S1)-(S7). The basic concept of the QTAIM approach is also described in the Supplementary Materials.
The behavior of E-E' (E, E' = S and Se) in 1-3 is expected to be related to that in the glutathione peroxidase (GPx) process. The dynamic and static nature of E-E' in 1-3 is elucidated by applying QTAIM-DFA. The same method is applied to E-E' in 4-6 and 7-9 to reexamine the nature of these bonds. We present the results of the theoretical elucidation of the nature of the E-E' bonds in 1-6 with QTAIM-DFA to better understand the role of E-E' in the antioxidant activity of GPx. Quantum chemical (QC) calculations are also applied to examine the structural features of 1-6. The E-E' bonds in 1-6 are classified and characterized by employing the criteria and the behavior of the bonds in 7-9 as references.

Methodological Details in Calculations
The structures were optimized employing the Gaussian 09 programs [46] unless otherwise noted. For each species 1-6, five conformers were optimized with the 6-311+G(3d) basis sets for S and Se, and with the 6-311++G(d, p) basis sets for O, N, C and H [47][48][49][50]. The basis set system is called BSS-A in this paper. The DFT level of M06-2X [51] is applied to the calculations. Before the final optimizations, the full conformer search with the Monte-Carlo method in Spartan 02 [52] was applied to each of 1-6. At least six thousand and five hundred conformers were generated for each of 1-3 with the Merck Molecular Force Field (MMFF) method [53]. The most stable thirty independent conformers from the Monte-Carlo method were optimized using the 3-21G basis sets at the B3LYP level [54,55] for each of 1-3. Next, the most stable fifteen conformers were optimized with M06-2X/6-31G(d) for each conformer, as predicted with B3LYP/3-21G of the Gaussian09 program. The most stable five conformers from M06-2X/6-31G(d) were further optimized with BSS-A at the M06-2X level (M06-2X/BSS-A). The final five optimized conformers were confirmed by the frequency analysis for each of 1-3. These five conformers are called a, b, c, d and e, where conformer a is the most stable among the five, followed by b, then c, then d and then e. In the case of 4-6, seven hundred and twenty conformers were generated for each species with the PM3 method [56]. Similar to the case for 1-3, the five most stable conformers (a-e) were determined for each of 4-6. The structures of 7 and 9 are optimized retaining the C 2 symmetry, while that of 8 is retaining the C 1 symmetry. The population analysis has also been performed by the natural bond orbital method [57] at M06-2X/BSS-A level of theory using natural bond orbital (NBO) program [58].
QTAIM functions were calculated using the Gaussian 09 program package at the same level of DFT theory (M06-2X/BSS-A), and the data were analyzed with the AIM2000 program [19,59]. Normal coordinates of internal vibrations (NIV) obtained by the frequency analysis were employed to generate the perturbed structures [38,39]. This method is called NIV and explained in Equation (1). The k-th perturbed structure in question (S kw ) was generated by the addition of the normal coordinates of the k-th internal vibration (N k ) to the standard orientation of a fully optimized structure (S o ) in the matrix representation. The coefficient f kw in Equation (1) controls the difference in structures between S kw and S o : f kw is determined to satisfy Equation (2) for an interaction in question, where r and r o show the interaction distances in question in the perturbed and fully optimized structures, respectively, with a o representing the Bohr radius (0.52918 Å) [12,13]. The perturbed structures with NIV correspond to those where r has been elongated or shortened by 0.05a o or 0.1a o , relative to r o , in the fully optimized structures [60]. The selected motion must be most effectively localized on the interaction in question among the zero-point internal vibrations. N k of five digits are used to predict S kw : for data from the five points of w = 0, ±0.05 and ±0.1 in Equation (2) in QTAIM-DFA. Each plot is analyzed using a regression curve of the cubic function as shown in Equation (3)
Species The r(S, S) values for 1a-1e are predicted to be larger than those of 7. The differences in r(S, S) for 1a-1e (∆r(S, S: 1x) = r(S, S: 1x) − r(S, S: 7), where x = a-e) are 0.02 Å < ∆r(S, S: 1x) < 0.03 Å for 1a-1c, ∆r(S, S: 1x) < 0.01 Å for 1d and ∆r(S, S: 1x) ≈ 0.20 Å for 1e. Similarly, the ∆r(E, E') values are less than or very close to 0.01 Å for na-ne (n = 2-6), except for ∆r(E, E') ≈ 0.015 Å for 2e, 4a, 5d and 5e with ∆r(E, E') ≈ 0.03 Å for 3d. There must be a specific reason for the unexpectedly large value of 0.20 Å for ∆r(E, E': 1e). Three S, S and O atoms align linearly in 1e, which is explained by assuming the formation of hypervalent interactions of S 2 O σ(3c-4e) of the n p (O)→σ*(S-S) type (see Figure 4). In this interaction, σ*(S-S) accepts electrons from n p (O). As a result, the S-S bond must be unexpectedly elongated relative to the usual length, and the O-S distance will be substantially shortened relative to the sum of the vdW radii. The O-S distance is predicted to be 2.7714 Å, which is shorter than the sum of the vdW radii by 0.55 Å. The S 2 O σ(3c-4e) model explains the predicted result for 1e, reasonably well.
In the case of φ(CEE'C) (= φ A ), the values for 1a-6e from the corresponding values of 7-9 are given by ∆φ , where n = 1-6; x = a-e; and E, E' = S and Se. The absolute values of φ A will be used to estimate ∆φ A . The magnitudes of the values are ∆φ A (E, E': nx) ≈ 58 • for 3d, 31 • < ∆φ A (E, E': nx) < 35 • for 1a-1c and 1e, ∆φ A (E, E': nx) ≈ 25 • for 2b, and 10 • ≤ ∆φ A (E, E': nx) ≤ 20 • for 1d, 3e, 4e and 6d. The magnitudes of ∆φ A (E, E': nx) are less than 10 • (−10 • ≤ ∆φ A (E, E': nx) ≤ 10 • ) for others, except for ∆φ A (E, E': nx) ≈ -12 • for 5d and 6d and ∆φ A (E, E': nx) ≈ −20 • for 2e, 4a and 5e. The results must be the reflection from the easy deformation in φ A . To evaluate the energy for the deformation in φ A , 7-9 were optimized assuming φ A = 0 • and 180 • , in addition to the fully optimized structures (85 • ≤ φ A ≤ 86 • ). They were optimized to be 7 (C 2v ), 8 (C s ) and 9 (C 2v ) at φ A = 0 • and 7 (C 2h ), 8 (C s ) and 9 (C 2h ) at φ A = 180 • . In the case of 9, the structures were further optimized with φ A fixed every 15 The results are summarized in Table S1 of the Supplementary Materials. Figure 2 shows the plot of the energies for the optimized structures versus φ A . The energy seems less than 15 kJ mol −1 for 45 The energy for the deformation of φ A in 7 and 8 seems comparable to that in 9. The very easy deformation in φ A is well demonstrated, exemplified by 7-9, which supports the results shown in Table 1. Such easy deformation in φ A is also reported for some dichalcogenides [62]. The r(S, S) values for 1a-1e are predicted to be larger than those of 7. The differences in r(S, S) for 1a-1e (∆r(S, S: 1x) = r(S, S: 1x) − r(S, S: 7), where x = a-e) are 0.02 Å < ∆r(S, S: 1x) < 0.03 Å for 1a-1c, ∆r(S, S: 1x) < 0.01 Å for 1d and ∆r(S, S: 1x) ≈ 0.20 Å for 1e. Similarly, the ∆r(E, E') values are less than or very close to 0.01 Å for na-ne (n = 2-6), except for ∆r(E, E') ≈ 0.015 Å for 2e, 4a, 5d and 5e with ∆r(E, E') ≈ 0.03 Å for 3d. There must be a specific reason for the unexpectedly large value of 0.20 Å for ∆r(E, E': 1e). Three S, S and O atoms align linearly in 1e, which is explained by assuming the formation of hypervalent interactions of S2O σ(3c-4e) of the np(O)→σ*(S-S) type (see Figure 4). In this interaction, σ*(S-S) accepts electrons from np(O). As a result, the S-S bond must be unexpectedly elongated relative to the usual length, and the O---S distance will be substantially shortened relative to the sum of the vdW radii. The O---S distance is predicted to be 2.7714 Å, which is shorter than the sum of the vdW radii by 0.55 Å. The S2O σ(3c-4e) model explains the predicted result for 1e, reasonably well. for 2e, 4a and 5e. The results must be the reflection from the easy deformation in φA. To evaluate the energy for the deformation in φA, 7-9 were optimized assuming φA = 0° and 180°, in addition to the fully optimized structures (85° ≤ φA ≤ 86°). They were optimized to be 7 (C2v), 8 (Cs) and 9 (C2v) at φA = 0° and 7 (C2h), 8 (Cs) and 9 (C2h) at φA = 180°. In the case of 9, the structures were further optimized with φA fixed every 15° for 0° ≤ φA ≤ 180°. The results are summarized in Table S1 of the Supplementary Materials. Figure 2 shows the plot of the energies for the optimized structures versus φA. The energy seems less than 15 kJ mol -1 for 45° ≤ φA ≤ 135° in 9. The energy for the deformation of φA in 7 and 8 seems comparable to that in 9. The very easy deformation in φA is well demonstrated, exemplified by 7-9, which supports the results shown in Table 1. Such easy deformation in φA is also reported for some dichalcogenides [62].

Structural
Feature of 1a-6a and 7-9 Figure 3 illustrates the molecular graphs of 1a-3a, drawn on the optimized structures, together with the optimized structures containing the non-covalent interactions. Figure 4 shows the molecular graphs of 1b-1e, drawn on the optimized structures. Molecular graphs of 2b-2e and 3b-3e are drawn in Figures S2 and S3 of the Supplementary Materials, respectively. Figure 5 illustrates the molecular graphs of 4a, 5a, 6a and 4b-4e, drawn on the optimized structures and the optimized structure. Molecular graphs of 5b-5e, 6b-6e and 7-9 are drawn in Figures S4-S6 of the Supplementary Materials, respectively.

Structural
Feature of 1a-6a and 7-9 Figure 3 illustrates the molecular graphs of 1a-3a, drawn on the optimized structures, together with the optimized structures containing the non-covalent interactions. Figure 4 shows the molecular graphs of 1b-1e, drawn on the optimized structures. Molecular graphs of 2b-2e and 3b-3e are drawn in Figures S2 and S3 of the Supplementary Materials, respectively. Figure 5 illustrates the molecular graphs of 4a, 5a, 6a and 4b-4e, drawn on the optimized structures and the optimized structure. Molecular graphs of 5b-5e, 6b-6e and 7-9 are drawn in Figures S4-S6 of the Supplementary Materials, respectively.
The structural features of 7-9 are described, first. Only classical chemical bonds are detected in the molecular graphs of 7-9, as shown in Figure S6 Figure 5). The numbers of the intramolecular non-covalent interactions are counted separately by the interaction types, which are differentiated by the colors. The results are collected in Table S2 of the Supplementary Materials. Indeed, the stability of the conformers is expected to relate to the numbers, but they must be stabilized by the total energy of the intramolecular non-covalent interactions. The C-H-X interactions are also detected, however, they are neglected, since they would not make a significant contribution to stabilizing the conformers. The structural features of 7-9 are described, first. Only classical chemical bonds are detected in the molecular graphs of 7-9, as shown in Figure S6 Figure 5). The numbers of the intramolecular non-covalent interactions are counted separately by the interaction types, which are differentiated by the colors. The results are collected in Table S2 of the Supplementary Materials. Indeed, the stability of the conformers is expected to relate to the numbers, but they must be stabilized by the total energy of the intramolecular non-covalent interactions. The C-H---X interactions are also detected, however, they are neglected, since they would not make a significant contribution to stabilizing the conformers.       The molecular graphs for 1a-3e are very complex, with many intramolecular non-covalent interactions. An effort is made to classify the interactions in 1a-3e, as in 4a-6e. The non-covalent interactions appearing in the molecular graphs of 1a-3e are similarly drawn on the optimized structures, as shown in Figure 3. The numbers of the intramolecular non-covalent interactions in 1a-3e are also counted separately based on the type of interaction. The results are collected in Table S2 of the Supplementary Materials. The numbers seem to be correlated to the stability of the conformers. However, it must be difficult to estimate numerically the stability of the conformers based on the numbers. The stability must be controlled by the total energy of the intramolecular interactions.
Nevertheless, it is very important to understand how E rel for the conformers are determined by the intramolecular non-covalent interactions in the conformers of 1a-3e, as a whole. How can E rel be evaluated based on the contributions from the non-covalent interactions? We searched for a method to evaluate the stability of the conformers based on the overall intramolecular non-covalent interactions. Then, we devised a method to evaluate the contributions, which is discussed next.

Factors Determining the Relative Energies of 1a-6a
The proposed method is explained in Scheme 3 with Equations (5)- (7). The evaluation process is as follows: (i) GEE'G is fully optimized; (ii) E and E' in the optimized GEE'G are replaced by H and H; (iii) The structural parameters for the two replaced H atoms are (partially) optimized, with other atoms fixed at the fully optimized geometry; (iv) The structural parameters of CEE'C are fixed at the fully optimized positions, and H atoms are added on each side of CEE'C in place of the organic ligands to give H 3 CEE'CH 3 ; (v) The structural parameters of the six H atoms are optimized: Molecules 2018, 23, 443 9 of 19 The molecular graphs for 1a-3e are very complex, with many intramolecular non-covalent interactions. An effort is made to classify the interactions in 1a-3e, as in 4a-6e. The non-covalent interactions appearing in the molecular graphs of 1a-3e are similarly drawn on the optimized structures, as shown in Figure 3. The numbers of the intramolecular non-covalent interactions in 1a-3e are also counted separately based on the type of interaction. The results are collected in Table S2 of the Supplementary Materials. The numbers seem to be correlated to the stability of the conformers. However, it must be difficult to estimate numerically the stability of the conformers based on the numbers. The stability must be controlled by the total energy of the intramolecular interactions.
Nevertheless, it is very important to understand how Erel for the conformers are determined by the intramolecular non-covalent interactions in the conformers of 1a-3e, as a whole. How can Erel be evaluated based on the contributions from the non-covalent interactions? We searched for a method to evaluate the stability of the conformers based on the overall intramolecular non-covalent interactions. Then, we devised a method to evaluate the contributions, which is discussed next.

Factors Determining the Relative Energies of 1a-6a
The proposed method is explained in Scheme 3 with Equations (5)- (7). The evaluation process is as follows: (i) GEE'G is fully optimized; (ii) E and E' in the optimized GEE'G are replaced by H and H; (iii) The structural parameters for the two replaced H atoms are (partially) optimized, with other atoms fixed at the fully optimized geometry; (iv) The structural parameters of CEE'C are fixed at the fully optimized positions, and H atoms are added on each side of CEE'C in place of the organic ligands to give H3CEE'CH3; (v) The structural parameters of the six H atoms are optimized: The method shown in Scheme 3 will evaluate the intramolecular G---G interactions and the deformation energies around CEE'C but not the steric factor around the E-E' moiety. Equation (5) shows the relationship between E(GEE'G) for the fully optimized structure and the energies evaluated by the proposed method, where α shows the errors in energy between E(GEE'G) and the components, which contains the steric factor around the E-E' moiety. The relationship for Erel(GEE'G) is shown in Equation (6), where E(CH4) disappears. As shown in Equation (7)  The method shown in Scheme 3 will evaluate the intramolecular G-G interactions and the deformation energies around CEE'C but not the steric factor around the E-E' moiety. Equation (5) shows the relationship between E(GEE'G) for the fully optimized structure and the energies evaluated by the proposed method, where α shows the errors in energy between E(GEE'G) and the components, which contains the steric factor around the E-E' moiety. The relationship for E rel (GEE'G) is shown in Equation (6), where E(CH 4 ) disappears. As shown in Equation (7), E rel (GEE'G) opt can be approximated as 2E rel (G fix -H opt ) + E rel [(H opt ) 3 C fix E fix -E' fix C fix (H opt ) 3 ] if α is almost constant. The E rel values are given as the values from the most stable conformers in 1a-6a, if applied to 1a-6e, respectively. The results of the calculations for 1a-6e are collected in Table S3 3 ] are abbreviated as E rel (2GH) p-opt and E rel (MeSSMe) p-opt , respectively (p-opt: partially optimizations). Figure 6 shows the plot of E rel (2GH + MeEE'Me) p-opt for 1a-3e, together with E rel (GEE'G) opt . The E rel (2GH + MeEE'Me) p-opt values seem to match the E rel (GEE'G) opt values for 1a-3e, except for 2b and 2d. Indeed, the relative stabilities of the conformers for the S-S and Se-Se species are explained well by the treatment, but they are not explained well by treatment for the S-Se species, especially for 2b and moderately for 2d. Other factors, such as the steric factor around the S-Se moiety, could be important in this case. The very large magnitude of φ A in 2b (110.1 • ) relative to other species (65.0-85.6 • ) is responsible for this result. The deviation in 2b, due to E rel (2GH + MeSSeMe) p-opt (−18.4 kJ mol −1 ) versus E rel (GSSeG) opt (1.0 kJ mol −1 ), is due to the high stability of E rel (2GH) p-opt (-22.7 kJ mol −1 ), which is due to the reflection of the C-H optimizations in 2GH from the unstable position of C-H by φ A for 2b (110.1 • ). The smaller magnitudes in E rel for 2a-2e may lead to greater deviations (see Table S3 of the Supplementary Materials). A similar plot for 4a-6e is shown in Figure S7 of the Supplementary Materials. The relationship between E rel (2CysH + MeEE'Me) p-opt and E rel (CysEE'Cys) opt seems unclear for 4a-6e. After clarification of the structural features of 1a-6e, the contour plots and negative Laplacians are examined next.
Molecules 2018, 23,443 10 of 19 Figure 6 shows the plot of Erel(2GH + MeEE'Me)p-opt for 1a-3e, together with Erel(GEE'G)opt. The Erel(2GH + MeEE'Me)p-opt values seem to match the Erel(GEE'G)opt values for 1a-3e, except for 2b and 2d. Indeed, the relative stabilities of the conformers for the S-S and Se-Se species are explained well by the treatment, but they are not explained well by treatment for the S-Se species, especially for 2b and moderately for 2d. Other factors, such as the steric factor around the S-Se moiety, could be important in this case. The very large magnitude of φA in 2b (110.1°) relative to other species (65.0°-85.6°) is responsible for this result. The deviation in 2b, due to Erel(2GH + MeSSeMe)p-opt (−18.4 kJ mol -1 ) versus Erel(GSSeG)opt (1.0 kJ mol -1 ), is due to the high stability of Erel(2GH)p-opt (-22.7 kJ mol -1 ), which is due to the reflection of the C-H optimizations in 2GH from the unstable position of C-H by φA for 2b (110.1°). The smaller magnitudes in Erel for 2a-2e may lead to greater deviations (see Table S3 of the Supplementary Materials). A similar plot for 4a-6e is shown in Figure S7 of the Supplementary Materials. The relationship between Erel(2CysH + MeEE'Me)p-opt and Erel(CysEE'Cys)opt seems unclear for 4a-6e. After clarification of the structural features of 1a-6e, the contour plots and negative Laplacians are examined next.  Figure 7 shows the contour plots of ρ(r), exemplified by 1a-1e, which are drawn on an SSC plane of 1a-e. The plots show that each BCP on E- * -E' exists at the three-dimensional saddle point of ρ(r).  Figure 7 shows the contour plots of ρ(r), exemplified by 1a-1e, which are drawn on an SSC plane of 1a-e. The plots show that each BCP on E- * -E' exists at the three-dimensional saddle point of ρ(r). Figure 8 illustrates the Negative Laplacian, exemplified by 1a-1e, similarly drawn on an SSC plane. All BCP on E- * -E' of 1a-1e exist in the red area of the plots, which means that the BCPs are all in the range of ∇ 2 ρ b (r c ) < 0. Therefore, E- * -E' of 1a-1e are classified by SS (shared shell) interactions. (See also Figure S8

Application of QTAIM-DFA to the E-E' Bonds in 1a-6e
QTAIM functions are calculated for 1a-6e and 7-9. Table 2 collects the results for 1a-3e and 7-9. Figure 9 shows the plots of Hb(rc) versus Hb(rc) − Vb(rc)/2 for 1a-1d and 7, where the data for 1e do not appear in the plotted area. Figure 10

Application of QTAIM-DFA to the E-E' Bonds in 1a-6e
QTAIM functions are calculated for 1a-6e and 7-9. Table 2 collects the results for 1a-3e and 7-9. Figure 9 shows the plots of H b (r c ) versus H b (r c ) − V b (r c )/2 for 1a-1d and 7, where the data for 1e do not appear in the plotted area. Figure 10 Table 2. QTAIM-DFA Parameters and QTAIM Functions at BCPs for the E- * -E bonds in 1a-3e and 7-9, 1 together with the frequencies (ν) and force constants (k f ) corresponding to the E- * -E species in question.   Table 2 and the perturbed structures of 1a-3e and 7-9 are plotted in Figure S9      The QTAIM-DFA parameters of (R, θ) and (θp, κp) are also collected in Table 2, together with the frequencies (ν) and force constants (kf) corresponding to the E- * -E′ bonds in question. While the data for 4a-4e are plotted in Figure S10, those for 5a-6e are in Figure S11 of the Supplementary Materials, together with those for 7-9. Similarly, the plots are analyzed to give the QTAIM-DFA parameters of (R, θ) and (θp, κp). The parameters are collected in Table S4 of the Supplementary Materials, together with the frequencies (νn) and force constants (kf) corresponding to the E- * -E′ bonds in question. The QTAIM-DFA parameters of (R, θ) and (θ p , κ p ) are also collected in Table 2, together with the frequencies (ν) and force constants (k f ) corresponding to the E- * -E bonds in question. While the data for 4a-4e are plotted in Figure S10, those for 5a-6e are in Figure S11 of the Supplementary Materials, together with those for 7-9. Similarly, the plots are analyzed to give the QTAIM-DFA parameters of (R, θ) and (θ p , κ p ). The parameters are collected in Table S4 of the Supplementary Materials, together with the frequencies (ν n ) and force constants (k f ) corresponding to the E- * -E bonds in question.

Nature of the E-E' Bonds in 1a-6e
The E- * -E bonds in 1a-6e are classified and characterized based on R, θ and θ p values, employing those of the standard interactions given in Scheme S2 of the Supplementary Materials as a reference. Before discussing the nature of the E- * -E' bonds in 1a-6e and 7-9, it would be instructive to survey the related criteria. Interactions will be classified as SS or CS interactions for θ > 180 • or θ < 180 • , respectively, which correspond to H b (r c ) − V b (r c )/2 < 0 and H b (r c ) − V b (r c )/2 > 0, respectively. The θ p values play an important role in characterizing the interactions. For SS interactions with θ > 180 • , θ p > 190 • is tentatively assigned, where θ p = 190 • corresponds to θ = 180 • for typical interactions. The covalent interactions will be sub-divided depending on the values of R. The (classical) covalent interactions will be called strong (Cov strong ) if R > 0.15 au; therefore, they should be called weak (Cov weak ) when R < 0.15 au.
The R value of 0.076 au (<0.15 au) is predicted for MeS- * -SMe (7), and those for S- * -S, S- * -Se and Se- * -Se in 1a-6e, 8 and 9, examined in this work, are all less than 0.076 au. Therefore, the Cov strong interactions are not detected in this work. As shown in Table 2, the (θ, θ p ) values for S- * -S in 1a-1e are (188.8-189.6 • , 197.2-197.6 • ) for 1a-1d with (181.9 • , 193.8 • ) for 1e. The value for 1e is apparently smaller than those for 1a-1d due to the elongation of S- * -S by the formation of S 2 O σ(3c-4e) in 1e. This means that S- * -S in 1e should be (much) weaker than those in 1a-1d. Nevertheless, the S- * -S interaction in 1e is classified by the SS interaction and characterized as Cov weak in nature (SS/Cov weak ). All S- * -S interactions in 1a-1d are, of course, predicted to have the nature of (SS/Cov weak ). The (θ, θ p ) values for S- * -Se in 2a-2e are (184.7-184.9 • , 188.1-188.7 • ). Therefore, the S- * -Se interactions are also classified by the SS interactions and characterized as Cov weak in nature (SS/Cov weak ), although the θ p values are slightly less than 190 • . In the case of Se- * -Se in 3a-3e, the (θ, θ p ) values are (185.9-186.4 • , 189.0-189.8 • ). The Se- * -Se interactions are predicted to have the nature of (SS/Cov weak ), similar to the cases of 1a-2e. Note that S- * -S in 1e is predicted to be weaker than S- * -Se in 2a-2e and Se- * -Se in 3a-3e by R and θ, although the trend is reversed for θ p . Indeed, Se- * -Se in 3a-3e is predicted to be stronger than S- * -Se in 2a-2e by θ and θ p , and the inverse trend is true for R. The E- * -E bonds in 4a-6e are all predicted to have the nature of (SS/Cov weak ), as are the bonds in 7-9, similar to the case for 1a-3e.

Factors that Stabilize the E-E' Bonds and the Conformers
The S-S bonds of 1a-1e are predicted to be less stable than that of 7. The S- * -S bond in 1a-1e is predicted to be weaker in the order shown in Equation (8), where 1e is significantly destabilized due to the elongation by the formation of S 2 O σ(3c-4e). The order for the strength of S- * -S seems to exhibit almost the reverse trend of stability compared with the conformers. A similar order is predicted for S- * -Se of 2a-2e with 8. Equation (9) shows the order for S- * -Se, where 2a seems substantially destabilized. On the other hand, the trend is not so clear for Se- * -Se in 3a-3e. The predicted order for Se- * -Se is given in Equation (10). The strength of E- * -E' seems to show a trend with almost inverse stability compared with the conformers containing E- * -E', as mentioned above. While the trend seems rather clear for 1a-1e with 7 and 2a-2e with 8, the trend seems unclear for 3a-3e with 9. This trend would be clear if the data were plotted in a narrow range, whereas it would not be clear if they were plotted in a wider range, although the mechanism is not clear: Se- * -Se in 3c > 9 > 3b > 3a > 3e > 3d S- * -S in 7 > 4b > 4c > 4d > 4e >> 4a (11) S- * -Se in 8 ≈ 5b > 5a ≥ 5c ≥ 5d > 5e (12) Se- * -Se in 9 > 6e > 6c > 6b > 6d > 6a (13) In the case of S-S in 4a-4e with 7, the S- * -S bond becomes less stable in the order shown in Equation (11). The order for the strength of S- * -S is almost the reverse of the stability of the conformers, with the species are divided into four groups of 7, 4b-4d, 4e and 4a. The order for the strength of Se- * -Se is also almost inverse, with the species divided into four groups of 9, 6b-6d, 6e and 6a, as shown in Equation (13). However, the trend in S- * -Se is not as clear for 5a-5e with 8, as predicted in Equation (12). As shown in Figures S10 and S11 of the Supplementary Materials, the data for 4a-4e with 7 are plotted in a narrow range, as are those for 6a-6e with 9, which seems to exhibit a clear trend. However, the data for 5a-5e with 8 are plotted over a wider range, and the trend seems unclear, similar to the cases of 1a-3e, although the mechanism remains unclear.
The trends shown in Equations (8)- (13) are also confirmed through the analysis of ρ b (r c ) and bond orders evaluated based on the natural atomic orbitals, for E-E in 1a-6e and 7-9. The plot of ρ b (r c ) versus the bond orders is shown in Figure S12 of the Supplementary Materials. See also Table S5 of the Supplementary Materials for bond orders of 1a-6e and 7-9.
While the results could be explained in a variety of ways, our explanation is as follows: the intramolecular attractive interactions in 1a-6c stabilize the species, but E-E' would be destabilized through the distortion. This is because the E-E' bonds operate to relax the excess deformation brought about by intramolecular attractive interactions such as HBs. This destabilization would increase the stability of the species. As a result, the E-E' bonds will be predicted to be less stable if they exist in more stable species. The E-E' bonds could be predicted to be rather stable, if the intramolecular attractive interactions do not affect the structures around the E-E' bonds.
The nature of the E-E' bonds in 1-9 is well-described by applying QTAIM-DFA based on the dynamic nature with (θ p , κ p ) and the static nature with (R, θ).

Conclusions
The dynamic and static nature of E-E' (E, E' = S and Se) in glutathione disulfide (1) and derivatives 2-3, respectively, is elucidated by applying QTAIM-DFA together with R-cystine and its derivatives (4-6) and MeEE'Me (7)(8)(9). Five conformers a-e for each of 1-6 are optimized with M06-2X/BSS-A. The conformers are called 1a-1e, which are defined to satisfy E(1a) < E(1b) < E(1c) < E(1d) < E(1e), for example. No intramolecular non-covalent interactions are detected in 7-9, while many such interactions operate to stabilize 1a-3e. Among such interactions, the formation of S 2 O σ(3c-4e) of the n p (O)→σ*(S-S) type detected in 1e elongates r(S, S) by 0.20 Å. The contribution from the intramolecular non-covalent interactions to the stability of the conformers of GEE'G is estimated by separately calculating the G-G part as 2G-H and the E-E' part (as MeE-E'Me), under suitable conditions. The E rel (2GH + MeEE'Me) values explain the E rel values well for 1a-1e and 3a-3e but not so for 2a-2e.
QTAIM-DFA is applied to E- * -E' in 1a-6e and 7-9 by plotting H b (r c ) versus H b (r c ) − V b (r c )/2 for the data from the fully optimized structures and the perturbed structures at BCPs. The QTAIM-DFA parameters of (R, θ) and (θ p , κ p ) are obtained by analyzing the plots. The (θ, θ p , R) values for S- * -S in 1a-1e are (188.8-189.6 • , 197.2-197.6 • , 0.0673-0.0744 au) for 1a-1d with (181.9 • , 193.8 • , 0.0343 au) for 1e. The values for 1e are apparently smaller than those for 1a-1d, due to the elongation of S- * -S by the formation of S 2 O σ(3c-4e) in 1e. Nevertheless, the S- * -S interactions in 1a-1e are all predicted to have the (SS/Cov weak ) nature. Similarly, the E-E' interactions in 2a-6e and 7-9 are all predicted to have the (SS/Cov weak ) nature. The S-S bonds of 1a-1e are predicted to be less stable than those of 7, and the S- * -S bond in 1a becomes weaker than those in 1b, 1d and 7, although 1a is the most stable among 1a-1e. This inverse trend between the stability of the conformers and the strength of E- * -E' is widely observed. The intramolecular non-covalent attractive interactions in 1a-6e stabilize the species but destabilize the E-E' bonds through distortion, where the E-E' bonds act to relax the excess deformation caused by the attractive interactions. These predictions of bond behavior help to understand the reactivity of E-E' in chemical and biological processes.
Supplementary Materials: The following are available online at www.mdpi.com/1420-3049/23/2/443/s1, Scheme S1: QTAIM-DFA: Plot of H b (r c ) versus H b (r c ) − V b (r c )/2 for weak to strong interactions, Scheme S2: Rough classification of interactions by θ and θ p , together with k b (r c ) (= V b (r c )/G b (r c )), Equations (S1)-(S7), Table S1: The energies for 7-9 of the optimized structures and partially optimized structures with φ A , fixed suitably, with M06-2X/BSS-A, Table S2 Table S4: QTAIM-DFA parameters and QTAIM functions at BCPs for the E-E' bonds in 4a-6e and 7-9, a together with the frequencies (ν) and force constants (k f ), corresponding to the E- * -E bonds in question, Table S5: NAO bond orders for 1a-6e and 7-8, Figure S1: Polar (R, θ) coordinate representation of H b (r c ) versus H b (r c ) − V b (r c )/2, with (θ p , κ p ) parameters, Figure S2: Molecular graphs of 2b-2e, drawn on the optimized structures, Figure S3: Molecular graphs of 3b-3e, drawn on the optimized structures, Figure S4: Molecular graphs of 5b-5e, drawn on the optimized structures, Figure S5: Molecular graphs of 6b-6e, drawn on the optimized structures, Figure S6: Molecular graphs of 7-9, drawn on the optimized structures, Figure S7: Plots of E rel of REE'R and (2R-H + MeEE'Me) for 4a-e (CysSSCys), 5a-e (CysSSeCys) and 6a-e (CysSeSeCys), evaluated with M06-2X/BSS-A, Figure S8: Trajectory plots of ρ b (r c ) drawn on the S-S-C planes of 1a-1e, similarly to the case of Figure 6 in the text. Color and marks are same as those in Figure 6, Figure S9: Plots of H b (r c ) versus H b (r c ) − V b (r c )/2 for 1a-3e and 7-9, Figure S10: Plots of H b (r c ) versus H b (r c ) -V b (r c )/2 for 4a-4e and 7, Figure S11: Plots of H b (r c ) versus H b (r c ) − V b (r c )/2 for 5a-6e and 7-8, Figure S12: Plots of ρ b (r c ) versus NAO bond orders for 1a-6e and 7-8 and Cartesian coordinates and energies of all the species involved in the present work.
Author Contributions: S.H. and W.N. formulated the project. Y.T. optimized all compounds 1a-6e after conformer search, applying the Monte-Carlo method. S.H., Y.T. and W.N. calculated the QTAIM functions and evaluated the QTAIM-DFA parameters and analyzed the data. S.H. and W.N. wrote the paper, while Y.T. organized the data to assist the writing.

Conflicts of Interest:
The authors declare no conflict of interest.