Identification of Some Glutamic Acid Derivatives with Biological Potential by Computational Methods

Glutamic acid is a non-essential amino acid involved in multiple metabolic pathways. Of high importance is its relationship with glutamine, an essential fuel for cancer cell development. Compounds that can modify glutamine or glutamic acid behaviour in cancer cells have resulted in attractive anticancer therapeutic alternatives. Based on this idea, we theoretically formulated 123 glutamic acid derivatives using Biovia Draw. Suitable candidates for our research were selected among them. For this, online platforms and programs were used to describe specific properties and their behaviour in the human organism. Nine compounds proved to have suitable or easy to optimise properties. The selected compounds showed cytotoxicity against breast adenocarcinoma, lung cancer cell lines, colon carcinoma, and T cells from acute leukaemia. Compound 2Ba5 exhibited the lowest toxicity, and derivative 4Db6 exhibited the most intense bioactivity. Molecular docking studies were also performed. The binding site of the 4Db6 compound in the glutamine synthetase structure was determined, with the D subunit and cluster 1 being the most promising. In conclusion, glutamic acid is an amino acid that can be manipulated very easily. Therefore, molecules derived from its structure have great potential to become innovative drugs, and further research on these will be conducted.


Introduction
Throughout history, cancer has been a major health problem. It has been shown that there is a positive correlation between cancer incidence and age [1][2][3]. The individual risk of cancer is also influenced by family history, genetic susceptibility or behaviour, and exposure to carcinogenic factors [4]. Furthermore, the Krebs cycle and amino acids are proven to significantly affect cancer metabolism. Thus, interfering with amino acid metabolic pathways is an active area of study in cancer metabolism [5].
Additional pathways involving amino acid transport suggest effective therapies. Tumour cells achieve high intracellular concentrations of glutamine primarily through the upregulation of glutamine transporters, including ASCT2 (alanine, serine, cysteine transporter 2 or SLC1A5) [5]. Pharmacological blockade of SLC1A5 can be a successful alternative in some types of cancer. V-9302, an SLC1A5 antagonist (Figure 3), elicited a marked anti-tumour response in preclinical tumour models [10,11]. It has blocked glutamine uptake in a broad spectrum of solid tumours (such as colorectal cancer cell lines) and several xenograft tumour models. This blocked glutamine uptake resulted in a profound alteration of tumour cell growth and survival [9,21]. It has been observed that V-9302 was more productive in inducing triple-negative breast cancer cell death in several human and mouse cell culture models [16]. The combination of CB-839 and V-9302 was also successful because of the dual inhibition of glutamine metabolism, resulting in a decrease in GSH levels and a lethal increase in the levels of free radicals. This resulted in severe DNA damage, especially in liver cancer cells [13]. Additional pathways involving amino acid transport suggest effective therapies. Tumour cells achieve high intracellular concentrations of glutamine primarily through the upregulation of glutamine transporters, including ASCT2 (alanine, serine, cysteine transporter 2 or SLC1A5) [5]. Pharmacological blockade of SLC1A5 can be a successful alternative in some types of cancer. V-9302, an SLC1A5 antagonist (Figure 3), elicited a marked anti-tumour response in preclinical tumour models [10,11]. It has blocked glutamine uptake in a broad spectrum of solid tumours (such as colorectal cancer cell lines) and several xenograft tumour models. This blocked glutamine uptake resulted in a profound alteration of tumour cell growth and survival [9,21]. It has been observed that V-9302 was more productive in inducing triple-negative breast cancer cell death in several human and mouse cell culture models [16]. The combination of CB-839 and V-9302 was also successful because of the dual inhibition of glutamine metabolism, resulting in a decrease in GSH levels and a lethal increase in the levels of free radicals. This resulted in severe DNA damage, especially in liver cancer cells [13]. Another therapeutic strategy could be inhibiting glutamate carboxypeptidase II (GCPII). This enzyme hydrolyses N-acetyl-aspartyl-glutamate (NAAG) to glutamate and N-acetyl aspartate. NAAG is a neurotransmitter in the brain and a glutamate provider to GCPIIpositive cancers if other sources do not produce enough glutamate. Therefore, inhibitors of GCPII can lead to cancer cell growth suppression by reducing glutamate concentrations [7]. Antagonists of metabotropic glutamate receptors are also promising anti-cancer alternatives without significant side effects. Metabotropic glutamate receptors (mGluRs) are G-protein coupled receptors (GPCRs) categorised into three groups based on their signal transduction pathways and pharmacological profiles. They seem to be more attractive therapeutic targets since they are not directly involved in excitotoxicity but intervene in modulating glutamate activity [22,23].  Another therapeutic strategy could be inhibiting glutamate carboxypeptidase II (GCPII). This enzyme hydrolyses N-acetyl-aspartyl-glutamate (NAAG) to glutamate and N-acetyl aspartate. NAAG is a neurotransmitter in the brain and a glutamate provider to GCPII-positive cancers if other sources do not produce enough glutamate. Therefore, inhibitors of GCPII can lead to cancer cell growth suppression by reducing glutamate concentrations [7]. Antagonists of metabotropic glutamate receptors are also promising anti-cancer alternatives without significant side effects. Metabotropic glutamate receptors (mGluRs) are G-protein coupled receptors (GPCRs) categorised into three groups based on their signal transduction pathways and pharmacological profiles. They seem to be more attractive therapeutic targets since they are not directly involved in excitotoxicity but intervene in modulating glutamate activity [22,23].
This article aims to identify new structural analogues of glutamic acid as potential candidates for anti-cancer therapy by computational methods. Several stages were followed: (1) analysis of recently published scientific data regarding the role of glutamate and its derivatives in the development of tumour cells; (2) identification of some new molecules with biological potential, starting with the structure of glutamic acid and the creation of a compound library; (3) conjugation of molecules of natural origin with glutamic acid residues to reduce glutamic acid toxicity and/or potentiate the anti-cancer effect; (4) selection of compounds with biological action and minimal toxicity according to the structural, physicochemical, pharmacokinetic, and pharmaco-toxicological properties determined by in silico methods; (5) evaluation of anti-tumour potential of selected molecules and the identification of possible mechanisms of action; (6) molecular dynamics simulation and molecular docking study to identify the binding site of a ligand molecule (with biological potential) on a known target.

Results and Discussion
The designed glutamic acid derivatives were classified by classes, groups, and subgroups (Table 1). Each one of the compounds received an ID code composed of the following elements: first digit-class; capital letter-group; small letter-subgroup; last digit-the compound's number in the subgroup; small letter at the end (if applicable)-a derivative of the lead-compound. The online software and test parameters that were used to obtain and characterise the compounds are mentioned in Table S1 (Supplementary This article aims to identify new structural analogues of glutamic acid as potential candidates for anti-cancer therapy by computational methods. Several stages were followed: (1) analysis of recently published scientific data regarding the role of glutamate and its derivatives in the development of tumour cells; (2) identification of some new molecules with biological potential, starting with the structure of glutamic acid and the creation of a compound library; (3) conjugation of molecules of natural origin with glutamic acid residues to reduce glutamic acid toxicity and/or potentiate the anti-cancer effect; (4) selection of compounds with biological action and minimal toxicity according to the structural, physicochemical, pharmacokinetic, and pharmaco-toxicological properties determined by in silico methods; (5) evaluation of anti-tumour potential of selected molecules and the identification of possible mechanisms of action; (6) molecular dynamics simulation and molecular docking study to identify the binding site of a ligand molecule (with biological potential) on a known target.

Results and Discussion
The designed glutamic acid derivatives were classified by classes, groups, and subgroups (Table 1). Each one of the compounds received an ID code composed of the following elements: first digit-class; capital letter-group; small letter-subgroup; last digit-the compound's number in the subgroup; small letter at the end (if applicable)-a derivative of the lead-compound. The online software and test parameters that were used to obtain and characterise the compounds are mentioned in Table S1 (Supplementary Materials). The structures of all obtained compounds and their computational descriptors are given in Table S2 (Supplementary Materials).
Pharmacokinetic properties were evaluated for each compound in terms of permeability (gastrointestinal absorption, blood-brain barrier permeability) and interactions with P-gp. In addition, we assessed the enzyme inhibitory effect on some isoforms of cytochrome P450 (Table S9; Supplementary Materials). Based on the previously calculated properties, we evaluated whether these compounds meet the "drug-likeness" criteria according to the Lipinski, Ghose, Veber, Egan, and Muegge rules. The number of rules violated by each molecule is shown in Table S10 (Supplementary Materials), along with the bioavailability score, drug-likeness score, lead-likeness score, and synthetic accessibility score.

The Elimination of Reactive and Toxic Compounds
The elimination of reactive and toxic compounds was carried out in several steps, as follows:

•
Step 1. In the first stage, compounds belonging to at least two toxicity classes are eliminated, as the risk of them causing severe adverse reactions is high.

•
Step 2. This step involves the removal of compounds that do not follow Lipinski and Veber's rules, and which have a CNS MPO score less than 4, as well as compounds with low solubility and/or an inhibitory effect on cytochrome P450 and/or gp-P enzymes.

•
Step 3. Compounds with medium toxicity, which fall into Class III (Cramer rules) and are positive for at least one toxicity criterion, are eliminated if the overall drug-likeness score does not exceed 0.90.

•
Step 4. Compounds that have violated all Ghose's rule criteria (four out of four) and belong to Cramer class III or II or overlap with the violation of at least one Muegge rule are eliminated.

•
Step 5. Compounds that have violated at least three Ghose criteria and at least two Muegge rules and belong to Cramer class III are eliminated.

•
Step 6. Removal of Cramer Class III compounds that violate at least one Ghose and Muegge rule, having an SA score below 2.

•
Step 7. Elimination of Class III Cramer compounds that violate at least one Ghose and Muegge rule, regardless of the SA score achieved.

•
Step 8. Removal of compounds that violate at least one Ghose and Muegge rule with a low GI absorption value.

•
Step 9. Compounds that violate at least one Ghose and Muegge rule with an SA score below 4, regardless of Cramer toxicity class, are eliminated.

•
Step 10. Elimination of Cramer Class III compounds that violate at least two Muegge criteria and have an SA score below 3 and/or overall drug-likeness score below 0.5.
Only nine compounds proved to have suitable properties or properties that can be easily optimised, representing 7.3% of the total. These selected compounds are presented in Table 3, along with their geometrical and isomer-conformation properties. Table 3. Structures of the nine "lead" compounds and their geometrical and isomer-conformation properties [105].

No. ID Code
Chemical Structure

1Aa7
Molecules 2023, 28, x FOR PEER REVIEW 9 of 29 Table 3. Structures of the nine "lead" compounds and their geometrical and isomer-conformation properties [105].

Characterisation of the "Lead" Compounds
The "lead" compounds were characterised by chemical structure, geometric isomers, isomerism, and conformations using the MarvinSketch platform [105] (Table 3). The platform automatically generated the conformations, and their number was limited to ten. The energy was calculated using force field methods, and the conformer with the lowest energy, i.e., having the highest stability, was chosen.
The main pathways of metabolism, bioactivity, action on cancer cells, mechanisms of action and possible adverse effects, and acute toxicity in rodents were further evaluated by in silico methods. For this, we used Toxtree [22] to assess the metabolism of the nine compounds (primary, secondary, tertiary, and quaternary sites of metabolism) and also SmartCyp and SOMP to determine the most reactive atom (involved in interactions with CYP3A4, CYP2D6, and CYP2C9) ( Tables 4 and 5). The algorithm used by the Smartcyp online platform requires a reactivity descriptor (E) and an accessibility descriptor (A). "E" estimates the energy required for a CYP to react at this position, and "A" is the relative topological distance of an atom from the centre of the molecule. The score is calculated for each atom according to the equation Score = E − 8*A − 0.04*SASA (where SASA is the solvent-accessible surface area). A lower score corresponds to an increased probability of being a site of metabolism [131].

Characterisation of the "Lead" Compounds
The "lead" compounds were characterised by chemical structure, geometric isomers, isomerism, and conformations using the MarvinSketch platform [105] (Table 3). The platform automatically generated the conformations, and their number was limited to ten. The energy was calculated using force field methods, and the conformer with the lowest energy, i.e., having the highest stability, was chosen.
The main pathways of metabolism, bioactivity, action on cancer cells, mechanisms of action and possible adverse effects, and acute toxicity in rodents were further evaluated by in silico methods. For this, we used Toxtree [22] to assess the metabolism of the nine compounds (primary, secondary, tertiary, and quaternary sites of metabolism) and also SmartCyp and SOMP to determine the most reactive atom (involved in interactions with CYP3A4, CYP2D6, and CYP2C9) ( Tables 4 and 5). The algorithm used by the Smartcyp online platform requires a reactivity descriptor (E) and an accessibility descriptor (A). "E" estimates the energy required for a CYP to react at this position, and "A" is the relative topological distance of an atom from the centre of the molecule. The score is calculated for each atom according to the equation Score = E − 8*A − 0.04*SASA (where SASA is the solvent-accessible surface area). A lower score corresponds to an increased probability of being a site of metabolism [131].  The bioactivity of the nine selected compounds was characterised using the following parameters: G protein-coupled receptor ligand, ion channel modulator, kinase inhibitor, nuclear receptor ligand, protease inhibitor, and enzyme inhibitor (Table 6). In addition, the most probable molecular targets and their identification data were determined using the SWISSTarget predictor (Table 7) [133].   Regarding the interpretation of the results from Table 6, a larger score value correlates with a higher probability for the particular molecule to be active. More explicitly, if the bioactivity score is more than 0.0, the compound is considered active; if the score is between −0.5 and 0.0, it exhibits moderate activity; if the bioactivity score is less than −0.5, then it is inactive [134].
The anticarcinogenic effect of the nine compounds was assessed using CLC-Pred software [135], predicting the most probable cell lines for which compounds exhibit cytotoxicity (Table 8). Table 8. Anticarcinogenic effect: most probable cell lines for which compounds exhibit cytotoxicity. Probability "to be active" (Pa) > Probability "to be inactive" (Pi) [135,136]. Possible mechanisms of action and adverse/toxic effects, lethal doses (LD50) in acute toxicity determined in rodents (intraperitoneal, intravenous, oral, and subcutaneous administration), and the classification of chemical compounds according to the OECD Project were also determined by in silico methods (Tables 9-11) [135][136][137][138].

No. ID Code
Based on the results of the bioactivity assessment by Molinspiration [134] (Table 6), molecular dynamics and docking studies were performed on compound 4Db6 and the bacterial GS enzyme from Salmonella typhimurium ( Figure S1; Supplementary Materials) [43,66,[139][140][141]. The Protein Data Bank (PDB) code for GS is 1lgr [142,143]. Table 9. Mechanisms of action and adverse/toxic effects (Pa > Pi) [137].  Table 10. Acute toxicity in rodents when administered intraperitoneally, intravenously, orally, and subcutaneously: LD50 in mg/kg [138]. The molecular dynamics simulation study was carried out using the UCSF Chimera 1.15 software [144,145]. Before the actual dynamics simulation, the chemical structure was processed according to the protocol established in the literature: hydrogen atoms were inserted, the protonation status corresponding to glutamic acid was used, and Gasteiger partial charges were assigned. The study was performed in water as solvent (SPCBOX, cube size 3 Å) with a density of 1024 g/cm 3 to simulate physiological conditions. In the neutralisation phase, we added Na/Cl counterions. The next step was the minimisation phase, whereby the system's energy would tend towards 0. Table 11. Acute toxicity in rodents. Classification of Chemicals according to the OECD Project [138]. In the equilibration phase, the temperature was set to 310 K (36.85 0 C, approximately physiological temperature) with a gradient of 10 K/ps. In the production phase, the following settings were made: Andersen barostat-pressure 1.0132 bar, relaxation time 1.5; Nose thermostat-emperature 310 K, relaxation time: 0.2. The entire simulation time was set to 100 ns. The energy values resulting from the molecular dynamics simulation for compound 4Db6 are included in Table 12. Geometry optimisation was performed following the Gaussian model, and we used the standard topology for non-protein molecules. Most biological processes involve, at the atomic scale, the recognition of one molecule by another. Estimation of such interactions at the molecular level is performed by docking methods [146]. In the molecular docking study, the interaction of the 4Db6 derivative with the GS enzyme was evaluated in comparison with phosphinothricin ((2S)-2-amino-4-(hydroxy-methyl-phosphoryl)butanoic acid), whose PDB code is PPQ [67,69,70,142,147]. Phosphinothricin, a GS inhibitor, shows the closest similarity (86.9%) to compound 4Db6, as scored by SwissSimilarity (Score = 0.869) [148]. The comparison was made to identify the most probable binding site in the enzyme structure [149].

No. ID Code Rat IP LD50 Classification Rat IV LD50 Classification Rat Oral LD50 Classification
The study was conducted using SwissDock [150][151][152], PatchDock [136,153,154], and AutoDockVina 1.1.2 [151,155]. In a study evaluating a crystalline structure of GS inhibited by phosphinothricin, the inhibitor molecule preferentially binds to the enzyme in the D subunit's active site. Phosphinothricin occupies the glutamate pocket and stabilises the Glu327 residue in a position that prevents glutamate from entering the active site [149]. This crystal structure (PDB code: 1FPY) was observed using the Mol* Viewer web app of RCSB PDB [142,156]. The preference for the D subunit was also confirmed by results obtained using the PatchDock app, which estimated the most probable binding site for the 4Db6 compound [136,153,154]. The top 10 best solutions are shown in Table 13. Figure 4 illustrates the first best result generated.   [136,153,154]; viewed with UCSF Chimera 1.15 [144,145].
For PPQ, SwissDock found 257 conformations. The most probable binding site was chosen according to the conformation with the lowest energy, having ΔG = −10.43 kcal/mol and a FullFitness value of −2192.23 kcal/mol [150,152,157]. The FullFitness parameter for a cluster is calculated using the average of 30% of the most favourable energies of its elements to lower the risk of inhibition of the entire cluster by some complexes. This energy is represented by the sum of the system's total energy and a solvation term [158]. For example, for compound 4Db6, SwissDock found 160 conformations. By comparing the PPQ binding site with the sites of the 160 conformations, we consider that clusters 1, 6 and 33 could bind to the same site in a relatively similar way (Table 14). The inhibition constant (Ki) was calculated using the following formula: Ki = e^((ΔG × 1000)/(R × T)), where e = 2.7182, R = 1.98719 cal/(mol × K) (Regnault constant) and T = 298.15 K = 25 °C [159]. It can be seen that cluster 1 shows the lowest energy according to the ΔG value, but Ki and the maximum FullFitness value belong to complex 33. Visualisation and processing of the results obtained in the molecular docking study ( Figure 5) were performed using UCSF Chimera 1.15 [144,145]. The grid sizes used in SwissDock for cluster 1 are (x, y, z) = (15.5, 15.5, 20.5) with centre coordinates (x, y, z) = (−98, 13.711, −87.161).  [136,153,154]; viewed with UCSF Chimera 1.15 [144,145].
However, the selected derivative does not bind to the active site. Thus, these derivatives will probably not show inhibitory activity towards the enzyme. Molecular docking was performed using SwissDock [134,150,152] and AutoDockVina 1.1.2 [151,155] to increase the accuracy of the study.
For PPQ, SwissDock found 257 conformations. The most probable binding site was chosen according to the conformation with the lowest energy, having ∆G = −10.43 kcal/mol and a FullFitness value of −2192.23 kcal/mol [150,152,157]. The FullFitness parameter for a cluster is calculated using the average of 30% of the most favourable energies of its elements to lower the risk of inhibition of the entire cluster by some complexes. This energy is represented by the sum of the system's total energy and a solvation term [158]. For example, for compound 4Db6, SwissDock found 160 conformations. By comparing the PPQ binding site with the sites of the 160 conformations, we consider that clusters 1, 6 and 33 could bind to the same site in a relatively similar way (Table 14). The inhibition constant (Ki) was calculated using the following formula: Ki = eˆ((∆G × 1000)/ (R × T)), where e = 2.7182, R = 1.98719 cal/(mol × K) (Regnault constant) and T = 298.15 K = 25 • C [159]. It can be seen that cluster 1 shows the lowest energy according to the ∆G value, but Ki and the maximum FullFitness value belong to complex 33. Visualisation and processing of the results obtained in the molecular docking study ( Figure 5) were performed using UCSF Chimera 1.15 [144,145]. The grid sizes used in SwissDock for cluster 1 are (x, y, z) = (15.5, 15.5, 20.5) with centre coordinates (x, y, z) = (−98, 13.711, −87.161). The inhibition constant (Ki) was calculated using the following formula: Ki = e^((ΔG × 1000)/(R × T)), where e = 2.7182, R = 1.98719 cal/(mol × K) (Regnault constant) and T = 298.15 K = 25 °C [159]. It can be seen that cluster 1 shows the lowest energy according to the ΔG value, but Ki and the maximum FullFitness value belong to complex 33. Visualisation and processing of the results obtained in the molecular docking study ( Figure 5) were performed using UCSF Chimera 1.15 [144,145]. The grid sizes used in SwissDock for cluster 1 are (x, y, z) = (15.5, 15.5, 20.5) with centre coordinates (x, y, z) = (−98, 13.711, −87.161).  To perform molecular docking using AutoDock Vina (a new version of the Webina online platform), the exhaustiveness of the search was set to 8 and the maximum energy difference to 3 kcal/mol. The space in which the test took place is represented by the volume of a cube (having the following dimensions: width = 20.4346, length = 27.864, height = 18.3759), and whose centre is defined by the coordinates x = −4.86256, y = −15.0503, z = −67.7222) [160]. Preparation for docking involves the insertion of hydrogen atoms on the chemical structure of both the ligand and the receptor molecule and the removal of the solvent. The protonation state corresponding to histidine was used, and Gasteiger partial charges were assigned ( Figure 6). To perform molecular docking using AutoDock Vina (a new version of the Webina online platform), the exhaustiveness of the search was set to 8 and the maximum energy difference to 3 kcal/mol. The space in which the test took place is represented by the volume of a cube (having the following dimensions: width = 20.4346, length = 27.864, height = 18.3759), and whose centre is defined by the coordinates x = −4.86256, y = −15.0503, z = −67.7222) [160]. Preparation for docking involves the insertion of hydrogen atoms on the chemical structure of both the ligand and the receptor molecule and the removal of the solvent. The protonation state corresponding to histidine was used, and Gasteiger partial charges were assigned ( Figure 6). Figure 6. Hydrogen bonds made between the ligand molecule (4Db6 compound) and the threonine residue of the receptor molecule. Visualisation of the ligand inserted into the "binding pocket" [144,145,[150][151][152]155,157].
The molecular docking results performed with AutoDock Vina are shown in Table  15, and the corresponding figures are presented in Figure S2 (Supplementary Materials). We chose to work further with model no.1 due to its low free energy (−6.3 kcal/mol) and root-mean-square deviation (RMSD) values that were below 2 Å. The 2 Å limit is often used as a criterion for predicting the correct binding site. The RMSD for two structures, a and b, of an identical molecule can be defined as follows: RMSDab = max(RMSD′ab, RMSD′ba), where rij represents the interatomic distance and the sum is over all N HA in structure a; the minimum is over all atoms in structure b with the same element type as the atom in structure a. RMSD is a measure of the distance between experimental and predicted structures that takes into account symmetry, partial symmetry (e.g., within a rotating Figure 6. Hydrogen bonds made between the ligand molecule (4Db6 compound) and the threonine residue of the receptor molecule. Visualisation of the ligand inserted into the "binding pocket" [144,145,[150][151][152]155,157].
The molecular docking results performed with AutoDock Vina are shown in Table 15, and the corresponding figures are presented in Figure S2 (Supplementary Materials). We chose to work further with model no.1 due to its low free energy (−6.3 kcal/mol) and root-mean-square deviation (RMSD) values that were below 2 Å. The 2 Å limit is often used as a criterion for predicting the correct binding site. The RMSD for two structures, a and b, of an identical molecule can be defined as follows:  (1) where r ij represents the interatomic distance and the sum is over all N HA in structure a; the minimum is over all atoms in structure b with the same element type as the atom in structure a. RMSD is a measure of the distance between experimental and predicted structures that takes into account symmetry, partial symmetry (e.g., within a rotating branch), and near-symmetry [160][161][162][163][164].
The main residues in the D subunit of the GS enzyme involved in interactions (within 1.49-2.81 Å) with the 4Db6 ligand are THR-223 (2 bonds) and GLU-129. Hydrogen bond connections play a key role in determining protein-ligand interactions [160,165]. In addition, the first conformation shows four active torsions: between C4 and P8, CA6 and C7, P8 and C9, and P8 and O11 [160].

Materials and Methods
Several series of analogous compounds (123 derivatives) have been theoretically designed based on the structure of glutamic acid to build a compound library of glutamic acid derivatives. From simple structure groups to more complex molecules, the chemical structures of the compounds were designed using BIOVIA Draw 21.1. [166]. The number of 123 compounds was reached after analysing the structure of glutamic acid to make as many specific structural modifications as possible. The classes of compounds and the structural changes made to the fundamental molecule were selected following the information found in the scientific literature. Our purpose was initially to design as many structural derivatives as possible because, after characterising and selecting these compounds based on wellestablished steps, we would be left with as many derivatives with optimal properties as possible to study further.
We also used the same software to generate the computational descriptors. To select suitable candidates for our purpose, we evaluated some properties of the molecules and their behaviour in the human organism. Physico-chemical characterisation of the desired compounds was carried out using SwissADME [92,103] and MarvinSketch [105]. Water solubility was tested using AquaSol [90], Chemicalize [91] and SwissADME [167,168]. Lipophilicity was analysed using SwissADME to determine the partition coefficients [169][170][171][172][173][174]. Toxicity was assessed using Toxtree [93] by applying the Cramer rules and the Kroess and Verhaar scheme, and GUSAR [175] was used to evaluate the acute toxicity in rodents.
Pharmacokinetic properties were analysed in terms of permeability and interactions with P-glycoprotein (P-gp) and some isoforms of cytochrome P450 using the SwissADME program. In addition, we evaluated the "drug-likeness" criteria according to Lipinski, Ghose, Veber, Egan, and Muegge rules using MarvinSketch, Chemicalize and DruLiTo [97].
The metabolism of the compounds was assessed using Toxtree, SmartCyp [131], and SOMP [132] and the bioactivity was evaluated using Molinspiration [134] and SWISSTarget prediction [133] (to predict the most probable molecular targets). The anticarcinogenic effect was assessed with the CLC-Pred software (Version 2.0) [135], which estimates in silico the cytotoxic effect based on the structural formula; the mechanism of action and adverse/toxic effects were tested using PASSonline [137].
Considering all the computed properties and their biological potential, "lead" compounds were selected.
We also attempted to validate our experimental procedures using positive and negative controls. Therefore, we chose methionine sulfoximine and phosphinothricin as positive controls for their proven activity of inhibiting glutamine synthetase [66,149]. As a negative control, we initially thought of glutamic acid, being the parent molecule for our derivatives [179]. However, it was interesting to observe that, according to the CLC-Pred software, it can show cytotoxic activity on four cell lines [135]. Therefore, in the end, we chose ampicillin as the negative control, which, according to the software, does not show cytotoxicity in any cancer cell line. All compounds were characterized using the previously described platforms and programs, passing through the same steps as the designed glutamic acid derivatives. Molecular docking was assessed using the ProteinsPlus online platform [180]. The results are presented in Tables S13-S18 and Figure S3 (Supplementary Materials).
To increase the accuracy of the study, molecular docking was carried out using several programs since they provided us with different information. PatchDock/ProteinPlus indicated the most probable binding sites in the protein's structure, calculated the surface area available for ligand binding, and generated the grid-box coordinates. Autodock Vina used these data and refined them, generating the values of ligand affinity for the target molecule and the distance from the RMSD lower bound and RMSD upper bound. It also showed the active torsions between atoms. Finally, SwissDock generated additional information, such as deltaG values and FullFitness, which were used to calculate the inhibition constant Ki.

Conclusions
Glutamic acid is an amino acid that can be manipulated very easily, and molecules derived from its structure have great potential to become innovative drugs. Of the 123 new GLA derivatives, 9 molecules proved to have biological potential, but more studies and optimisation are needed. The selected compounds show cytotoxicity against breast adenocarcinoma, lung cancer cell lines, colon carcinoma, and T cells from acute leukaemia. Compound 2Ba5 exhibited the lowest toxicity, while derivative 4Db6 exhibited the most intense bioactivity and could act like an ion channel modulator, protease inhibitor or enzyme inhibitor. A molecular docking study determined the binding site of the 4Db6 compound in the GS structure, D subunit, and found cluster 1 to be the most promising, having the lowest free energy value. Since compounds 5Aa1-5Ea3 were eliminated due to their increased toxicity, it is most probable that a single glutamic acid residue bound to the parent molecule cannot reduce the side effects or increase its biological activity. The toxicity of these compounds did not change significantly compared with the parent molecules, except for 7-hydroxynuciferine derivatives, which showed a higher risk of irritation, negative effects on the reproductive system, genotoxic carcinogenicity, tumorigenesis, and a higher risk of mutagenicity compared with 7-hydroxynuciferine. On the other hand, GLA-lycorine and GLA-dehydrolycorine complexes were less irritating to the skin than lycorine and dehydrolycorine, according to data provided by Toxtree and OSIRIS (Table S12; Supplementary Material). Further studies can be performed using these plant-derived molecules combined with more glutamic acid residues or poly-L glutamic acid to obtain more favourable results.
Based on the results provided by Molinspiration and CLC-Pred, further studies can be performed on other enzymes, ion channels, or proteases specific to the colon HCT-116 carcinoma cell line to simulate an interaction with the tumour itself. By marking isotopes at carbon 9 (bonded to the phosphorus atom) in the structure of 4Db6, the molecule can be analysed as a radiopharmaceutical compound (radioligand) as a potential candidate for anti-cancer therapy.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/molecules28104123/s1, Table S1: Programs and tested parameters; Table S2: Chemical structures, ID codes, and computational descriptors of glutamic acid derivatives obtained with Biovia Draw; Table S3: Structural and physicochemical properties: protonation and electric charge; Table S4: Water solubility; Table S5: Lipophilicity-partition coefficients; Table S6: Toxicity-I: Cramer rules, Kroess and Verhaar scheme; Table S7: Toxicity-II. Carcinogenic (genotoxic and non-genotoxic) and mutagenic effects evaluated using two different apps (Toxtree and OSIRIS); Table S8: Toxicity-III. Irritant/corrosive effect on the skin and eyes, effect on the reproductive system, biodegradability, and protein and DNA binding alerts, as assessed using Toxtree and OSIRIS; Table S9: Permeability and interactions with P-gp. Enzyme inhibitory effect on isoforms of cytochrome P450; Table S10: The number of broken rules, according to Lipinski, Ghose, Veber, Egan, and Muegge and the bioavailability score, the drug-likeness score, the lead-likeness score and the synthetic accessibility score; Figure S1: Homododecameric structure of the bacterial GS enzyme and D subunit; Figure S2: Molecular docking results visualised using UCSF Chimera and AutoDock Vina (Webina); Table S11: Chemical structures of colchicine, neferine, 7-hydroxynuciferine, lycorine, and 5,6-dehydrolycorine; Table S12: Toxicity comparison of vegetal compounds and their complexes with glutamic acid; Table S13: Characterization of phosphinothricin, methionine sulfoximine, glutamic acid, and ampicillin; Table S14: Molecular dynamics simulation results for phosphinothricin, methionine sulfoximine, glutamic acid, and ampicillin; Table S15: Molecular docking results for phosphinothricin, methionine sulfoximine, glutamic acid, and ampicillin; Figure S3: Molecular docking results. Interaction with glutamine synthetase of (a) hosphinothricin, (b) methionine sulfoximine, (c) glutamic acid, and (d) ampicillin; Table S16: Grid sizes used in Swissdock and energetic values of the most probable ligand-receptor complexes for phosphinothricin, methionine sulfoximine, glutamic acid, and ampicillin; Table S17: Molecular docking results obtained using AutoDock Vina for phosphinothricin, methionine sulfoximine, glutamic acid, and ampicillin; Table S18: Grid sizes used in AutoDock Vina for phosphinothricin, methionine sulfoximine, glutamic acid and ampicillin.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.
Sample Availability: Not applicable.