New Insights on Glutathione’s Supramolecular Arrangement and Its In Silico Analysis as an Angiotensin-Converting Enzyme Inhibitor

Angiotensin-converting enzyme (ACE) inhibitors are one of the most active classes for cardiovascular diseases and hypertension treatment. In this regard, developing active and non-toxic ACE inhibitors is still a continuous challenge. Furthermore, the literature survey shows that oxidative stress plays a significant role in the development of hypertension. Herein, glutathione’s molecular structure and supramolecular arrangements are evaluated as a potential ACE inhibitor. The tripeptide molecular modeling by density functional theory, the electronic structure by the frontier molecular orbitals, and the molecular electrostatic potential map to understand the biochemical processes inside the cell were analyzed. The supramolecular arrangements were studied by Hirshfeld surfaces, quantum theory of atoms in molecules, and natural bond orbital analyses. They showed distinct patterns of intermolecular interactions in each polymorph, as well as distinct stabilizations of these. Additionally, the molecular docking study presented the interactions between the active site residues of the ACE and glutathione via seven hydrogen bonds. The pharmacophore design indicated that the hydrogen bond acceptors are necessary for the interaction of this ligand with the binding site. The results provide useful information for the development of GSH analogs with higher ACE inhibitor activity.


Introduction
Cardiovascular disease (CD) is still the leading cause of mortality worldwide [1]. It is one of the costliest diseases for governments and healthcare systems (up to USD 320 billion/year) [2]. In 2019, according to a report from the World Health Organization (WHO), 17.9 million people died from CD, including clinical disorders of the heart and blood vessels [3]. There are several risk factors for CD, such as high blood pressure, diabetes mellitus, smoking, dyslipidemia, and being overweight/obesity [4,5]. Regarding increased blood pressure, its physiopathology is multifactorial and based on salt intake, obesity/insulin resistance, the influence of the sympathetic nervous system, and the renin-angiotensin-aldosterone system [6]. Among these factors, the renin-angiotensinaldosterone system exerts a controlling function on blood pressure and hypertension via the synthesis of angiotensin II (a potent vasoconstrictor) [7] catalyzed by the angiotensinconverting enzyme (ACE, EC 3.4.15.1) [8]. ACE inhibitors comprise the first line of medicines employed in hypertension therapy, heart failure, myocardial infarction, and diabetic nephropathy [9].

Solid State Analysis
The crystal structure of the GSH is a ubiquitous thiol-containing tripeptide (L-γglutamyl-L-cysteinylglycine) that exists in its zwitterionic form. This structure is also reflected in the IR spectrum of the compounds recorded in KBr. Above the 3000 cm −1 regions, several associated OH, NH, NH 2 , and NH 3 + (3346, 3249, and 3126) bands can be seen. A band at 2524 cm −1 can be associated with the cysteine thiol (SH) group. The 1713 cm −1 band corresponds to the νC=O band of the protonated carboxyl groups, while the strong bands at 1537 and 1385 cm −1 are related to the asymmetric and symmetric stretches, respectively, of the deprotonated CO 2 group [29]. The low molecular weight sulfur (thiols) present in GSH are easily oxidized and can be regenerated rapidly; these characteristics allow them to play an essential role in cell biology, such as protecting cells via an antioxidant process [30]. The polymorphism has been observed only under ambient pressure conditions, revealing that their polymorphs crystallize in the orthorhombic space group P2 1 2 1 2 1 . The differences in the geometry, crystal data, and structure refinement details for GSHA and GSHB are summarized in Table 1. The ORTEP diagrams of the GSHA and GSHB can be found in Figure 1, as well as the overlap of these polymorphs. The results obtained for the GSH conformer molecular geometries were compared to the crystallographic model through the mean absolute deviation percentage formula, where χ DFT and χ XRD indicate the theoretical and experimental bond length or angle, respectively. The exchange and correlation functional M06-2X captures medium-range electronic correlation and can satisfactorily describe scattering interactions such as noncovalent interactions [31,32]. Other functionals were tested, such as the hybrid functional B3LYP, the highly parameterized empirical functional M06-HF [31], and the double-hybrid functional B2-PLYP [33]. However, the MAPD values showed that the functional M06-2X can satisfactorily describe scattering interactions, such as non-covalent interactions [31]. The graph presented in Figure 2 shows the results of the MAPD values obtained for the different levels of theory for GSH.  The results obtained for the GSH conformer molecular geometries were compared to the crystallographic model through the mean absolute deviation percentage formula, where and indicate the theoretical and experimental bond length or angle, respectively. The exchange and correlation functional M06-2X captures medium-range electronic correlation and can satisfactorily describe scattering interactions such as noncovalent interactions [31,32]. Other functionals were tested, such as the hybrid functional B3LYP, the highly parameterized empirical functional M06-HF [31], and the doublehybrid functional B2-PLYP [33]. However, the MAPD values showed that the functional M06-2X can satisfactorily describe scattering interactions, such as non-covalent interactions [31]. The graph presented in Figure 2 shows the results of the MAPD values obtained for the different levels of theory for GSH.  The MAPD values obtained for bond lengths and angles, at the M06-2X/6-311++G(d,p) level of theory, were 1.994 and 1.530%, respectively, where the Pearson correlation coefficients were 0.9790 and 0.9816. Figure 3 presents the graphs comparing  The MAPD values obtained for bond lengths and angles, at the M06-2X/6-311++G(d,p) level of theory, were 1.994 and 1.530%, respectively, where the Pearson correlation coefficients were 0.9790 and 0.9816. Figure 3 presents the graphs comparing the geometric properties obtained experimentally and theoretically. Regarding the lengths, stretching of the N 2 -C 4 , N 3 -C 8 , O 1 -C 5 , O 6 -C 10 , and C 10 -C 9 and the compression of the N 1 -C 1 and O 2 -C 5 bonds were observed. In the GSH crystal, the molecules are found in the zwitterionic form, observed in the glutamic acid portion. In this structure, the O 1 -C 5   The zwitterionic form of GSH is predominant in the crystalline state. However, thermodynamic calculations (Table 2) showed that, when it is kept isolated, the lowest energy state for the molecule is the neutral form, and 134.03 kcal/mol is more stable than the ionized form. In addition, all thermodynamic parameters (internal energy, enthalpy, free energy, entropy, etc.) calculated for both forms resulted in values lower than their ionized state, except in the case of entropy, which is 18.47 cal/mol•K higher for the neutral form. A relaxed scan calculation was carried out to verify the change in the total energy of the molecular system in which the proton starts from the N1 atom, bound at 1.0 Å, towards the O1 atom. From the scan plot shown in Figure 4, it is possible to verify that the system's total energy decays as H + approaches the conjugated carbonyl base. The transit of the proton from the primary amino group to the conjugate base also compresses the N 1 -C 1 bond by 2.8%, leaving the electron pair of the N atom free to resonate along the bond. The calculations also showed that the N 2 -C 4 and N 3 -C 8 peptide bonds are stretched by 3.2 and 3.5% in the free form of the molecule. Finally, the O 6 -C 10 bond of the glycine portion is stretched by 2.6% in the isolated form of the molecule. This effect is explained by the strong O 6 -H· · · O 3 interaction in the crystalline state, causing the proton to approach the O 6 atom, stretching O 6 -C 10 .
The zwitterionic form of GSH is predominant in the crystalline state. However, thermodynamic calculations (Table 2) showed that, when it is kept isolated, the lowest energy state for the molecule is the neutral form, and 134.03 kcal/mol is more stable than the ionized form. In addition, all thermodynamic parameters (internal energy, enthalpy, free energy, entropy, etc.) calculated for both forms resulted in values lower than their ionized state, except in the case of entropy, which is 18.47 cal/mol·K higher for the neutral form. A relaxed scan calculation was carried out to verify the change in the total energy of the molecular system in which the proton starts from the N 1 atom, bound at 1.0 Å, towards the O 1 atom. From the scan plot shown in Figure 4, it is possible to verify that the system's total energy decays as H + approaches the conjugated carbonyl base. - * Di f f = ε n − ε z , where ε n is the thermochemical property of the neutral form and ε z is the ionized form. The zwitterionic form of GSH is predominant in the crystalline state. However, thermodynamic calculations (Table 2) showed that, when it is kept isolated, the lowest energy state for the molecule is the neutral form, and 134.03 kcal/mol is more stable than the ionized form. In addition, all thermodynamic parameters (internal energy, enthalpy, free energy, entropy, etc.) calculated for both forms resulted in values lower than their ionized state, except in the case of entropy, which is 18.47 cal/mol•K higher for the neutral form. A relaxed scan calculation was carried out to verify the change in the total energy of the molecular system in which the proton starts from the N1 atom, bound at 1.0 Å, towards the O1 atom. From the scan plot shown in Figure 4, it is possible to verify that the system's total energy decays as H + approaches the conjugated carbonyl base.  Deviations greater than 2.0% also occurred in the bond angles. In the carboxyl group of the glutamic acid portion, a 4.3% increase in the O 2 -C 5 -C 1 angle and a 2.3% decrease in the O 1 -C 5 -O 2 angle were observed; the angles formed by C 1 were also increased by the proton transition in the molecular structure (3.2% in N 1 -C 1 -C 2 and 3.4% in N 1 -C 1 -C 5 ). On the other hand, the C 5 -C 1 -C 2 angle decreased by 5.2%. However, the distinction between the molecular structures of GSH in both polymorphs lies in the torsion of the tripeptide's carboxyl group of the glycine portion. The overlap of the peptide bond, O 4 =C 8 -N 3 -C 9 , showed that the dihedral angles are −6.93 • and 1.23 • , respectively, in GSHA and GSHB, and the -COOH groups are at 73.63 • and 77 • .51 • out of the plane of the peptide bond in the respective polymorphs. Figure 5 shows the overlapping of the molecular structures of GSHA and GSHB, starting from this dihedral, where the RMS was 0.0292 Å. Torsions in the other portions of the molecule were also observed. Namely, the dihedral angle of C 8 -C 6  showed that the dihedral angles are −6.93° and 1.23°, respectively, in GSHA and GSHB, and the -COOH groups are at 73.63° and 77°.51° out of the plane of the peptide bond in the respective polymorphs. Figure 5 shows the overlapping of the molecular structures of GSHA and GSHB, starting from this dihedral, where the RMS was 0.0292 Å. Torsions in the other portions of the molecule were also observed. Namely, the dihedral angle of C8-C6-C7-S (cysteine portion) is −66.29° and −58.05°, respectively, while the dihedral angle of the glutamic acid portion is 145.94° and 112.44°. Figure 5. Torsion of the carboxyl group, -COOH, of the glycine portion in the glutathione molecule, and the overlapping of their molecular structures from the peptide bond that joins the glycine chain to the cysteine chain.

Molecular Modeling Analysis
The frontier molecular orbitals, HOMO and LUMO, are represented in Figure 6, and their values are shown in Table 3. These orbitals resulted in a very high energy gap (242.6 kcal/mol), calculated by GAP = − , and indicating the high kinetic stability of the GSH molecule. Its antioxidant power is related to protection against reactive oxygen species and electrophilic species produced in cellular oxidative processes; therefore, the values found for the energies of the frontier orbitals agree with this statement, indicating that the GSH molecule is a reductant ( < 0), but not an oxidant ( > 0). HOMO has a lone pair and pure p orbital characteristics, whose occupancy is 1.85558e. LUMO, on

Molecular Modeling Analysis
The frontier molecular orbitals, HOMO and LUMO, are represented in Figure 6, and their values are shown in Table 3. These orbitals resulted in a very high energy gap (242.6 kcal/mol), calculated by GAP = E LUMO − E HOMO , and indicating the high kinetic stability of the GSH molecule. Its antioxidant power is related to protection against reactive oxygen species and electrophilic species produced in cellular oxidative processes; therefore, the values found for the energies of the frontier orbitals agree with this statement, indicating that the GSH molecule is a reductant (E HOMO < 0), but not an oxidant (E LUMO > 0). HOMO has a lone pair and pure p orbital characteristics, whose occupancy is 1.85558e. LUMO, on the other hand, is an σ * antibonding orbital located on the longitudinal axis of the C 10 -O 5 bond, formed by the contribution of 65.02% of the C 10 atom sp 1.86 and 34.98% of the O 5 atom sp 1.37 ; it has an occupancy of 0.01811e.  The sulfhydryl group (-SH) of the cysteine portion is highly polarizable, characterizing it as a good nucleophile [34,35]. The bonding orbital of the S-H bond has an occupancy of 1.  The sulfhydryl group (-SH) of the cysteine portion is highly polarizable, characterizing it as a good nucleophile [34,35]. The σ bonding orbital of the S-H bond has an occupancy of 1.99027e and is formed by the contribution of 57.66% of the sp 5.61 d 0.04 orbital of S and 42.34% of the s orbital of H (Figure 7a). The lone pair form the η 1 bonding orbitals (sp 0,47 ), with an occupancy of 1.98974e (Figure 7b), and η 2 bonding orbital (p), with an occupancy of 1.95353e (Figure 7c). The calculated energies of these orbitals were −406.650, −460.529, and −213.606 kcal/mol, explaining the process of oxidation of GSH to glutathione disulfide through the activity of the glutathione oxidase enzyme, glutathione peroxidase, and glutathione reductase.  The sulfhydryl group (-SH) of the cysteine portion is highly polarizable, characterizing it as a good nucleophile [34,35]. The bonding orbital of the S-H bond has an occupancy of 1.99027e and is formed by the contribution of 57.66% of the sp 5.61 d 0.04 orbital of S and 42.34% of the s orbital of H (Figure 7a). The lone pair form the bonding orbitals (sp 0,47 ), with an occupancy of 1.98974e (Figure 7b), and bonding orbital (p), with an occupancy of 1.95353e (Figure 7c). The calculated energies of these orbitals were −406.650, −460.529, and −213.606 kcal/mol, explaining the process of oxidation of GSH to glutathione disulfide through the activity of the glutathione oxidase enzyme, glutathione peroxidase, and glutathione reductase. From the energies of the HOMO and LUMO, the chemical reactivity descriptors of the GSH were determined as the potential chemical, chemical hardness From the energies of the HOMO and LUMO, the chemical reactivity descriptors of the GSH were determined as the potential chemical, chemical hardness and the global electrophilicity index, where E is the energy of the system, N is the number of electrons, υ is the external potential generated by nuclei, I ≈ −E HOMO is the ionization potential, and A ≈ −E LUMO is the electron affinity. The high values found for the energy gap, as well as η, indicate that the GSH molecule is kinetically stable, having a low electron affinity and high electron transfer power during chemical processes. Compared to other organic compounds, the results show that the GSH molecule has a nucleophilic character (ω < 0.93 a.u.) [36]. Oxygen atoms of the carboxyl groups and the carbonyl groups of GSH have a high charge density, showing the behavior of a Lewis base. This can be seen by the red color on the electronic isodensity surface of the MEP map represented in Figure 8. On the other hand, regions of lower charge density appear in blue and show the behavior of a Lewis acid.
Compared to other organic compounds, the results show that the GSH molecule has a nucleophilic character ( < 0.93 a.u.) [36]. Oxygen atoms of the carboxyl groups and the carbonyl groups of GSH have a high charge density, showing the behavior of a Lewis base. This can be seen by the red color on the electronic isodensity surface of the MEP map represented in Figure 8. On the other hand, regions of lower charge density appear in blue and show the behavior of a Lewis acid. To determine the local electrophilicity, the Fukui function [37,38] was used to predict the regions of nucleophilic, To determine the local electrophilicity, the Fukui function [37,38] was used to predict the regions of nucleophilic, electrophilic, and radical attacks, Radical reactions caused by pathological processes, by the administration of oxidizing drugs, or even by physical activities can cause oxidative transformations of phospholipids, proteins, and deoxyribonucleic acid (DNA) in cell membranes. According to Fukui index calculations, the O 1 and O 2 atoms of the carboxyl group and the N 1 atom of the amine group in the glutamic acid portion (Figure 9a) are favorable to the capture of free radicals produced in these processes, which explains the antioxidant character of GSH. Furthermore, the S atom in the cysteine portion can also undergo radical attacks resulting in the oxidized form of GSH, glutathione disulfide (GSSG), by the intervention of the enzymes glutathione oxidase or glutathione peroxidase. The calculations of the Fukui indices also indicated that GSH has active sites favorable to electrophilic attacks, identified by the isosurface regions in Figure 9b. All the oxygen and nitrogen atoms along the GSH molecule favor this type of attack. This can be explained by the mesomeric effect caused by the delocalized electrons of the carboxyl and amide groups in the tripeptide molecule. In addition, C atoms are favorable to electrophilic attacks, resulting from the inductive effect caused by the presence of O and N atoms, which reduces the charge density in the carbon chain skeleton.
ces also indicated that GSH has active sites favorable to electrophilic attacks, identified by the isosurface regions in Figure 9b. All the oxygen and nitrogen atoms along the GSH molecule favor this type of attack. This can be explained by the mesomeric effect caused by the delocalized electrons of the carboxyl and amide groups in the tripeptide molecule. In addition, C atoms are favorable to electrophilic attacks, resulting from the inductive effect caused by the presence of O and N atoms, which reduces the charge density in the carbon chain skeleton.

Supramolecular Arrangement
The molecular topology of the crystals by HS mapped with helps us to understand better the crystallographic forces driving the molecular arrangement and the interactions reported in Table 4. The contacts that are shorter than the sum of the van der Waals radii are represented by the red spots on a predominantly blue surface [39]. The crystal packing of GSHA is stabilized by the bifurcated intermolecular interactions S-H⋯O1 and S-H⋯O2 observed in the cysteine portion, involving the thiol group forming a motif [40] R (4), and, for GSHB, S-H⋯O1 was observed forming a motif C (11) and the

Supramolecular Arrangement
The molecular topology of the crystals by HS mapped with d norm helps us to understand better the crystallographic forces driving the molecular arrangement and the interactions reported in Table 4. The contacts that are shorter than the sum of the van der Waals radii are represented by the red spots on a predominantly blue surface [39]. The crystal packing of GSHA is stabilized by the bifurcated intermolecular interactions S-H· · · O 1 and S-H· · · O 2 observed in the cysteine portion, involving the thiol group forming a motif [40] R 2 1 (4), and, for GSHB, S-H· · · O 1 was observed forming a motif C 1 1 (11) and the intermolecular interactions N 1 -H· · · O 2 [C 1 1 (5)] and N 2 -H· · · O 2 [C 1 1 (8)]; these contacts form bifurcates in each compound. The red convex area above the glutamic acid correlates with the hydrogen bond interactions N 1 -H· · · O 4 [C 1 1 (5)] and N 1 -H· · · O 1 [C 1 1 (5)]. Furthermore, the interactions O 6 -H· · · O 3 [C 1 1 (10)] and N 3 -H· · · O 5 [C 1 1 (5)] involve part of the glycine. The presence of weak non-classical H-bond interactions is also evident in the GSHA, which is C 2 -H· · · O 4 [C 1 1 (8)] and, for GSHB, C 3 -H· · · O 2 [C 1 1 (6)]. The fingerprint plots ( Figure 10) allow us to analyze the differences in the intermolecular patterns of the contacts and quantitatively evaluate the contributions among atoms [41]. The contacts involving H· · · O accumulated a percentage of 50% and are viewed as a distinct pair of spikes, evidence that the H-bonds are dominant in the crystalline environment [42]. The H· · · H contacts, close to 34%, with a decrease of 0.8% for GSHB, contribute to the overall crystal packing. In contrast, the proportion of S· · · H increased from 6.3% for GSHA to 9.6% for GSHB. In addition, weak C· · · H contacts decreased from 3% in GSHA to 2% in GSHB, indicating that the GSHA forms more H-bonds than GSHB.
The topological parameters for the intermolecular interactions in the molecular arrangements of the GSH polymorphs are presented in Table 4, and the molecular graphs are represented in Figure 11. In QTAIM, the observable properties of the molecular system are contained in electron density ρ(r), in which the Laplacian electron density, ∇ 2 ρ, determines the depletions or peaks of the electron charge concentration between nuclear attractors, indicating the location of bond critical points (BCPs) [43][44][45]. In other words, ∇ 2 ρ corresponds to the concentration of the electronic charge in the intranuclear region of two attractors. If the electron density is accumulated in the intranuclear region, its value is negative in the BCP, and the interaction is shared so that the attractors are covalently attached. On the other hand, if ∇ 2 ρ > 0, the electron density is concentrated in the attractors, a closed-shell interaction in which the attractors are linked by weak interactions.
The low values of the electron density (ρ < 0.1) and the positive Laplacian found on the BCP indicate closed-shell interactions [43][44][45][46]. Furthermore, when the values for the ratio |v|/G < 1.0 and the total energies are very low, h(r) ≈ 0, this indicates the intermolecular interactions in the supramolecular arrangement of GSH are of low intensity, configuring hydrogen bonding. In the O 6 -H· · · O 3 and N 1 -H· · · O 2 interactions, the values found for the ratio |v|/G = 1.0, with h < 0; however, these values do not indicate that the hydrogen bonds in these interactions have any covalent character, as the values of h(r) are minimal. ence of weak non-classical H-bond interactions is also evident in the GSHA, which is C2-H⋯O4 [C (8)] and, for GSHB, C3-H⋯O2 [C (6)].
The fingerprint plots ( Figure 10) allow us to analyze the differences in the intermolecular patterns of the contacts and quantitatively evaluate the contributions among atoms [41]. The contacts involving H⋯O accumulated a percentage of 50% and are viewed as a distinct pair of spikes, evidence that the H-bonds are dominant in the crystalline environment [42]. The H⋯H contacts, close to 34%, with a decrease of 0.8% for GSHB, contribute to the overall crystal packing. In contrast, the proportion of S⋯H increased from 6.3% for GSHA to 9.6% for GSHB. In addition, weak C⋯H contacts decreased from 3% in GSHA to 2% in GSHB, indicating that the GSHA forms more H-bonds than GSHB.
(a) (b) The topological parameters for the intermolecular interactions in the molecular arrangements of the GSH polymorphs are presented in Table 4, and the molecular graphs are represented in Figure 11. In QTAIM, the observable properties of the molecular system are contained in electron density ( ), in which the Laplacian electron density, ∇ , determines the depletions or peaks of the electron charge concentration between nuclear attractors, indicating the location of bond critical points (BCPs) [43][44][45]. In other words, ∇ corresponds to the concentration of the electronic charge in the intranuclear region of two attractors. If the electron density is accumulated in the intranuclear region, its value is negative in the BCP, and the interaction is shared so that the attractors are covalently attached. On the other hand, if ∇ > 0, the electron density is concentrated in the attractors, a closed-shell interaction in which the attractors are linked by weak interactions. The low values of the electron density ( < 0.1) and the positive Laplacian found on the BCP indicate closed-shell interactions [43][44][45][46]. Furthermore, when the values for the ratio | |/ < 1.0 and the total energies are very low, ℎ( ) ≈ 0, this indicates the intermolecular interactions in the supramolecular arrangement of GSH are of low intensity, configuring hydrogen bonding. In the O6-H⋯O3 and N1-H⋯O2 interactions, the values found for the ratio   Symmetry codes:  Closed-shell interactions are weak compared to shared interactions. The analysis of the NBO calculations showed that the hydrogen bonds in the supramolecular arrangements of GSH are stabilized by the hyperconjugation of lone pairs (Lewis type) with σ * antibonding orbitals (non-Lewis type) of the interacting region of the neighboring molecule at higher energies. The most significant results of the NBO calculations are found in Tables S1 and S2 (Supplementary Materials). They show that the higher hyperconjugation energies stabilize the molecules in the GSHB polymorph more than in the GSHA polymorph. As an example, we cite the interactions that are commonly found in both polymorphs with high E (2) energies. The O 6 -H· · · O 3 interaction is stabilized by the hyperconjugation of the lone pairs of the O 3 atom with the σ * antibonding orbitals of the O 6 The hyperconjugation of the η 1 (O 3 ) orbital (occupancy 1.98e, sp 0.7 hybrid) resulted in energy equal to 8.52 kcal/mol in GSHA and 7.80 kcal/mol in GSHB (ratio 1.1:1); on the other hand, the hyperconjugation of the η 2 (O 3 ) orbital (occupancy 1.88, pure p) resulted in energy equal to 9.45 kcal/mol in GSHA and 22.01 kcal/mol in GSHB (ratio 1:2.3). This means that the O 6 -H· · · O 3 interaction is more cohesive in the GSHB polymorph. Two other essential interactions in both polymorphs were N 2 -H· · · O 2 and N 3 -H· · · O 5 . The first is stabilized by the hyperconjugation of the two lone pairs of O 2 , η 1 (occupancy 1.98e, sp 0.6 hybrid) and η 2 (occupancy 1.89e, pure p) with the σ * antibonding orbital of the N 2 -H bond in GSHA, with energies of 2.01 and 2.93 kcal/mol, but only by the lone pair η 2 (occupancy 1.96e, pure p) in GSHB, with an energy of 6.80 kcal/mol. In the N 3 -H· · · O 5 interaction, only the η 2 lone pair of O 5 (occupancy 1.86e, pure p) contributes to its stabilization in GSHA, while the two lone pairs of O 5 , η 1 (occupancy 1.97e, hybrid sp 0.7 ) and η 2 (occupancy 1.87e, pure p) hyperconjugate with the σ * antibonding orbital of the N 3 -H bond in GSHB; however, both energies in this polymorph further stabilize the interaction.

Molecular Docking Analysis
The ACE (EC 3.4.15.1), a Zn metalloproteinase, plays an essential role in cardiovascular function by converting the decapeptide angiotensin I to angiotensin II (vasopressor octapeptide). There are two isoforms of ACE (somatic and testis) that are transcribed from the same gene in a tissue-specific manner [47]. The somatic form consists of two homologous domains (N and C domain), each containing an active site with a conserved zinc-binding motif. Despite the high sequence similarity between the two domains, they differ in substrate and inhibitor specificity and their activation by chloride ions [48]. The C domain seems to be the dominant angiotensin-converting site [47]. Docking calculations indicate that lisinopril had the higher affinity for ACE_C. [48]. Recently, the binding activity of several peptides with ACE inhibitory activity was modeled by in silico approaches, such as the quantitative structure-activity relationship (QSAR) and molecular docking. The methods provide details illuminating the interaction mechanism between the receptor and the ligands [48][49][50][51][52].
In the present work, using redocking analysis, molecular docking was validated for a co-crystallized ligand (lisinopril, a commercial ACE inhibitor) with the ACE. The same model employed in redocking lisinopril was also used for GSH. Seven of the ten poses presented RMSD values less than 2.0 Å. As shown in Figure 12, the residues involved in the interaction and the types of non-bonding interactions involved between the studied complexes were Gln281, Ala353, Ala354, Ala 356, Tyr520, Tyr523, His387, and Phe512 (Conventional hydrogen bond, Unfavorable Acceptor-Acceptor, Pi-Sulfur, van der Waals). Some compounds of the herbal species Cucurbita pelo L. presented potential interaction in silico with the ACE, and the His353 also interacted by a hydrogen bond. Thus, it seems that His353 represents a critical interaction point. The number and the distances of hydrogen bonds play an essential role in the potential biological interaction between the ligand and the binding site. According to Figure 13, the hydrogen bond distances ranged from 1.91 to 2.60 Å [53,54].  Furthermore, the same carboxyl group interacting with Gln281 also appears in the pharmacophore model (Figures 12c and 14). GSH is a well-known endogenous antioxidant [55,56]. This substance is synthesized in the liver and is involved in various organic processes such as antioxidant defense, metabolism, and regulation. Moreover, after drug treatment, GSH reduces oxidative stress in the heart of diabetic rabbits and also has a protective effect in mice from sepsis by inhibiting the inflammation process [57,58]. Thus, GSH may have protective cardiovascular effects associated with its antioxidant properties, which could provide vascular protection. Furthermore, GSH presents a potential impact on antihypertension, and its antioxidant properties can reduce vascular damage [59].  Furthermore, the same carboxyl group interacting with Gln281 also appears in the pharmacophore model (Figures 12c and 14). GSH is a well-known endogenous antioxidant [55,56]. This substance is synthesized in the liver and is involved in various organic processes such as antioxidant defense, metabolism, and regulation. Moreover, after drug treatment, GSH reduces oxidative stress in the heart of diabetic rabbits and also has a protective effect in mice from sepsis by inhibiting the inflammation process [57,58]. Thus, GSH may have protective cardiovascular effects associated with its antioxidant properties, which could provide vascular protection. Furthermore, GSH presents a potential impact on antihypertension, and its antioxidant properties can reduce vascular damage [59]. In addition, when oxidative stress creates vascular damage by promoting cell growth, extracellular matrix protein deposition, endothelial dysfunction, and increased vascular tone, these features contribute to the vascular phenotype in hypertension [60]. In the study performed by Bessa, S. S. and colleagues (2009), it was observed that in hypertensives, both systolic and diastolic pressures were negatively correlated with glutathione S-transferase in Egyptian patients [61]. Beyond oxidative stress, our results suggest that glutathione might interact with some regulatory points in the cardiovascular system, such as the angiotensin-converting enzyme, playing an essential role in blood pressure regulation.

Computational Methods
The crystallographic information files of GSHA and GSHB were obtained from the Cambridge Crystallographic Data Centre (CCDC) [62] under codes 762195 and 1500579, respectively. Theoretical calculations were carried out for both conformations [63] using the highly parameterized empirical exchange-correlation functional M06-2X [31], combined with the triple-basis set 6-311++G(d,p), in gas phase, implemented in the Gauss-ian16 [64] software package. Frontier molecular orbitals (FMOs) [63], the highest occupied molecular orbital (HOMO), and the lowest unoccupied molecular orbital (LUMO) were calculated and the chemical reactivity and kinetic stability of both molecules were obtained. To verify the electronic charge distribution on the GSH molecular surface, the MEP [65] map was obtained through electrostatic potential ( ) [66], In addition, when oxidative stress creates vascular damage by promoting cell growth, extracellular matrix protein deposition, endothelial dysfunction, and increased vascular tone, these features contribute to the vascular phenotype in hypertension [60]. In the study performed by Bessa, S. S. and colleagues (2009), it was observed that in hypertensives, both systolic and diastolic pressures were negatively correlated with glutathione S-transferase in Egyptian patients [61]. Beyond oxidative stress, our results suggest that glutathione might interact with some regulatory points in the cardiovascular system, such as the angiotensin-converting enzyme, playing an essential role in blood pressure regulation.

Computational Methods
The crystallographic information files of GSHA and GSHB were obtained from the Cambridge Crystallographic Data Centre (CCDC) [62] under codes 762195 and 1500579, respectively. Theoretical calculations were carried out for both conformations [63] using the highly parameterized empirical exchange-correlation functional M06-2X [31], combined with the triple-ζ basis set 6-311++G(d,p), in gas phase, implemented in the Gaussian16 [64] software package. Frontier molecular orbitals (FMOs) [63], the highest occupied molecular orbital (HOMO), and the lowest unoccupied molecular orbital (LUMO) were calculated and the chemical reactivity and kinetic stability of both molecules were obtained. To verify the electronic charge distribution on the GSH molecular surface, the MEP [65] map was obtained through electrostatic potential V(r) [66], where Z α is the charge of nuclei α at point r α and ρ(r ) is the charge density at the point r ' . This function has been used for predicting nucleophilic and electrophilic regions and is valuable in studying intermolecular interactions. To predict the local electrophilicity, we used the Fukui function [67] where N is the number of electrons present in the system and the constant term v in the partial derivative is external potential. Fukui indices have been widely used to predict the reactive sites of organic molecules, so the function values are higher in these regions than in regions with no reactivity. All calculations were carried out by DFT, implemented in the Gaussian 16 [64] software package at the M06-2X/6-311++G(d,p) level of theory in the gas phase, and the results were obtained by the GaussView 6 [68] software.

Supramolecular Arrangement
HS analysis using the CrystalExplorer program [69] was employed to study the intermolecular interactions [39]. The normalized distance (d norm ) is constructed by the distance from the surface to the nearest nucleus inside the surface (d i ) and outside the surface (d e ), the distance from the surface to the nearest inner atom, and the van der Waals radii of the internal and external atoms (r vdW i and r vdW e ) [39], In addition, an analysis of intermolecular contacts and their contributions to the packaging of crystals, based on the combination of d e vs. d i plots, generate fingerprints, which summarize the percentage contribution to the nature and type of intermolecular interaction present in the molecule [41]. The topological parameters of GSH molecular systems were obtained by the QTAIM [43] by Multiwfn software [70], and the stabilities of the interaction were verified by NBO [71] calculations; a second-order perturbation theory formula provides the hyperconjugation energies: where < σ F σ > 2 or F 2 ij is the Fock matrix element between the natural bond orbitals i and j; ε σ * is the energy of the antibonding orbital σ * and ε σ is the energy of the bonding orbital σ; and n σ stands for the population occupation of the σ donor orbital. Theoretical calculations were carried out at the M06-2X/6-311G++(d,p) level of theory, in which the atom positions were fixed in their crystallographic positions.

Molecular Docking Analysis
ACE (PDB ID: 1O86) was downloaded from Protein Data Bank (https://www.rcsb.org; accessed on 12 September 2022) for molecular docking simulations. The enzyme was prepared by adding hydrogen atoms, and water molecules were removed. Next, the GSH molecule was saved in mol2 format. After that, the binding site of the macromolecule was selected, and the grid box was determined (X, Y, and Z coordinates). The downloaded file was saved in PDB format before being subjected to GOLD Suite 5.7.0 (Mark Thompson and Planaria Software LLC) [72,73] to locate the binding site and determine the binding affinity. Biovia Discovery Studio 3.5 software was employed to obtain the 2D interaction figures, and PyMOL Molecular Graphics System 2.0 was used to obtain the 3D interaction images. Redocking with the co-crystallized ligand (lisinopril) was carried out to validate the models obtained. The CHEMPLP score function was used.

Pharmacophore Design
Binding DB (https://www.bindingdb.org/bind/index.jsp; accessed on 7 September 2022) was used to identify the compounds with antagonist activity against the ACE. The five molecules with the lowest half-maximal inhibitory concentration (IC 50 ) were selected to obtain pharmacophore models. PharmaGist web server [73,74] (https://bioinfo3d. cs.tau.ac.il/PharmaGist/; accessed on 10 September 2022) was employed to obtain the 3D pharmacophore model from the set of ACE antagonists and GSH. The algorithm determines potential pharmacophores and selects the highest-scoring ones. The suggested pharmacophore hits were determined after the spatial alignment of the GSH with the five EGFR antagonists. Pharmacophore analysis was used to search the possible critical features of ACE antagonists responsible for the interaction of the binding site. Hydrogen bond acceptors, donors, hydrophobic groups, and aromatic rings were selected to generate the pharmacophore models. A minimum of three to a maximum of six elements was considered. The following weights were used as default in the program (3.0 for aromatic rings, 1.0 for charges, 1.5 for hydrogen bond donors or acceptors, and 0.3 for hydrophobic groups). The images were generated in Biovia Discovery Studio 3.5

Conclusions
A comparison of the crystalline structures of GSH showed that the torsion of the carboxyl group in the glycine portion did not change the crystalline characteristics of the polymorphs, so they are in the same space group. Furthermore, the theoretical calculations showed that, when kept isolated, the molecule tends to remain in the neutral form, while in the crystalline environment, the zwitterionic form is predominant. The calculations also proved the antioxidant power of GSH by the frontier molecular orbitals and showed the regions where electrophilic and radical attacks occur in chemical processes. In supramolecular arrangements, some variations were observed in the contacts that form the intermolecular interactions, so they are more stabilized in GSHB compared to GSHA. In addition, the molecular docking analysis showed interactions between the active site residues of the ACE and glutathione through seven hydrogens. Moreover, these hydrogen bond acceptors are essential for binding the GSH with the ACE. The observed results suggest that glutathione analogs would be attractive dual-acting (antioxidant plus ACE inhibitor) antihypertensive agents and provide useful information for developing GSH analogs with higher ACE inhibitor activity.