Focused Design of Novel Cyclic Peptides Endowed with GABARAP-Inhibiting Activity

(1) Background: Disfunctions in autophagy machinery have been identified in various conditions, including neurodegenerative diseases, cancer, and inflammation. Among mammalian autophagy proteins, the Atg8 family member GABARAP has been shown to be greatly involved in the autophagy process of prostate cancer cells, supporting the idea that GABARAP inhibitors could be valuable tools to fight the progression of tumors. (2) Methods: In this paper, starting from the X-ray crystal structure of GABARAP in a complex with an AnkirinB-LIR domain, we identify two new peptides by applying in silico drug design techniques. The two ligands are synthesized, biophysically assayed, and biologically evaluated to ascertain their potential anticancer profile. (3) Results: Two cyclic peptides (WC8 and WC10) displayed promising biological activity, high conformational stability (due to the presence of disulfide bridges), and Kd values in the low micromolar range. The anticancer assays, performed on PC-3 cells, proved that both peptides exhibit antiproliferative effects comparable to those of peptide K1, a known GABARAP inhibitor. (4) Conclusions: WC8 and WC10 can be considered new GABARAP inhibitors to be employed as pharmacological tools or even templates for the rational design of new small molecules.


Introduction
Autophagy plays a fundamental role in cellular, tissue, and organismal homeostasis. Metabolic stress (starvation and hypoxia) or the presence of dangerous cellular components, including dysfunctional organelles, intracellular microbes, and pathogenic proteins, can activate this pathway. Briefly, a multistep process, starting with the assembly of the phagophore, mediates the sequestration of organelles and proteins into the autophagosome. Its subsequent fusion with a lysosome leads to the formation of the autolysosome, in which the autophagosome content is degraded by lysosomal hydrolases [1].
More than 50 proteins (called Atgs) are involved in the autophagy machinery, but those responsible for the formation of the autophagosome and for cellular trafficking are members of the Atg8 family. In mammals, Atg8 proteins (mAtg8) are further divided into two subfamilies: GABARAP (GABA-A receptor-associated protein) and MAP1LC3 (microtubule-associated protein 1 light chain 3), or simply LC3. The former comprises GABARAP, GARAPL1, and GABARAPL2, while the latter includes LC3A (with the two splicing variants LC3Aα and LC3Aβ), LC3B, LC3B2, and LC3C [2]. Proteins belonging to we are confident that this study will open new avenues to identify the chemical entities endowed with significant therapeutic effects on this widespread disease.

Results and Discussion
Initially, the GABARAP/AnkB-LIR computational model was generated, starting from the three-dimensional data of the complex, available in the PDB. The simulation, refined by molecular modeling techniques (see Section 3 for details), was then used to design new peptides endowed with high affinity for GABARAP. The following procedure was adopted: 1.
Identification of the minimal AnkB-LIR sequence (core) responsible for the interaction with GABARAP; 2.
Mutation of the core sequence aimed at improving the theoretical affinity of the resulting peptides; 3.
Rigidification of the most promising peptides by adding disulfide bonds; 4.
Optimization of the peptide sequence by the addition of residues potentially occupying additional GABARAP binding pockets; 5.
Estimation of the binding affinity of the peptides by biophysical experiments; 6.
Evaluation of the killing effect on prostate cancer cells, exerted by the most promising candidates.

Setup of the GABARAP Computational Model and the Identification of the AnkB-LIR Core Sequence
The GABARAP/AnkB-LIR complex model ( Figure 1A) was retrieved from the Protein Data Bank (PDB) (accession code 5YIR) [11] and then refined by energy minimization and molecular dynamics (MDs) simulations, following the procedure reported in Section 3. The AnkB-LIR peptide rapidly reached the geometrical stability over the 500 ns long MD simulation, as demonstrated by the protein Cα RMSD plot ( Figure 1B). we are confident that this study will open new avenues to identify the chemical entities endowed with significant therapeutic effects on this widespread disease.

Results and Discussion
Initially, the GABARAP/AnkB-LIR computational model was generated, starting from the three-dimensional data of the complex, available in the PDB. The simulation, refined by molecular modeling techniques (see Section 3 for details), was then used to design new peptides endowed with high affinity for GABARAP. The following procedure was adopted: 1. Identification of the minimal AnkB-LIR sequence (core) responsible for the interaction with GABARAP; 2. Mutation of the core sequence aimed at improving the theoretical affinity of the resulting peptides; 3. Rigidification of the most promising peptides by adding disulfide bonds; 4. Optimization of the peptide sequence by the addition of residues potentially occupying additional GABARAP binding pockets; 5. Estimation of the binding affinity of the peptides by biophysical experiments; 6. Evaluation of the killing effect on prostate cancer cells, exerted by the most promising candidates.

Setup of the GABARAP Computational Model and the Identification of the AnkB-LIR Core Sequence
The GABARAP/AnkB-LIR complex model ( Figure 1A) was retrieved from the Protein Data Bank (PDB) (accession code 5YIR) [11] and then refined by energy minimization and molecular dynamics (MDs) simulations, following the procedure reported in Section 3. The AnkB-LIR peptide rapidly reached the geometrical stability over the 500 ns long MD simulation, as demonstrated by the protein Cα RMSD plot ( Figure 1B). representation of the GABARAP/AnkB-LIR complex, as derived from the X-ray structure (PDB accession code 5YIR). The protein surface is colored depending on the atomic partial charges of the protein residues: blue for positive and red for negative charges. The AnkB-LIR peptide is represented as cyan sticks. (B) Plot of the protein and ligand (AnkB-LIR) Cα atoms RMSD over the simulation time.
As expected, and verified by inspecting the MD trajectory frames, the LIR domain (sequence WVIV of AnkB-LIR) created numerous contacts with GABARAP (Supplementary Materials Figure S1). Interestingly, the SDEE residues were also involved in productive contacts, including the electrostatic interactions between the side chains of the peptide glutamates and the positively charged area of GABARAP close to K46 and R47. Conversely, the remaining C-terminal residues were mainly involved in internal contacts, stabilizing the α helix shaped by the DEEIEEARQKA sequence. representation of the GABARAP/AnkB-LIR complex, as derived from the X-ray structure (PDB accession code 5YIR). The protein surface is colored depending on the atomic par-tial charges of the protein residues: blue for positive and red for negative charges. The AnkB-LIR peptide is represented as cyan sticks. (B) Plot of the protein and ligand (AnkB-LIR) Cα atoms RMSD over the simulation time.
As expected, and verified by inspecting the MD trajectory frames, the LIR domain (sequence WVIV of AnkB-LIR) created numerous contacts with GABARAP (Supplementary Materials Figure S1). Interestingly, the SDEE residues were also involved in productive contacts, including the electrostatic interactions between the side chains of the peptide glutamates and the positively charged area of GABARAP close to K 46 and R 47 . Conversely, the remaining C-terminal residues were mainly involved in internal contacts, stabilizing the α helix shaped by the DEEIEEARQKA sequence.
Then, to exactly define the minimal portion of AnkB-LIR with the highest affinity for GABARAP, the peptides AnkB-LIR and AnkB-core (sequence WVIVSDEE) were subjected to MD simulations and MM-GBSA calculations to estimate their binding free energy.
Desmond and Prime tools of Maestro were employed for these computations, which predicted ∆G* values of −135.1 and −107.9 kcal/mol for AnkB-LIR and AnkB-core, respectively (Table 1). This result indicates that the 8 amino acids belonging to the AnkB-core contribute 75% of the overall interaction energy of the full AnkB-LIR peptide (composed of 20 residues). Then, to further define the contact area, MD simulations and MM-GBSA calculations were performed on the GABARAP/WVIV complex model, in which only the LIR motif (AnkB-wviv peptide) was simulated. As reported in Table 1, residues WVIV contribute 65% of the overall binding free energy. This outcome confirms that the LIR motif, shared by all proteins involved in the autophagy machinery, displays the highest complementarity with the GABARAP-binding site and is responsible for the most significant protein-protein contacts. Subsequently, the same computational protocol (MD simulations and MM-GBSA calculations) was applied to study the interaction of peptide K1. Considering that its experimentally determined K d lies in the nanomolar range (390 nM), this peptide could be considered as a reference inhibitor of GABARAP, together with AnkB-LIR. By our computations, the predicted ∆G* for peptide K1 on GABARAP was −118.9 kcal/mol (Table 1), a value slightly higher than that of AnkB-LIR (−135.1 kcal/mol). This result is in line with the K d values reported for the two peptides (0.27 and 390 nM, respectively).

Computational Design of AnkB-Core Analogs
Considering that the WVIVSDEE (AnkB-core) sequence accounts for 75% of the GABARAP/AnkB-LIR contacts, we proceeded to design small peptides endowed with high affinity for GABARAP using AnkB-core as a template. In this attempt, only the residues of the LIR domain (WVIV, positions 2-5 of AnkB-LIR) of AnkB-core (WVIVSDEE) were mutated, because of their direct involvement in the interaction with GABARAP. In this challenging effort, we tried to optimize the peptide sequence, also shared by other GABARAP binders, to improve the ligand/protein complementarity and selectivity. To this end, the "affinity maturation protocol", implemented into the Prime module of the Maestro software, was utilized to mutate the VIVS residues into all possible natural amino acid combinations. To avoid the combinatorial explosion, the Monte Carlo optimization option was selected. By this option, 2000 peptides were randomly generated by Monte Carlo algorithm and the peptides with a maximum of three simultaneous mutations were accepted to create the output file containing 100 solutions. Then, the Prime module was also employed to establish whether the mutations led to a more favorable interaction with the biological counterpart, by calculating the mutant peptides binding free-energy (∆Affinity) values [13].
At the end of these calculations, we visually inspected the results for the first 100 peptides with the highest predicted affinity. Among the predicted peptides, we noted that 7 of them displayed ∆Affinity values lower than 2 kcal/mol with respect to the initial template (AnkB-core). In these peptides, position 2 was substituted by Arg, Glu, or Ile; positions 3 and 4 contained only Ile; while position 5 included only alkaline residues, such as His and Arg. Among them, only one candidate, WEIIHDEE, named Pep-sol4, was further investigated by MD simulations and MM-GBSA calculations, because it presented an interesting Glu in position 2. Through this acidic amino acid, the peptide could potentially interact with the positive area shaped by GABARAP-K 46 and GABARAP-R 67 (two conserved residues among Atg8 proteins). Moreover, GABARAP-K 46 is considered to be a universal gate-keeper, regulating the entrance of ligands interacting through the LIR motif [14]. The structural alignment of the GABARAP/Pep-sol4 complex to the GABARAPL2/UBA5 NMR structure (PDB accession code 6H8C) [15] confirmed that the glutamate in position 2 of Pep-sol4 could reproduce the interaction displayed by E 15 (GAMEIIHEDNEWGIELVSE) of the "ubiquitin-like modifier activating enzyme 5" (UBA5) with GABARAP-K 46 .
By applying the computational protocol previously adopted for the reference inhibitors, the binding free-energy value of Pep-sol4 was calculated to be slightly lower than that of AnkB-core (−103.3 vs. −101.4 kcal/mol), suggesting that the new peptide could bind GABARAP with a similar binding affinity ( Table 2). investigated by MD simulations and MM-GBSA calculations, because it presented an interesting Glu in position 2. Through this acidic amino acid, the peptide could potentially interact with the positive area shaped by GABARAP-K46 and GABARAP-R67 (two conserved residues among Atg8 proteins). Moreover, GABARAP-K46 is considered to be a universal gate-keeper, regulating the entrance of ligands interacting through the LIR motif [14]. The structural alignment of the GABARAP/Pep-sol4 complex to the GABARAPL2/UBA5 NMR structure (PDB accession code 6H8C) [15] confirmed that the glutamate in position 2 of Pep-sol4 could reproduce the interaction displayed by E15 (GAMEIIHEDNEWGIELVSE) of the "ubiquitin-like modifier activating enzyme 5" (UBA5) with GABARAP-K46. By applying the computational protocol previously adopted for the reference inhibitors, the binding free-energy value of Pep-sol4 was calculated to be slightly lower than that of AnkB-core (−103.3 vs. −101.4 kcal/mol), suggesting that the new peptide could bind GABARAP with a similar binding affinity ( Table 2). Moreover, the visual inspection of the GABARAP/Pep-sol4 MD trajectory and the root mean square fluctuation (RMSF) plot of the ligand atoms over the simulation time suggested that the C-terminal portion of the peptide was not firmly bound to the GABA-RAP surface, thus preventing a stable and productive interaction with the protein ( Figure  2). Consequently, considering that the side chains of I3 and D8 were spatially close in the binding mode adopted by Pep-sol4, we designed a cyclic peptide in which I3 and D8 were mutated into two Cys residues bound by a disulfide bond. This modification aimed at reducing the conformational flexibility of the ligand, generating a more stable binding mode on the GABARAP surface. The resulting peptide (named Pep-sol4cc, WECIHDEC) was again analyzed in the complex with GABARAP by MD simulations and MM-GBSA calculations. At the end of these computational procedures, the estimated ΔG* of Pep-sol4cc was −103.7 kcal/mol, a value comparable to that of Pep-sol4 (−103.3 kcal/mol). This information led us to conclude that the structural rigidification did not affect the affinity of the peptide; however, as demonstrated by the ligand RMSF plot (Figure 2), an improvement of the conformational stability was successfully achieved. investigated by MD simulations and MM-GBSA calculations, because it presented an interesting Glu in position 2. Through this acidic amino acid, the peptide could potentially interact with the positive area shaped by GABARAP-K46 and GABARAP-R67 (two conserved residues among Atg8 proteins). Moreover, GABARAP-K46 is considered to be a universal gate-keeper, regulating the entrance of ligands interacting through the LIR motif [14]. The structural alignment of the GABARAP/Pep-sol4 complex to the GABARAPL2/UBA5 NMR structure (PDB accession code 6H8C) [15] confirmed that the glutamate in position 2 of Pep-sol4 could reproduce the interaction displayed by E15 (GAMEIIHEDNEWGIELVSE) of the "ubiquitin-like modifier activating enzyme 5" (UBA5) with GABARAP-K46. By applying the computational protocol previously adopted for the reference inhibitors, the binding free-energy value of Pep-sol4 was calculated to be slightly lower than that of AnkB-core (−103.3 vs. −101.4 kcal/mol), suggesting that the new peptide could bind GABARAP with a similar binding affinity (Table 2). Moreover, the visual inspection of the GABARAP/Pep-sol4 MD trajectory and the root mean square fluctuation (RMSF) plot of the ligand atoms over the simulation time suggested that the C-terminal portion of the peptide was not firmly bound to the GABA-RAP surface, thus preventing a stable and productive interaction with the protein ( Figure  2). Consequently, considering that the side chains of I3 and D8 were spatially close in the binding mode adopted by Pep-sol4, we designed a cyclic peptide in which I3 and D8 were mutated into two Cys residues bound by a disulfide bond. This modification aimed at reducing the conformational flexibility of the ligand, generating a more stable binding mode on the GABARAP surface. The resulting peptide (named Pep-sol4cc, WECIHDEC) was again analyzed in the complex with GABARAP by MD simulations and MM-GBSA calculations. At the end of these computational procedures, the estimated ΔG* of Pep-sol4cc was −103.7 kcal/mol, a value comparable to that of Pep-sol4 (−103.3 kcal/mol). This information led us to conclude that the structural rigidification did not affect the affinity of the peptide; however, as demonstrated by the ligand RMSF plot (Figure 2), an improvement of the conformational stability was successfully achieved. investigated by MD simulations and MM-GBSA calculations, because it presented an interesting Glu in position 2. Through this acidic amino acid, the peptide could potentially interact with the positive area shaped by GABARAP-K46 and GABARAP-R67 (two conserved residues among Atg8 proteins). Moreover, GABARAP-K46 is considered to be a universal gate-keeper, regulating the entrance of ligands interacting through the LIR motif [14]. The structural alignment of the GABARAP/Pep-sol4 complex to the GABARAPL2/UBA5 NMR structure (PDB accession code 6H8C) [15] confirmed that the glutamate in position 2 of Pep-sol4 could reproduce the interaction displayed by E15 (GAMEIIHEDNEWGIELVSE) of the "ubiquitin-like modifier activating enzyme 5" (UBA5) with GABARAP-K46. By applying the computational protocol previously adopted for the reference inhibitors, the binding free-energy value of Pep-sol4 was calculated to be slightly lower than that of AnkB-core (−103.3 vs. −101.4 kcal/mol), suggesting that the new peptide could bind GABARAP with a similar binding affinity (Table 2). Moreover, the visual inspection of the GABARAP/Pep-sol4 MD trajectory and the root mean square fluctuation (RMSF) plot of the ligand atoms over the simulation time suggested that the C-terminal portion of the peptide was not firmly bound to the GABA-RAP surface, thus preventing a stable and productive interaction with the protein ( Figure  2). Consequently, considering that the side chains of I3 and D8 were spatially close in the binding mode adopted by Pep-sol4, we designed a cyclic peptide in which I3 and D8 were mutated into two Cys residues bound by a disulfide bond. This modification aimed at reducing the conformational flexibility of the ligand, generating a more stable binding mode on the GABARAP surface. The resulting peptide (named Pep-sol4cc, WECIHDEC) was again analyzed in the complex with GABARAP by MD simulations and MM-GBSA calculations. At the end of these computational procedures, the estimated ΔG* of Pep-sol4cc was −103.7 kcal/mol, a value comparable to that of Pep-sol4 (−103.3 kcal/mol). This information led us to conclude that the structural rigidification did not affect the affinity of the peptide; however, as demonstrated by the ligand RMSF plot (Figure 2), an improvement of the conformational stability was successfully achieved. investigated by MD simulations and MM-GBSA calculations, because it presented an interesting Glu in position 2. Through this acidic amino acid, the peptide could potentially interact with the positive area shaped by GABARAP-K46 and GABARAP-R67 (two conserved residues among Atg8 proteins). Moreover, GABARAP-K46 is considered to be a universal gate-keeper, regulating the entrance of ligands interacting through the LIR motif [14]. The structural alignment of the GABARAP/Pep-sol4 complex to the GABARAPL2/UBA5 NMR structure (PDB accession code 6H8C) [15] confirmed that the glutamate in position 2 of Pep-sol4 could reproduce the interaction displayed by E15 (GAMEIIHEDNEWGIELVSE) of the "ubiquitin-like modifier activating enzyme 5" (UBA5) with GABARAP-K46. By applying the computational protocol previously adopted for the reference inhibitors, the binding free-energy value of Pep-sol4 was calculated to be slightly lower than that of AnkB-core (−103.3 vs. −101.4 kcal/mol), suggesting that the new peptide could bind GABARAP with a similar binding affinity (Table 2). Moreover, the visual inspection of the GABARAP/Pep-sol4 MD trajectory and the root mean square fluctuation (RMSF) plot of the ligand atoms over the simulation time suggested that the C-terminal portion of the peptide was not firmly bound to the GABA-RAP surface, thus preventing a stable and productive interaction with the protein ( Figure  2). Consequently, considering that the side chains of I3 and D8 were spatially close in the binding mode adopted by Pep-sol4, we designed a cyclic peptide in which I3 and D8 were mutated into two Cys residues bound by a disulfide bond. This modification aimed at reducing the conformational flexibility of the ligand, generating a more stable binding mode on the GABARAP surface. The resulting peptide (named Pep-sol4cc, WECIHDEC) was again analyzed in the complex with GABARAP by MD simulations and MM-GBSA calculations. At the end of these computational procedures, the estimated ΔG* of Pep-sol4cc was −103.7 kcal/mol, a value comparable to that of Pep-sol4 (−103.3 kcal/mol). This information led us to conclude that the structural rigidification did not affect the affinity of the peptide; however, as demonstrated by the ligand RMSF plot (Figure 2), an improvement of the conformational stability was successfully achieved. Moreover, the visual inspection of the GABARAP/Pep-sol4 MD trajectory and the root mean square fluctuation (RMSF) plot of the ligand atoms over the simulation time suggested that the C-terminal portion of the peptide was not firmly bound to the GABARAP surface, thus preventing a stable and productive interaction with the protein (Figure 2). Consequently, considering that the side chains of I 3 and D 8 were spatially close in the binding mode adopted by Pep-sol4, we designed a cyclic peptide in which I 3 and D 8 were mutated into two Cys residues bound by a disulfide bond. This modification aimed at reducing the conformational flexibility of the ligand, generating a more stable binding mode on the GABARAP surface. The resulting peptide (named Pep-sol4cc, WECIHDEC) was again analyzed in the complex with GABARAP by MD simulations and MM-GBSA calculations. At the end of these computational procedures, the estimated ∆G* of Pep-sol4cc was −103.7 kcal/mol, a value comparable to that of Pep-sol4 (−103.3 kcal/mol). This information led us to conclude that the structural rigidification did not affect the affinity of the peptide; however, as demonstrated by the ligand RMSF plot (Figure 2), an improvement of the conformational stability was successfully achieved.

Computational Design of the WC8 and WC10 Peptides
Then, with the aim of improving the theoretical binding affinity of Pep-sol4cc, H5 was mutated into a Phe, supposing that it could better fill the hydrophobic area sized by P52, L55, and Q59. The resulting peptide (WC8, sequence WECIFDEC) was analyzed by MD simulations and MM-GBSA calculations, which led to a ΔG* value of 12 kcal/mol, lower

Computational Design of the WC8 and WC10 Peptides
Then, with the aim of improving the theoretical binding affinity of Pep-sol4cc, H 5 was mutated into a Phe, supposing that it could better fill the hydrophobic area sized by P 52 , L 55 , and Q 59 . The resulting peptide (WC8, sequence WECIFDEC) was analyzed by MD simulations and MM-GBSA calculations, which led to a ∆G* value of 12 kcal/mol, lower than that of the originator ( Table 2). The cluster analysis performed on the MD trajectory frames displayed that, in the structure representative of the most populated cluster of GABARAP/WC8 conformations (accounting for 79% of conformational ensembles explored), the ligand was stably bound on the GABARAP surface (see Figure 2 for the RMSF plot), forming numerous interactions (Table 3 and Figure 3A).  Figure 2. RMSF plots of AnkB-core analogs. Backbone atoms were considered in these calculations. The residues shared by all peptides are highlighted by capital letters.

Computational Design of the WC8 and WC10 Peptides
Then, with the aim of improving the theoretical binding affinity of Pep-sol4cc, H5 was mutated into a Phe, supposing that it could better fill the hydrophobic area sized by P52, L55, and Q59. The resulting peptide (WC8, sequence WECIFDEC) was analyzed by MD simulations and MM-GBSA calculations, which led to a ΔG* value of 12 kcal/mol, lower than that of the originator ( Table 2). The cluster analysis performed on the MD trajectory frames displayed that, in the structure representative of the most populated cluster of GABARAP/WC8 conformations (accounting for 79% of conformational ensembles explored), the ligand was stably bound on the GABARAP surface (see Figure 2 for the RMSF plot), forming numerous interactions (Table 3 and Figure 3A).   In detail, several H-bonds were established, and salt bridges formed between WC8-E 2 and the side chains of GABARAP-K 46 and -R 67 , and between the C ter of WC8-C 8 and the side chain of GABARAP-R 28 .
Regarding the hydrophobic contacts, the indole ring of WC8-W 1 was positioned in a pocket formed by residues I 23 , I 32 , P 30 , K 48 , and F 104 , while the side chain of WC8-I 4 pointed toward another pocket delimited by Y 49 , V 51 , F 60 , L 63 , and I 64 , establishing van der Waals (vdW) interactions. Finally, similar hydrophobic contacts were also observed between WC8-F 5 and the GABARAP area shaped by P 52 , L 55 , and L 63 . The complete 2D representation of the interaction network between WC8 and GABARAP is shown in the Supplementary Materials Figure S2A.
WC8 exhibited an estimated ∆G* value close to that of K1; hence, with the aim of designing a more potent peptide, we included two additional N-terminal residues. This hypothesis was supported by the fact that the AnkB-LIR peptide (EEWVIVSDEEIEEARQKA), used as a template, contains two glutamate residues before the AnkB-core (WVIVSDEE). For this reason, we speculated that the homologation of WC8 on the N ter could lead to a more potent peptide, considering that the new atoms could create an additional bond network. Our objective was to reach the region sized by I 32 , Y 5 , and K 47 , close to the W site, on the GABARAP surface (the yellow area in Supplementary Materials Figure S3A). Therefore, to find the optimal N-terminal sequence, two glycines were initially added to WC8 (GGWECIFDEC); then, the application of the "affinity maturation protocol" on the first Gly residue led to the identification of Tyr (YC10, Table 2) and Trp (WC10, Table 2) as the most suitable substitutions. In this attempt, the glycine in position 2 was not mutated to allow a certain conformational mobility on the N-terminal tail of the new peptide. Interestingly, the N-terminal residues (WG) and the Glu in position 4 (E 4 ) of WC10 (WGWECIFDEC) reproduced the interactions displayed by UBA5 (GAMEIIHEDNEWGIELVSE) in the complex with GABARAP and GABARAPL2 [15] (Supplementary Materials Figure S3B).
MD simulations and MM-GBSA calculations on the GABARAP/YC10 and GABARAP/ WC10 complexes revealed that the latter possessed the highest affinity, with a predicted ∆G* value almost 7 kcal/mol lower than that of WC8 (Table 2). Cluster analysis was then performed on the GABARAP/WC10 MD trajectory frames; the representative structure of the most populated cluster of conformations (which accounts for the 88% of total conformational ensemble explored) is represented in Figure 3B. Notably, the visual inspection of the GABARAP/WC10 most representative structure highlighted that the side chain of the newly added residue (W 1 ) formed a cation-π interaction with GABARAP-K 46 , while the carbonyl group of the same residue established an H-bond with the side chain of GABARAP-K 48 ( Figure 3B). Surprisingly, the side chain of W 1 did not occupy the expected region of GABARAP, but the new additional cation-π interaction greatly contributed to the calculated binding affinity of the peptide. In addition, WC10 (1) shares all the interaction networks established by the GABARAP/WC8 complex, (2) is able to orientate GABARAP-K 48 in order to establish additional cation-π interactions with WC10-W 3 , and (3) is able to form an additional H-bond interaction between the I 6 (C=O) and GABARAP-R 67 (=NH 2 + ) ( Figure 3). The 2D representation of the interaction network between WC10 and GABARAP is showed in the Supplementary Materials Figure S2B.
To conclude, we designed two new cyclic peptides (WC8 and WC10) endowed with a reduced conformational mobility and calculated binding free-energy values in a lower range than those estimated for the reference peptides, AnkB-core and K1. In light of these data, WC8 and WC10 could exhibit higher experimental affinities compared to the reference peptides. Nevertheless, it must be considered that our computations did not account for the entropic contributions to the binding free energy; hence, they should be regarded as a starting point for further experimental studies.

Experimental Validation and Biophysical Experiments
Based on the results of the computational study, the K1, AnkB-core, WC8, and WC10 peptides were purchased by Proteogenix for the experimental investigations. In detail, MST and SPR assays were conducted on the peptides displaying a sufficient stability in water and PBS buffer. Then, the anticancer potential of the most promising candidates was investigated. Initially, we verified that the peptides were water soluble and stable in the buffer in which the recombinant GABARAP protein was solved. Unfortunately, AnkB-core was not soluble in water; thus, it was impossible to use it as a reference peptide. Conversely, K1, WC8, and WC10 displayed an excellent stability in water and PBS, so they were tested by biophysical methods.
In detail, MST and SPR experiments were carried out with the aim of measuring their K d values on GABARAP. As a preliminary step, the K d of the reference peptide K1 was determined in order to (1)    Surprisingly, the Kd measured for K1 was close to 3 µM, a value 7 times higher than the one reported in the literature. However, all the techniques employed in this study agreed on this value. The data obtained for WC8 revealed a Kd of 22 µM ( Figure 5A,B), consistent among the different biophysical approaches. Remarkably, WC10 displayed a Kd in the same range of the reference peptide K1, with a value close to 4 µM, obtained by both MST and SPR ( Figure 5C,D). Surprisingly, the K d measured for K1 was close to 3 µM, a value 7 times higher than the one reported in the literature. However, all the techniques employed in this study agreed on this value. The data obtained for WC8 revealed a K d of 22 µM ( Figure 5A,B), consistent among the different biophysical approaches. Remarkably, WC10 displayed a K d in the same range of the reference peptide K1, with a value close to 4 µM, obtained by both MST and SPR ( Figure 5C,D). Because the Kd value of peptide K1 proved to be higher than the one reported in the literature, we decided to validate our data by repeating the MST-binding affinity experiments using another Monolith instrument (Monolith NT.115Pico), located in a different laboratory. The new results confirmed our previous findings, with all the peptides displaying Kd values consistent with those determined earlier (Supplementary Materials Figure S4).
According to the theoretical predictions, WC10 should be more active than K1 (ΔG* = −122.0 vs. −118.9 kcal/mol, respectively), and WC8 less active than the others (ΔG* = −115.7 kcal/mol), as shown in Table 2. Considering the confidence range of the experimental Kd and the omission of the entropic contribution to the estimated binding freeenergy values, the computations predicted the affinity trend of the selected peptides well. Because the K d value of peptide K1 proved to be higher than the one reported in the literature, we decided to validate our data by repeating the MST-binding affinity experiments using another Monolith instrument (Monolith NT.115Pico), located in a different laboratory. The new results confirmed our previous findings, with all the peptides displaying K d values consistent with those determined earlier (Supplementary Materials Figure S4).
According to the theoretical predictions, WC10 should be more active than K1 (∆G* = −122.0 vs. −118.9 kcal/mol, respectively), and WC8 less active than the others (∆G* = −115.7 kcal/mol), as shown in Table 2. Considering the confidence range of the experimental K d and the omission of the entropic contribution to the estimated binding freeenergy values, the computations predicted the affinity trend of the selected peptides well.

Biological Experiments
Finally, K1, WC8, and WC10 were tested in vitro on PC-3 prostate cancer cells, to evaluate their potential antitumor effects ( Figure 6). Prostate cancer is the second most commonly diagnosed malignancy in men worldwide. Considering that the probability of developing the disease during a man's lifetime is 15% and that prostate tumor cells can also spread to the lungs and bones via angiogenesis [16], we decided to evaluate the biological activity of the peptides in vitro on a prostate cancer model. PC-3 cells were chosen for the screening due to their highly metastatic nature, effectively mimicking an aggressive form of the disease. Notably, it has recently been demonstrated that prostate cancer models show a significant upregulation of autophagy [6,8,17]. Therefore, the biological activity of different concentrations of K1, WC8, and WC10 (from 0.5 to 10 µM) was evaluated with an MTS cell viability assay on PC-3 cells and noncancerous PNT2 prostate cells ( Figure 6). The results reported in Figure 6A show that, 96 h post-treatment, none of the tested samples displayed a significant cytotoxicity (cell availability > 90%), confirming the excellent biocompatibility and potential pharmacological selectivity for tumor cells. Indeed, as shown in Figure 6B, a reduction in cell viability (expressed as percentage % of viable cells) was observed in PC-3 cells treated with K1, WC8, and WC10 compared to the untreated control. Interestingly, the treatments of PC-3 cells with WC8 and WC10 (from 1 to 10 µM) display high efficacy, when compared to Paclitaxel ( Figure 6B). The in vitro data demonstrate that the compounds exhibited a considerable anticancer activity, especially at the highest tested concentration (cell viability 27.16% for K1, 24.06% for WC8, and 22.5% for WC10).
The biological data on PC-3 cells indicate that all peptides possess IC50 values close to 5 µM, consistent with the Kd estimated by the biophysical experiments. Surprisingly, WC8 displayed a better activity profile than the reference peptide K1. Based on this finding, we may speculate that some other biochemical mechanism or additional activity on different mAtg8 subfamilies could improve the activity of the new peptides [8]. Nevertheless, since the work presented here is a proof-of-concept study, the peptides have been preliminary tested in a non-cancerous and subsequently in a cancer cell line, to exclude possible off-target cytotoxicity and perform an initial pilot study to evaluate the in vitro anti-cancer efficacy. However, we are planning to extend the screening to other cancer cell lines in the upcoming further evaluation of the peptides and their antineoplastic mechanism of action. Furthermore, to shed light on the possible secondary targets, we have planned biological and biophysical experiments on LC3B to evaluate if our peptides could show any affinity to it. Moreover, we should remember that the data on PC-3 cells represent a preliminary assessment that merely suggests the potential application of GABA-RAP inhibitors as anticancer agents. Further biological assays are needed to unveil the mechanism by which these peptides trigger cell death. Therefore, the biological activity of different concentrations of K1, WC8, and WC10 (from 0.5 to 10 µM) was evaluated with an MTS cell viability assay on PC-3 cells and non-cancerous PNT2 prostate cells ( Figure 6). The results reported in Figure 6A show that, 96 h post-treatment, none of the tested samples displayed a significant cytotoxicity (cell availability > 90%), confirming the excellent biocompatibility and potential pharmacological selectivity for tumor cells. Indeed, as shown in Figure 6B, a reduction in cell viability (expressed as percentage % of viable cells) was observed in PC-3 cells treated with K1, WC8, and WC10 compared to the untreated control. Interestingly, the treatments of PC-3 cells with WC8 and WC10 (from 1 to 10 µM) display high efficacy, when compared to Paclitaxel ( Figure 6B). The in vitro data demonstrate that the compounds exhibited a considerable anticancer activity, especially at the highest tested concentration (cell viability 27.16% for K1, 24.06% for WC8, and 22.5% for WC10).
The biological data on PC-3 cells indicate that all peptides possess IC 50 values close to 5 µM, consistent with the K d estimated by the biophysical experiments. Surprisingly, WC8 displayed a better activity profile than the reference peptide K1. Based on this finding, we may speculate that some other biochemical mechanism or additional activity on different mAtg8 subfamilies could improve the activity of the new peptides [8]. Nevertheless, since the work presented here is a proof-of-concept study, the peptides have been preliminary tested in a non-cancerous and subsequently in a cancer cell line, to exclude possible offtarget cytotoxicity and perform an initial pilot study to evaluate the in vitro anti-cancer efficacy. However, we are planning to extend the screening to other cancer cell lines in the upcoming further evaluation of the peptides and their antineoplastic mechanism of action. Furthermore, to shed light on the possible secondary targets, we have planned biological and biophysical experiments on LC3B to evaluate if our peptides could show any affinity to it. Moreover, we should remember that the data on PC-3 cells represent a preliminary assessment that merely suggests the potential application of GABARAP inhibitors as anticancer agents. Further biological assays are needed to unveil the mechanism by which these peptides trigger cell death.

Materials and Methods
3.1. The GABARAP Computational Model, General Protocol for MD Simulations, and ∆G* Estimation The computational model utilized in this study was built using the 3D coordinates of the GABARAP/AnkB-LIR complex, available in the PDB, accession code 5YIR [11]. The chains A (GABARAP) and G (AnkB-LIR) were chosen for the creation of the starting complex model. The resulting model was then optimized by the "protein preparation wizard" module implemented in the Maestro software (release 2019-4, Schrödinger, LLC, New York, NY, USA, 2017). This tool permitted us to (1) check the protonation state of the residues at pH 7.4, (2) verify the residue completeness, (3) eliminate atomic clashes, and (4) assign the OPLS3e force field [18]. Then, a cubic box of water molecules (almost 7000, represented by the TIP3P model [19]) was built around the protein-ligand complex and energy minimized by the Desmond algorithm implemented in Maestro. A single run of 500 ns MD simulations was performed, again utilizing the Desmond algorithm, and the "simulation interactions diagram" tool was employed to evaluate the stability of the peptide interacting with GABARAP (see Figure 1B for the GABARAP/AnkB-LIR Cα atoms RMSD variation over the simulation time).
This protocol was also applied for the model setup and MD simulations on the GABARAP complexes resulting from the variation of the length or the mutation of the AnkB-LIR peptide (see Supplementary Materials Figure S5, for the Cα atom RMSD plots of all the other simulated systems). The calculations of the binding free energy values were accomplished by the Prime algorithm [20] available on the Maestro platform, by applying the MM-GBSA algorithm. In these calculations, the single-trajectory approach was applied and the entropy contributions to the binding free energy was neglected. For this reason, the estimated binding free-energy values are termed ∆G* by us. The affinity maturation functionality, implemented in the Bioluminate/Prime module of Maestro, estimated the change in the affinity (∆Affinity) between the mutant peptides and GABARAP. Finally, the peptides achieving the highest gain in ∆Affinity were additionally refined by MD simulations and MM-GBSA calculations, by applying the previously described protocol. The cluster analysis of the MD trajectories was accomplished by the "Desmond trajectory clustering" tool of Maestro, setting the creation of at least 5 clusters. In these calculations, the RMSD of the backbone atoms was used to create the matrix, which was then analyzed through the affinity propagation clustering method [21].

Microscale Thermophoresis (MST)
The binding affinity (K d ) between the target protein (6-His-tagged recombinant GABARAP, produced by Abcam, Cambridge, UK) and peptide ligands (K1, WC8 and WC10, acquired by ProteoGenix, Schiltigheim, France) was measured by MST [22]. Briefly, histidine-tagged GABARAP was labeled by using a His-tag-specific dye for 30 min at room temperature, using either the conventional Monolith His-Tag Labeling Kit RED-tris-NTA (MO-L008) or the newer Monolith His-Tag Labeling Kit RED-tris-NTA 2nd Generation (MO-L018), both purchased from NanoTemper Technologies (GmbH, München, Germany). MST experiments were performed on both a Monolith NT.115 and a Monolith NT.115Pico instrument (NanoTemper Technologies, München, GmbH), in order to further increase the statistical significance of the attained results.
A fixed concentration of the labeled GABARAP protein was mixed with 16 1:1 serial dilutions of K1, WC8, and WC10. The protein and the peptide were incubated for 60-90 min at room temperature. The MST analysis was performed using standard capillaries, with the following experimental settings: an MST power of 40% (to create a temperature gradient) and a different excitation power in relation to the considered system and instrument (see Supplementary Materials for details, Table S1). K d values were calculated from com-pound concentration-dependent changes of normalized fluorescence (Fnorm). In all the experiments, both the protein and the peptides were dissolved in PBS-T buffer (phosphatebuffered saline + 0.05% Tween™ 20; NanoTemper Technologies, GmbH) and 2.5% DMSO. The auto-fluorescence of each peptide was evaluated before proceeding to the determination of the K d . At least two independent experiments were performed to calculate the K d values. Data were processed using the NanoTemper MO.Affinity Analysis software (version 2.3), and the fitting was performed by using the K d model.

Conclusions
In this study, starting from the GABARAP/AnkB-LIR X-ray crystal structure, we created an affordable GABARAP/AnkB-LIR complex computational model. This was utilized to investigate the role played by different regions of the AnkB-LIR sequence. Then, by applying integrated computational techniques (MD simulations, MM-GBSA calculations, and affinity maturation), we designed two peptides (WC8 and WC10) endowed with theoretical affinities in line with the ones predicted for the reference peptides AnkB-core and K1. The experimental measurement of the affinity values led us to prove that WC10 (2 residues shorter and more rigid than K1) displays a biological activity similar to that of K1. MST, SPR, and in vitro assays on PC-3 cells confirmed this observation. In our opinion, this study has the potential to open new avenues of research towards the design of novel anticancer compounds, employing WC10 as a template. Our results confirm again that a suitable interference with the autophagy process of cancer cells can represent an innovative and viable therapeutic strategy. Consequently, we are confident that the discovery of new potent and specific autophagy modulators will become increasingly important in the near future [23].