Computational Insights into Novel Inhibitor N-(3-(tert-Butylcarbamoyl)-4-methoxyphenyl)-indole and Ingliforib Specific against GP Isoenzyme Dimers Interaction Mechanism

The high conservation of the three subtypes of glycogen phosphorylase (GP) presents significant challenges for specific inhibitor studies targeting GP. Our prior screening revealed that compound 1 exhibited unequal inhibitory activity against the three GP subtypes, with a noticeable effect against brain GP (PYGB). The commercially available ingliforib demonstrated potent inhibitory activity specifically against liver GP (PYGL). To guide the further design and screening of high-specificity inhibitors, the possible reasons for the differential inhibitory activity of two compounds against different GP subtypes were analyzed, with ingliforib as a reference, through molecular docking and molecular dynamics simulations. Initially, the study predicted the binding modes of ligands with the three GP receptor subtypes using molecular docking. Subsequently, this was validated by molecular dynamics experiments, and possible amino acid residues that had important interactions were explored. The strong correlation between the calculated interaction free energies and experimental inhibitory activity implied the reasonable binding conformations of the compounds. These findings offer insight into the different inhibitory activity of compound 1 and ingliforib against all three GP subtypes and provide guidance for the design of specific target molecules that regulate subtype selectivity.


Introduction
Glycogen phosphorylase (GP) is a critical rate-limiting enzyme that catalyzes the breakdown of glycogen into glucose-1-phosphate during glycogenolysis. This reaction is vital in maintaining glucose homeostasis, which is necessary to sustain normal cellular function. Therefore, GP was initially considered a potential target for the treatment of type 2 diabetes (T2D) [1]. Recently, GP inhibitors have also been found to have a positive effect on myocardial ischemia, cerebral ischemia and cancer [2][3][4]. In mammals, GP is a family of three isozymes, named brain GP (PYGB), liver GP (PYGL) and muscle GP (PYGM), according to the tissue in which they are expressed, and they exist as active or inactive homodimers ( Figure 1) [5]. Studies have shown that the different GP isoforms usually prioritize different types of physiological activity; for example, the inhibition of PYGL lowers blood glucose [6] and the inhibition of PYGB facilitates myocardial protection [7], while PYGM is associated with muscle GP deficiency (McArdle disease), schizophrenia and cancer [8]. However, these isozymes were found to be 80% homologous [9], and cross-reactivity may occur when non-selective GP inhibitors are used [10]. Thus, one of the challenges in developing inhibitors of GP is to achieve sufficient selectivity for enzyme isoforms.

Binding Model Studies by Molecular Docking
To investigate the associations between compound 1 and various GP subtypes, molecular docking was carried out. Initially, the GP protein structure was acquired from the protein database (refer to Materials and Methods) considering the following criteria: ① human origin; ② conformational resolution ≤ 2.5 Å, minimized to the greatest extent; ③ maximum possible conformational sequence coverage; ④ crystalline pH value as close to the normal physiological range of the human body as feasible. Ensuring that the above four conditions were satisfied to the greatest extent possible, the GP protein structures with PDB IDs of 5IKO (PYGB), 2ZB2 (PYGL) and 1Z8D (PYGM) were chosen for the molecular docking experiments in this study, following appropriate adjustments. Upon completion of molecular docking, the outcomes were examined based on the subsequent aspects.

Spatial Position and Orientation
After docking, the conformation with the lowest binding energy in GP was selected as the probable binding conformation. The spatial position and orientation of the ligand ingliforib and compound 1 in the three GP subtype receptors were observed. As shown in Figure 2A, ingliforib binds within the central cavity consisting of two identical subunits. In the complexes of the three subtypes, the indole, benzyl and dihydroxypyrrolidine moieties of the ligand molecule point towards three different pockets ( Figure 2C-E). Specifically, the indole moiety of ingliforib is deeply buried within pocket 1, a closed hydrophobic cavity [19]. The dihydroxypyrrolidine moiety extends towards a wide pocket 2, while the benzyl moiety points towards pocket 3. The ligand molecule adopts a "Y"-shaped conformation at the dimer interface. Comparison of the three subtype complexes using Pymol indicates that their spatial conformations are similar. The differences are mainly observed in dihydroxypyrrolidine, located in pocket 2, which has a different spatial orientation.
For compound 1, it also binds in a central cavity ( Figure 2B) [12,20]. Alignment of the three subtype complexes by Pymol reveals that the spatial location and orientation of the With advances in structural biology and computational modeling providing new insights into the mechanisms of GP inhibition, this could inform the design of new inhibitors with better potency and selectivity. The crystal structure of PYGM was first resolved in the early 1970s [11]. Many structures of this isoenzyme were then obtained in the presence of different metastable effectors and/or drugs to better understand the metastable regulation of PYGM. Later, in 2000, the crystallographic structure determination of PYGL revealed the structural basis for the differences in the regulation of these isozymes [12]. In 2016, the crystal structure of human PYGB was obtained for the first time by Mathieu et al. (PDB ID: 5IKO) [13]. Comparison of the structures of the three isozymes showed that the structure of PYGB is highly similar to the active state of PYGM [14]. This may be the reason that our previous compound 2, with its IC 50 value for PYGB, was closer to that of PYGM compared to the IC 50 value for PYGL [3]. Meanwhile, they also show some differences, such as the tower helix (helix 7 of each monomer), which is a major component of the dimerization interface in GP isozymes, and these two helices show antiparallel binding and control the dimerization and activation of the enzyme. They show typical crossover angles of 75 • , 84 • and 45 • in active PYGB, PYGM and PYGL, respectively [14]. Freeman et al. explored the potential selectivity of GPi688 and GPi819 against the liver and muscle GP isoforms and showed that the sensitivity of the inhibitors was dependent on the activation state of the enzyme [15]. To avoid cross-reactivity between GP isoforms, Konkimalla et al. used pentacyclic triterpenes as a ligand and explored the design of selective inhibitors for the liver and muscle isoforms via molecular docking techniques [16]. Pfizer reported a series of GP inhibitors containing indole structures acting on the new allosteric site. Ingliforib (CP 368296; GPi 296) is one of them and has entered phase II clinical trials [17]. Experiments have shown that ingliforib has IC 50 values of 52, 352 and 150 nM for the liver, muscle and brain GP, respectively. However, further investigation is needed to elucidate the selective mechanisms of ingliforib for the three GP subtypes.
According to our previous report, when compound 1 acted on the three protein subtypes, a good distinction between the three GP subtypes was obtained compared to other screened compounds [3]. Compound 1 could be considered as a potential lead compound for the design of effective and selective GP inhibitors. Compound 1 mainly consists of two parts: the N-(3-(tert-butylcarbamoyl)-4-methoxyphenyl) moiety and the indole moiety. In 2000, indole-2-carboxamide analogs were first reported to act on the novel conformational site of HLGPa; they attracted extensive attention due to their unique structural properties [18]. Consequently, the indole unit was introduced into compound 1 for the exploration of specific inhibitors. However, compound 1 s specific mechanism of action with the three GP subtypes was also not clear. To explore the mechanism of action, in this study, compound 1 and ingliforib were used as the ligands to investigate the binding mode in the three GP isoforms by molecular docking, and the behavior of the complexes was also observed via molecular dynamics simulations. This work is expected to provide guidance for the further structural design and optimization of the highly selective compound.

Binding Model Studies by Molecular Docking
To investigate the associations between compound 1 and various GP subtypes, molecular docking was carried out. Initially, the GP protein structure was acquired from the protein database (refer to Materials and Methods) considering the following criteria: 1 human origin; 2 conformational resolution ≤ 2.5 Å, minimized to the greatest extent; 3 maximum possible conformational sequence coverage; 4 crystalline pH value as close to the normal physiological range of the human body as feasible. Ensuring that the above four conditions were satisfied to the greatest extent possible, the GP protein structures with PDB IDs of 5IKO (PYGB), 2ZB2 (PYGL) and 1Z8D (PYGM) were chosen for the molecular docking experiments in this study, following appropriate adjustments. Upon completion of molecular docking, the outcomes were examined based on the subsequent aspects.

Spatial Position and Orientation
After docking, the conformation with the lowest binding energy in GP was selected as the probable binding conformation. The spatial position and orientation of the ligand ingliforib and compound 1 in the three GP subtype receptors were observed. As shown in Figure 2A, ingliforib binds within the central cavity consisting of two identical subunits. In the complexes of the three subtypes, the indole, benzyl and dihydroxypyrrolidine moieties of the ligand molecule point towards three different pockets ( Figure 2C-E). Specifically, the indole moiety of ingliforib is deeply buried within pocket 1, a closed hydrophobic cavity [19]. The dihydroxypyrrolidine moiety extends towards a wide pocket 2, while the benzyl moiety points towards pocket 3. The ligand molecule adopts a "Y"-shaped conformation at the dimer interface. Comparison of the three subtype complexes using Pymol indicates that their spatial conformations are similar. The differences are mainly observed in dihydroxypyrrolidine, located in pocket 2, which has a different spatial orientation.
For compound 1, it also binds in a central cavity ( Figure 2B) [12,20]. Alignment of the three subtype complexes by Pymol reveals that the spatial location and orientation of the ligand in PYGL are different from those of the other two subtypes. However, after rotation, PYGL can be found to be in essentially the same position as the other two ligands. This leads to the inference that the ligand in PYGL may have entered another similar pocket in the dimer cavity with reverse symmetry. This may be reasonable as the site is located at the interface of a dimer composed of two identical subunits. Previous X-ray crystallography studies have demonstrated that two indole molecules of the GP inhibitor CP526423 bind identically within the solvent cavity at the dimer interface, within 6Å of each other [21]. In the three subtype complexes, the conformation of compound 1 in the pocket shows an approximate "Y" shape too ( Figure 2F-H), i.e., the indole, tert-butylamide group and methoxy point to three different pockets. The indole portion of compound 1 is positioned deep within pocket 1. The tert-butyl amide fragment points to a wide pocket 2. The methoxy fragment is in pocket 3. Comparison of the GP subtype complexes shows that the methoxy fragment of pocket 3 is not identical, with the ligand molecule in PYGB occupying the S1 region and the ligands in PYGL and PYGM occupying mainly the S2 region ( Figure 2I-K), which is also different from ingliforib. This may be due to the spatial or conformational constraints of the three GP isoforms. In addition, the ligand molecule in PYGB has room for further growth in pocket 3. In short, the ligand molecules exhibit a "Y"-shaped conformation at the dimer interface and there is a slight difference in pocket 3 occupied by methoxy. The difference in ligand conformation aroused our interest in their binding modes. Therefore, the binding mode was analyzed afterwards.

Binding Mode
Next, to illustrate the interaction mechanism, the binding modes of the ingliforib and compound 1 to the three isoforms were analyzed. The molecular docking results showed that the docking scores of ingliforib with PYGB, PYGL and PYGM were −10.782, −10.854 and −9.872 kcal/mol, respectively. Based on the molecular docking figures, although the binding modes of the three complexes were similar, some differences were also observed. In PYGB, the indole NH of ingliforib formed a hydrogen bond with Glu190, the indole formed a cation-π interaction with Lys191, the amide NH of the linker formed a hydrogen bond with Thr38′, the phenyl ring formed a π-π interaction with His57′, the carbonyl oxygen of the side chain amide formed a hydrogen bond with Ala192 and another hydroxyl group formed a hydrogen bond with Asn187 ( Figure 3A). The rich receptor-ligand inter-

Binding Mode
Next, to illustrate the interaction mechanism, the binding modes of the ingliforib and compound 1 to the three isoforms were analyzed. The molecular docking results showed that the docking scores of ingliforib with PYGB, PYGL and PYGM were −10.782, −10.854 and −9.872 kcal/mol, respectively. Based on the molecular docking figures, although the binding modes of the three complexes were similar, some differences were also observed. In PYGB, the indole NH of ingliforib formed a hydrogen bond with Glu190, the indole formed a cation-π interaction with Lys191, the amide NH of the linker formed a hydrogen bond with Thr38 , the phenyl ring formed a π-π interaction with His57 , the carbonyl oxygen of the side chain amide formed a hydrogen bond with Ala192 and another hydroxyl group formed a hydrogen bond with Asn187 ( Figure 3A). The rich receptor-ligand interactions endowed ingliforib with good inhibitory activity against PYGB. In PYGL, ingliforib also had abundant receptor-ligand interactions and additionally formed a hydrogen bond with Arg60, resulting in the best inhibitory activity against PYGL ( Figure 3B). However, in PYGM, ingliforib exhibited relatively poor inhibitory activity due to the lack of a cation-π interaction with Lys191 or a hydrogen bond interaction with Arg60, although it had similar receptor-ligand interactions as in PYGB ( Figure 3C). bonds were consistent with PYGB. Furthermore, the indole formed a cation-π interaction with Arg60, but its alkyl and π-alkyl interactions with some amino acid residues were weakened, and His57′ formed a π-π T-shaped interaction with the phenyl ring ( Figure  3F). This may be an important reason for the activity differences observed among the ligands in the three GP subtypes [3].
In comparing the complex of compound 1 with that of ingliforib, it was found that for the same subtype, the differences in their structures resulted in certain variations in the ways in which they interacted with amino acid residues. For example, in PYGL, the dihydroxy-pyrrolidine group of ingliforib in pocket 2 formed hydrogen bonds with Glu190′ and Ser192′, whereas this was not observed in compound 1. Additionally, the ways in which they interacted with the same amino acid residue could also be different. For instance, in PYGB, the tail benzene ring of compound 1 formed a π-alkyl interaction with ALA192. Such differences might explain why the two compounds exhibited distinct selectivity towards the three GP subtypes. In summary, the results suggested that the conformation of ligand compound 1 was stabilized in the three GP isoforms via hydrophobic interactions, hydrogen bonding and π-cation stacking, which is consistent with previous findings [16,20]. The indole ring as well as the two amide parts were considered to be the key parts of the activity. The docking scores of compound 1 with PYGB, PYGL and PYGM were −8.329, −7.244 and −8.266 kcal/mol, respectively. In PYGB, compound 1 formed hydrogen bonds with Glu190 and Thr38 , and a carbonyl oxygen of its side chain amide formed a hydrogen bond with Asn187 . Additionally, Lys191 formed a cation-π interaction with the indole ( Figure 3D). In PYGL, compound 1 lost the hydrogen bond with Asn187 and the carbon hydrogen bond interaction with Gly186 ( Figure 3E). In PYGM, except for the loss of a carbon-hydrogen bond interaction with the isobutylamide group, all other hydrogen bonds were consistent with PYGB. Furthermore, the indole formed a cation-π interaction with Arg60, but its alkyl and π-alkyl interactions with some amino acid residues were weakened, and His57 formed a π-π T-shaped interaction with the phenyl ring ( Figure 3F). This may be an important reason for the activity differences observed among the ligands in the three GP subtypes [3].
In comparing the complex of compound 1 with that of ingliforib, it was found that for the same subtype, the differences in their structures resulted in certain variations in the ways in which they interacted with amino acid residues. For example, in PYGL, the dihydroxy-pyrrolidine group of ingliforib in pocket 2 formed hydrogen bonds with Glu190 and Ser192 , whereas this was not observed in compound 1. Additionally, the ways in which they interacted with the same amino acid residue could also be different. For instance, in PYGB, the tail benzene ring of compound 1 formed a π-alkyl interaction with ALA192. Such differences might explain why the two compounds exhibited distinct selectivity towards the three GP subtypes. In summary, the results suggested that the conformation of ligand compound 1 was stabilized in the three GP isoforms via hydrophobic interactions, hydrogen bonding and π-cation stacking, which is consistent with previous findings [16,20]. The indole ring as well as the two amide parts were considered to be the key parts of the activity.

Molecular Dynamics Simulations
Based on the molecular docking results of ingliforib and compound 1 with the different subtypes of GP, molecular dynamics simulation was employed to explore the stability of their interactions and the key amino acid residues involved. The simulation was set for 50 ns with the initial structure as a reference, and the RMSD values were calculated to reflect the stability of the receptor-ligand complex conformation during the simulation. The hydrogen bond occupancy rate was calculated based on the stable trajectory to assess the stability of the hydrogen bond interactions. Additionally, the binding free energy was calculated, and energy decomposition was performed to determine which amino acid residues contributed significantly to the binding of compound 1 and ingliforib.
In Figure 4A,

Molecular Dynamics Simulations
Based on the molecular docking results of ingliforib and compound 1 with the different subtypes of GP, molecular dynamics simulation was employed to explore the stability of their interactions and the key amino acid residues involved. The simulation was set for 50 ns with the initial structure as a reference, and the RMSD values were calculated to reflect the stability of the receptor-ligand complex conformation during the simulation. The hydrogen bond occupancy rate was calculated based on the stable trajectory to assess the stability of the hydrogen bond interactions. Additionally, the binding free energy was calculated, and energy decomposition was performed to determine which amino acid residues contributed significantly to the binding of compound 1 and ingliforib.
In Figure 4A,    Figure 4D). Similarly, the PYGL-compound 1 complex also remained stable ( Figure 4E). After 10 ns, the protein predominantly fluctuated around 2.0 Å, while compound 1 fluctuated between 1.25 and 2.25 Å. For the PYGM-compound 1 complex ( Figure 4F), equilibrium was reached after 10 ns, and the protein fluctuated mainly between 2.0 and 2.5 Å, while compound 1 fluctuated between 2.0 and 3.0 Å. Overall, the fluctuation range of the compound 1-GP complex was narrower than that of the ingliforib-protein complex.
Considering the memory required to calculate hydrogen bond occupancy rates, the stable molecular dynamics simulation trajectories between 45 and 50 ns were selected to calculate the hydrogen bond occupancy rates. Comparing the hydrogen bond occupancy rates between the three receptor-ligand complexes (Table 1), it can be observed that both PYGB-compound 1 and PYGL-compound 1 maintained a stable binding conformation with the molecular docking and formed a stable hydrogen bond with Glu190. Additionally, in PYGL-compound 1, the compound formed stable hydrogen bond interactions with Thr38 and His57 . However, the hydrogen bond between the compound and Glu190 was not maintained in PYGM-compound 1, but compound 1 formed a stable hydrogen bond with His57, with an occupancy rate of 99.85%, indicating that the conformation of the compound may have undergone some flipping. In PYGB-compound 1 and PYGMcompound 1, the intramolecular hydrogen bonds of compound 1 were relatively stable, with hydrogen bond occupancy rates of 59.03% and 91.72%, respectively, suggesting that the tert-butyl carbamoyl amide may have undergone some flipping. Although there were also intramolecular hydrogen bonds in PYGL-compound 1, the occupancy rate was only 11.74%. These trajectories reflected that compound 1 could obtain relatively stable binding conformations after binding to the three subtypes of GP. The stable molecular dynamics simulation trajectories between 45 and 50 ns of the ingliforib-protein complexes were selected to calculate the hydrogen bond occupancy and binding free energy. The results of the hydrogen bond occupancy are presented in Table 1. In PYGB-ingliforib, stable hydrogen bonds were formed between Glu190 and Lys191 and ingliforib, with hydrogen bond occupancy of 96.86% for Glu190. Ingliforib also showed some intramolecular hydrogen bonding interactions. In PYGL-ingliforib, ingliforib formed stable hydrogen bonds with Thr38 , Glu190 and Ser192, with high hydrogen bond occupancies. The hydrogen bond occupancy with Ser192 reached 76.07%. In PYGMingliforib, ingliforib primarily formed stable hydrogen bonds with Glu190, Thr38 and Lys191, along with some intramolecular hydrogen bonds. Comparing the three systems, it was found that the stable hydrogen bond interaction with Glu190 was present in all three systems. The hydrogen bond interaction with Lys191 was present in PYGB and PYGM, but not in PYGL. The hydrogen bond interaction with Thr38 was absent in PYGB, suggesting that this hydrogen bond may not be necessary for the inhibitory activity of the compound against PYGB. The stable hydrogen bond interaction with Ser192 was only present in PYGL, indicating that it may enhance the selectivity of the compound for PYGL.

Calculation of Binding Free Energy Based on MM-PBSA
On the basis of the MD simulation, the MM-PBSA method was employed for the calculation of the binding free energies of compound 1 with the three subtypes of GP (Table 2), which were found to be −38.26 ± 1.61, −37.61 ± 2.08 and −35.10 ± 2.57 kcal/mol for compound 1 with PYGB, PYGM and PYGL, respectively. This indicates that compound 1 has the strongest binding affinity with PYGB, followed by PYGM and then PYGL (IC 50 of PYGB, PYGL and PYGM, which are 0.11 ± 0.01 µM, 0.35 ± 0.02 µM and 0.93 ± 0.01 µM, respectively) [3]. The strong correlation observed between the computed interaction free energies and the inhibitory activity obtained experimentally indicates that the binding conformations of the inhibitors are reasonable. Then, the binding free energy was further decomposed to identify the amino acid residues that contribute significantly to compound 1 s binding between 45 and 50 ns, focusing on residues whose contribution values are greater than −0.01 kcal/mol, with values exceeding −1.00 kcal/mol considered important. In the PYGB-compound 1 complex ( Figure 5A), Arg60, Glu190, Thr38 , Val40 , Phe53 and His57 were found to make a significant contribution to compound 1 s binding with PYGB. In the PYGL-compound 1 complex ( Figure 5B), Arg60, Val64, Glu190, Lys191, Thr38 , Val40 , His57 and Pro188 were identified as important residues for compound 1 s binding with PYGL. In the PYGM-compound 1 complex ( Figure 5C), His57, Arg60, Trp189, Thr38 , Leu39 , Phe53 and Pro188 were found to have significant contributions to compound 1 s binding with PYGM. These results indicate that amino acid residues on both monomers of the protein dimer contribute significantly to compound 1 s binding. In comparison, Pro188 had a relatively small contribution value in the PYGB-compound 1 system, but it was an important residue in both the PYGL-compound 1 and PYGM-compound 1 systems. Additionally, Glu190 played an important role in the binding of compound 1 to both PYGB and PYGL, while it was not an important residue in the PYGM-compound 1 system, despite the substantial contribution of Trp189. Therefore, we concluded that maintaining interactions with Arg60, Glu190 and Thr38 , as well as interacting with Leu39 , Val40 , Phe53 and His57 on the other monomer, while reducing interactions with Trp189 and Pro188 , can lead to improved selectivity for PYGB.

Pharmacophore Modeling
To provide better theoretical guidance for the design of new inhibitors, here, we provide some references for the design of new compounds. Based on the experimental results, maintaining the interactions between Arg60, Glu190 and Lys191, while forming interactions with Thr38′ and Leu39′ of the other half of the protein monomer, are some key factors required for the inhibitory activity of the compound against GP. For PYGL, given that the amino acid residue 192 is Ala in both PYGB and PYGM but Ser in PYGL, the formation of a strong interaction with Ser192 would facilitate the selectivity of the compound for PYGL. Conversely, strengthening the interactions between the compound and Ala192 is also ben- As demonstrated in Table 2, the strongest binding was between ingliforib and PYGL, followed by PYGB, and the weakest was with PYGM. The energy decomposition analysis revealed the importance of certain amino acid residues in ingliforib binding. For the PYGB-ingliforib complex ( Figure 5D), Arg60, Glu190, Lys191, Ala192, Thr38 , Leu39 and Phe53 were found to make significant contributions to the binding. For the PYGL-ingliforib complex ( Figure 5E), Arg60, Glu190, Lys191, Ser192, Arg193, Thr38 , Leu39 and Val40 were identified as key contributors. In the PYGM-ingliforib complex ( Figure 5F), Arg60, Lys191, Ala192, Thr38 and Leu39 were crucial for the binding. It can be seen that the amino acid residues discovered through energy decomposition make significant contributions to the binding of the compound. Comparing the three proteins, maintaining mutual interactions between Arg60, Glu190 and Lys191, while interacting with Thr38 and Leu39 of the other half of the protein, is necessary for the inhibitory activity of the compound against all three proteins. Moreover, hydrogen bonding with Ser192 can increase the selectivity of the compound for PYGL, since Ala192 is present in both PYGB and PYGM. Lys191 also plays a critical role in the structure of the compound-PYGL complex, with its contribution being twice that in PYGB and PYGM, suggesting that enhancing the interaction between Lys191 and the compound would improve its selectivity against PYGL. In addition, Arg193 is also important for the structure of the complex, while its contribution in PYGB and PYGM is not significant; thus, increasing the interaction between the compound and Arg193 would also be advantageous for its selectivity against PYGL.

Pharmacophore Modeling
To provide better theoretical guidance for the design of new inhibitors, here, we provide some references for the design of new compounds. Based on the experimental results, maintaining the interactions between Arg60, Glu190 and Lys191, while forming interactions with Thr38 and Leu39 of the other half of the protein monomer, are some key factors required for the inhibitory activity of the compound against GP. For PYGL, given that the amino acid residue 192 is Ala in both PYGB and PYGM but Ser in PYGL, the formation of a strong interaction with Ser192 would facilitate the selectivity of the compound for PYGL. Conversely, strengthening the interactions between the compound and Ala192 is also beneficial to improve its selectivity against PYGB and/or PYGM. Therefore, introducing corresponding functional groups or structures into the compound to enhance the interactions with Ser192 could be considered. In PYGB, a strategy may be to enhance the interaction of the compound with the amino acid residues around the S1 region that make up pocket 3 through structural modifications.
In terms of compound structure, for the indole portion of pocket 1, the indole unit is constrained in a relatively narrow sub-pocket (pocket 1). Although the scope for the growth of the indole unit is limited, it may be a better strategy to try to substitute it according to the bioelectronic isomer rule. For the methoxy fragment in pocket 3, introducing electronwithdrawing groups to enhance its interaction may be effective due to the electronegative region at the bottom of the pocket ( Figure 6A), which may not be valid for PYGL and PYGM ( Figure 6B,C). As shown in Figure 6D-F, attempting to introduce more hydrophilic groups while maintaining the original hydrogen bonds in the ligand portion of pocket 2 may be effective in further improving the overall activity (forming surface colors ranging from hydrophilic blue to hydrophobic brown based on the hydrophobicity of the receptor residues). The above findings may offer some contribution to the design of GP inhibitors specifically against GP isoenzymes.
In summary, in the GP-compound 1 or GP-ingliforib complexes, there was some variability in the amino acids involved and their contributions. Therefore, modifications should be made to the compound according to specific situations for different protein systems, to enhance or weaken its interactions with certain amino acid residues. Meanwhile, the spatial structure of the compound should be further optimized to better fit the active site of the target protein and increase its non-covalent interactions with the target protein.
Molecules 2023, 28, x FOR PEER REVIEW 12 of 16 active site of the target protein and increase its non-covalent interactions with the target protein.

Molecular Docking
Firstly, molecular docking was performed using the Glide module of Schrö-dinger2021. The PDB files 5IKO [13,22], 2ZB2 [23,24] and 1Z8D [25,26] were used for PYGB, PYGL and PYGM, respectively. The dimer was directly downloaded from the Uni-Prot database. For Tetramer, delete the redundant monomer to obtain the dimer. Although only 2ZB2 has a complex crystal structure of PYGL with a small molecule, information about the binding sites of small molecules in PYGB and PYGM can be obtained by superimposing them based on the sequence homology. According to the crystal structures, small molecules bind to the middle site of the protein dimer. Consequently, a grid file of the binding site was generated using Glide Grid Generation. The Ligpre module was utilized to process small molecules with the OPLS4 force field setting, hydrogenation, ionization and energy optimization. Finally, molecular docking was conducted using the Glide XP (extra precision) mode. The docking was semi-flexible, i.e., the ligand conformation was flexible while the pocket was fixed. The grid file for the binding site was defined by a box centered on the crystal ligand with a similar size. No constraints were included during grid generation. The Glide XP docking score was used to rank the ligand

Molecular Docking
Firstly, molecular docking was performed using the Glide module of Schrödinger2021. The PDB files 5IKO [13,22], 2ZB2 [23,24] and 1Z8D [25,26] were used for PYGB, PYGL and PYGM, respectively. The dimer was directly downloaded from the UniProt database. For Tetramer, delete the redundant monomer to obtain the dimer. Although only 2ZB2 has a complex crystal structure of PYGL with a small molecule, information about the binding sites of small molecules in PYGB and PYGM can be obtained by superimposing them based on the sequence homology. According to the crystal structures, small molecules bind to the middle site of the protein dimer. Consequently, a grid file of the binding site was generated using Glide Grid Generation. The Ligpre module was utilized to process small molecules with the OPLS4 force field setting, hydrogenation, ionization and energy optimization. Finally, molecular docking was conducted using the Glide XP (extra precision) mode. The docking was semi-flexible, i.e., the ligand conformation was flexible while the pocket was fixed. The grid file for the binding site was defined by a box centered on the crystal ligand with a similar size. No constraints were included during grid generation. The Glide XP docking score was used to rank the ligand positions. Protonation of the residues was automatically determined under the neutral condition using the Protein Preparation Wizard workflow. Results were visualized using Pymol and Discovery Studio 4.5 Visualizer.

Molecular Dynamics Simulation
The molecular dynamics simulation utilized the GROMACS 2019.6 package [27], with the CHARMM 36 force field [28][29][30][31] the TIP4P water model and suitable ion parameters [32,33] as well as ligand force field parameters generated with the SwissParam webserver [34]. The protein-ligand complex was first energy-minimized in a vacuum, followed by immersion in a dodecahedral periodic box with three-dimensional boundary conditions, at a minimum distance of 12 Å from the box frame. To simulate a physiological environment, the simulation box was filled with TIP4P model water molecules, and sodium and chloride ions were added to balance the system charge at an ion concentration of 150 mM. After solvation, the protein's energy was minimized again with positional restraints applied to its backbone atoms, while allowing the solvent to freely diffuse. Equilibration of the system continued with a 1 ns NVT run at 300 K (using a heat bath based on the v-rescale algorithm with EnerPres correction), followed by a 1 ns NPT run at 1 atm (using the Parrinello Rahman the pressure controller), both still with positional restraints. NPT collection runs with no positional restraints were 50 ns long. The LinCS method constrained all bonds. Particle Mesh Ewald (PME) was used for long-range electrostatic interactions, with a cutoff value of 10 Å. The program provided with GROMACS 2019.6 was used to determine the root mean square deviation (RMSD). The complexes for each GP subtype in complex with the compounds obtained by molecular docking were used as the initial conformations for molecular dynamic simulations and also used as the reference structures to calculate the RMSD values. The results were visualized using Pymol and Discovery Studio 4.5 Visualizer.

Calculation of Binding Free Energy Based on MM-PBSA Method
The estimation of the interaction free energy was performed using the molecular mechanics Poisson Boltzmann surface area (MM-PBSA) approach. The binding free energy was computed through the MM-PBSA method with g_mmpbsa [35,36]. The system's enthalpy was evaluated using the molecular mechanics (mm) approach, and the contributions of the polar and non-polar components of the solvent effect to the free energy were assessed by resolving the Poisson Boltzmann (PB) equation and computing the molecular surface area (SA). The vacuum potential energy, polar solvation energy and non-polar solvation energy were ascertained using g_mmpbsa computation. The Python script included in the g_mmpbsa package was employed to determine the mean binding energy and standard deviation.

Conclusions
In this study, the interaction mechanisms between ingliforib and compound 1 and PYGB, PYGL and PYGM were elucidated through molecular docking and molecular dynamics simulations. Observations revealed that certain dissimilarities existed in the three protein isoforms with regard to their spatial conformations, and corresponding differences were also evident in the binding modes of the ligands. These conformations were stabilized by interactions such as hydrophobic interactions, hydrogen bonds and π-cation stacking. The rich receptor-ligand interactions endowed ingliforib and compound 1 with good specific inhibitory activity against PYGL and PYGB, respectively. Subsequently, these results were validated through molecular dynamics simulations. The RMSD values of the protein-ligand complexes during 50 ns of molecular dynamics simulation remained stable within a specific range. The hydrogen bond occupancy and binding free energy were further calculated using the MM-PBSA method, and some key amino acid residues in each system that significantly contributed to ligand binding were identified. These findings provide a theoretical basis for the optimization or design of new GP inhibitors to achieve better inhibitory activity against specific target proteins.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/molecules28134909/s1, Figure S1: Predicted binding mode of molecular docking between compound 1 and PYGB; Figure S2: Predicted binding mode of molecular docking between compound 1 and PYGL; Figure S3: Predicted binding mode of molecular docking between compound 1 and PYGM; Figure S4: Predicted binding mode of molecular docking between Ingliforib and PYGB; Figure S5: Predicted binding mode of molecular docking between Ingliforib and PYGL; Figure S6: Predicted binding mode of molecular docking between Ingliforib and PYGM; Figure S7: Atomic number of compound 1; Figure S8: Atomic number of ingliforib.