In silico panning for a non-competitive peptide inhibitor

Background Peptide ligands have tremendous therapeutic potential as efficacious drugs. Currently, more than 40 peptides are available in the market for a drug. However, since costly and time-consuming synthesis procedures represent a problem for high-throughput screening, novel procedures to reduce the time and labor involved in screening peptide ligands are required. We propose the novel approach of 'in silico panning' which consists of a two-stage screening, involving affinity selection by docking simulation and evolution of the peptide ligand using genetic algorithms (GAs). In silico panning was successfully applied to the selection of peptide inhibitor for water-soluble quinoprotein glucose dehydrogenase (PQQGDH). Results The evolution of peptide ligands for a target enzyme was achieved by combining a docking simulation with evolution of the peptide ligand using genetic algorithms (GAs), which mimic Darwinian evolution. Designation of the target area as next to the substrate-binding site of the enzyme in the docking simulation enabled the selection of a non-competitive inhibitor. In all, four rounds of selection were carried out on the computer; the distribution of the docking energy decreased gradually for each generation and improvements in the docking energy were observed over the four rounds of selection. One of the top three selected peptides with the lowest docking energy, 'SERG' showed an inhibitory effect with Ki value of 20 μM. PQQGDH activity, in terms of the Vmax value, was 3-fold lower than that of the wild-type enzyme in the presence of this peptide. The mechanism of the SERG blockage of the enzyme was identified as non-competitive inhibition. We confirmed the specific binding of the peptide, and its equilibrium dissociation constant (KD) value was calculated as 60 μM by surface plasmon resonance (SPR) analysis. Conclusion We demonstrate an effective methodology of in silico panning for the selection of a non-competitive peptide inhibitor from small virtual peptide library. This study is the first to demonstrate the usefulness of in silico evolution using experimental data. Our study highlights the usefulness of this strategy for structure-based screening of enzyme inhibitors.


Background
According to market research, the potential of peptide therapeutics has recently intensified [1][2][3]. Worldwide, there are more than 40 marketed peptides, with about 270 peptides in clinical phase testing, and about 400 peptides in advanced preclinical phases [1]. Natural peptides, such as insulin, vancomycin, oxytocin, and cyclosporine, and synthetically produced peptides, such as Fuzeon (enfuvirtide) and Integrilin (eptifibatide), are among the approved peptide-based drugs. Compared to low-molecular-mass chemical drugs, peptide drugs offer several advantages, such as high specificity, minimization of drug-drug interactions, lower accumulation in tissues, lower toxicity, and biological diversity. However, peptides also have some disadvantages, which include low oral bioavailability, lower stability, higher risk of immunogenic effects, difficulties associated with delivery due to rapid clearance from the body, and costly synthesis. Recently, several novel and interesting approaches to deliver protein-based drugs through the skin have been reported [4]. Since peptides require costly synthesis, highthroughput screening (HTS) of numerous peptides from combinatorial libraries is inefficient. Therefore, novel procedures that require less effort for the screening of peptide ligands are required. From this point of view, structurebased computational drug design is an effective methodology. Recent advances in protein structure determination, achieved either through X-ray crystallography or NMR, are providing informative data related to the design of useful drugs based on these proteins. The identification of the binding sites on these newly determined protein structures have led to the development of a variety of docking strategies. There are numerous reports of drug discovery from small molecule ligand libraries [5,6], although it is difficult to calculate the docking energies of all the peptide sequence patterns, as they show enormous diversity. Therefore, we focused on the use the genetic algorithms (GAs) to reduce the redundancy of the selection procedure.
GAs represent a class of algorithms that mimic some of the major characteristics of Darwinian evolution [7,8]. GAs are based on the process of genetic evolution observed in biological systems, in which three successive operations, selection, crossover, and mutation, are performed on a set of strings. GAs provide an effective means of exploring the conformational space of flexible molecules. GAs also provide an effective approach to protein folding [9], identification of the biomolecular conformation space [10], docking methodology [11,12], optimization of lead compounds [8,13], chemical evolution of combinatorial chemistry [14], and identification of receptor-ligand binding sites [15].
We have previously reported the application of GAs to select a peptide inhibitor [16], an α-helix-forming peptide [17], and a DNA aptamer with higher-order structure [18][19][20]. From the result of those studies, it is clear that GAs are useful for the efficient selection of molecules that have a desired property or function, since we can reduce the number of rounds of evaluation. In the present study, we have focused on the application of GAs for effective peptide ligand selection from a docking simulation. Belda et al. [21] have also reported a combination of computational docking and combinatorial experimental screening but have not provided experimental data. We propose an effective approach to derive peptide ligands, which we call 'in silico panning'. By combining the docking study and GAs, we are able to identify promising peptide ligands from a small virtual peptide library with less effort ( Figure  1).
To demonstrate in silico panning, we chose the water-soluble quinoprotein glucose dehydrogenase (PQQGDH) from Acinetobacter calcoaceticus as the target protein.
Mainly due to its high catalytic activity and non-dependence on oxygen as an electron acceptor, PQQGDH has replace glucose oxidase (GOD) as the major enzyme used in glucose sensor systems [22]. The only aspects in which PQQGDH is inferior to GOD are substrate specificity and operational stability. PQQGDH shows high activity not only for glucose, but also for disaccharides, such as lactose and maltose, as substrates. Improvements in substrate specificity are expected to lead to the development of more-sensitive glucose sensors. Therefore, PQQGDH engineering has been carried out in our research group [22,23]. We have already reported mutants of PQQGDH (Glu277Lys, Asn452Thr, Asp167Glu, and Asp167Glu/ Asn452Thr) that show improved substrate specificities [24,25].
We have also been proposing quaternary structure engineering of proteins, which is the control of protein function using an artificial subunit. In nature, several enzymes are composed of subunits that form an active quaternary structure, which endows a higher order of function. The strategy of quaternary structure engineering is to mimic the native enzyme using a peptide ligand as an artificial subunit with a novel function. From a phage display peptide library, we previously identified a peptide with seven amino acids (Thr-Thr-Ala-Thr-Glu-Tyr-Ser) that narrowed the substrate specificity of PQQGDH without significant loss of enzyme activity [26]. In the present study, in silico peptide evolution was performed to investigate a new methodology for protein modification. Our ultimate objective is the selection of a non-competitive peptide ligand that does not interrupt glucose binding but inhibits the interaction between disaccharides, such as maltose and lactose. The disaccharides appear to access a large pocket structure next to the catalytic site of PQQGDH, which comprises Arg148, Arg406, Arg408, Arg45, and Lys28 ( Figure 2). We defined this large pocket as the docking field. From size of this pocket, tetra peptide was designed to be the most suitable size as initial tetra peptide library. Then, the peptide ligand was selected through a combination of the docking simulation and GAs on the basis of binding indicators.
We chose the Molecular Operating Environment of docking (MOE-Dock) system with simulated annealing method as a starter system [27,28]. Simulated annealing is a global optimization technique that is based of the Monte Carlo method. It explores various states of a configuration space by generating small random changes in the current state and then accepting or rejecting each new state according to the Metropolis criterion [29].

Docking energy transition in combination with in silico evolution
One of the key steps in docking is how to define the userspecified three-dimensional docking box (3D docking box). To select the non-competitive peptide ligand, we designated the large pocket next to the glucose-binding site ( Figure 2). Then the peptide ligand was selected by in silico panning; a combination of the docking simulation and GAs on the basis of binding indicators. In all, four rounds of selection were carried out from the beginning of ten virtual peptide library on the computer. The docking energy (fitness) results for each selection round are listed in Table 1. In MOE-Dock, electrostatic energies and Van der Waals energies between the target and ligand are calculated by simulated annealing. The docking energy values were calculated as the sum of the electrostatic, Van der A schematic diagram of the in silico peptide evolution system Figure 1 A schematic diagram of the in silico peptide evolution system. The docking calculation program and genetic algorithms (GAs) were combined to evolve the peptide ligand on the computer. For the docking program, we used the MOE-Dock software, which is based on the simulated annealing method. GAs were used to evolve the peptide to produce the next generation. Four rounds of peptide evolution were performed in this study.
Waals energies and the flexibility of the ligand itself. Low docking energy indicates high binding ability. The distributions of the docking energies in each generation are illustrated in Figure 3; each docking energy decreased gradually as peptide evolution progressed. After 4 th round of peptide evolution, the docking energies of most of the peptides were lower than those in the 1 st -3 rd round. Peptide evolution and the direction of the selection proceeded successfully using the GA.

Homologies of highly ranked peptide ligands from the docking simulation
In comparisons of the homologies of high-scoring peptides, the sequences appeared to converge ( Table 2). The peptide homologies converged almost completely after the 4 th round of peptide evolution. A negatively charged residue was revealed at the third position with a probabil-ity of 80%. Furthermore, 60% of the peptides contained positively charged residues at the fourth position. On the other hand, lysine residues did not appear at the fourth position during the four rounds of selection due to deflection of evolution. For this reason, we designed additional peptides, and re-calculated their docking energies. Finally, the top three peptides and negative controls were synthesized for further analysis. Figure 4 displays an image of PQQGDH interacting with peptide GEKD, which was derived by MOE-Dock.

Effects of synthetic peptide ligands on PQQGDH activity
As shown in Figure 5A, 20 µM of peptide SERG showed the strongest inhibitory effect on the target enzyme. The same concentration of GEKD or GERD also showed higher inhibition than the control peptide DDDD. The velocity at maximal concentrations of substrate(V max ) Defining a docking box next to the glucose-binding site Figure 2 Defining a docking box next to the glucose-binding site. The 3D structure of water-soluble quinoprotein glucose dehydrogenase from Acinetobacter calcoaceticus (PQQGDH: PDB entry code 1CQ1) is shown in gray. PQQ and glucose are represented in yellow and green, respectively. The large pocket next to the glucose-binding site of PQQGDH comprises Arg148, Arg406, Arg408, Arg45 (red) and Lys28 (orange). The docking box was set as the target area for this pocket. The docking box is shown as a yellow rectangular solid. The size of the docking box was fixed as a 37 Å The distribution of docking energies in each generation Figure 3 The distribution of docking energies in each generation. The ranks of the peptide ligands in each generation are shown on the X-and Y-axes. The Z-axis shows the docking energies, which were calculated as follows: docking total energy = electrostatic energies + Van der Waals energies + the energy of the (flexible) ligand. Low docking energy means high binding ability of the ligand.  Figure 5B), SERG was deemed to show non-competitive inhibition of the substrate. Non-competitive inhibitors do not affect the combination of the substrate with the enzyme, but it does affect the velocity. In a pure non-competitive system, the substrate has an identical affinity for both the Enzyme(E)-Inhibitor(I) complex and enzyme. Unlike the Enzyme(E)-Substrate(S) complex, the E-I-S complex cannot convert the substrate to product. Therefore, the K m value is unchanged while V max is lowered. The data show that SERG decreases the rate constant for product formation, while K m shows a constant value. It follows that SERG is not a glucose competitor. On the other hand, GERD and GEKD showed 'mixed-type inhibition' (data not shown).
The peptide SERG showed dose-dependent inhibition ( Figure 5C) and its enzyme inhibition constant(Ki) value was calculated as 20 µM from the Dixon plot ( Table 3). All of the Ki values were calculated using the KaleidaGraph software (Synergy Software, Boston, MA, USA). Although the selected peptides showed potent inhibition, they did not show significantly decreased substrate specificities for disaccharides (data not shown).

Binding parameters of selected peptides by surface plasmon resonance (SPR)
During the SPR experiments, the binding assay for the biotinylated peptides (GEKD and SERG) and PQQGDH was processed by the BIAcoreX instrument (BIAcore AB, Uppsala, Sweden). The equilibrium dissociation constant (K D ) value used to evaluate the enzyme-peptide binding affinity was determined by Scatchard plot (Table 3). For each trial, the signal was corrected for the control surface response. We were able to confirm peptide binding to the enzyme, and the K D values of peptides GEKD and SERG were calculated as 155 µM and 60 µM, respectively.

Discussion
In this paper, we introduce the idea of in silico panning, which is a novel strategy to select peptide inhibitors by combining a docking simulation and GA. Setting the docking field next to the substrate binding site, we were Peptide GEKD docked in PQQGDH Figure 4 Peptide GEKD docked in PQQGDH. GDH, PQQ, and glucose are represented in purple, green, and pink, respectively. The GEKD peptide is displayed in red. The docking box is shown as a yellow rectangular solid. GEDG -99 10 GGGG -98 The homologies of selected peptides ranked in the top 10 after the 4 th round of evolution were listed to the left. Additional peptides that contained a lysine (K) residue at the third position were also listed on the middle. The final rank and the synthetic peptide sequences were shown on the right. The DDDD peptide was used as a negative control.
PQQGDH activity for glucose in the presence of individual peptides  able to obtain a non-competitive peptide inhibitor from a small virtual peptide library. In this study, the non-competitive peptide inhibitor was successfully and with less effort selected from a small virtual peptide library using in silico panning. This methodology can be used to screen an allosteric binding peptide inhibitor at an early stage.
Information derived from X-ray crystallography, NMR spectroscopy, and homology modelling greatly facilitates the rational design of selective and potent inhibitors. In addition, rapid identification of lead molecules and optimization of inhibitors are facilitated by large combinatorial libraries and high-throughput screening. However, successful virtual screening of chemical libraries in the drug discovery process requires a sufficiently large and chemically diverse library of compounds. In addition, the selection of promising peptide inhibitors from such large libraries involves significant cost and effort [1]. Therefore, we focused on the use of the GAs to reduce redundancy during the selection procedure.
In in silico panning, all of the selection processes are carried out by the computer. First, we calculate the affinity of the peptide ligand from docking simulation, and then evolve the sequence by GA, which mimic Darwinian evolution. In the present study, the initial group involved ten virtual tetrapeptides designed randomly from seven different amino acid residues (Arg, Lys, Asp, Glu, Ser, Pro, and Gly). The Arg, Lys, Asp, Glu, and Ser residues were mainly chosen to form the electrostatic and hydrogen interaction. In addition, hydrophobic Pro residues were appeared frequently in a phage display peptide library in a previous study [26], so we chose this residue. Gly was chosen to increase the flexibility of the ligand. The large pocket next to the glucose binding site of PQQGDH is mainly composed of hydrophilic residues, including Arg148, Arg406, Arg408, Arg45, and Lys28. In this study, to simplify the experimental method and concentrate on the hydrophilic interactions, we chose mainly hydrophilic residues for designing the initial peptide library. In the near future, a more complex library, containing several hydrophobic residues, will be used. Totally, four rounds of selection were carried out, and 45 types of tetrapeptide encompassing seven different amino acid residues were evolved by the computer. Since this covered only 18.7% of the 2401 possible tetrapeptide combinations, we conclude that the peptide ligands were efficiently selected from this small library by in silico panning. This is the first study to demonstrate the usefulness of in silico evolution using experimental data.
We have reported previously the identification from a phage display peptide library of a 7-mer peptide that causes PQQGDH substrate specificity towards disaccharides to decrease significantly [26]. We expected that some of the evolved tetrapeptides would improve the substrate specificity of PQQGDH. Although we were not able to obtain a peptide that improved substrate specificity, we obtained a non-competitive peptide inhibitor of the target protein. The selected peptide inhibits the enzyme activity not only for glucose, but also for disaccharides (data not shown). Thus, the substrate specificity does not change in the presence of the peptide ligand. To overcome these problems, we have to select peptide ligands that bind to the target pocket without inhibiting glucose binding. It is possible to calculate the binding of glucose or disaccharides to target area after derivation of a candidate peptide using the docking simulation. Thus, we may be able to choose a ligand that decreases only the reactivities for disaccharides.
In this study, we present a valuable method for the selection of a non-competitive peptide inhibitor. For target enzyme that are highly homologous, particularly at the catalytic site, target-specific inhibition is possible through the use of an allosteric or non-competitive inhibitor, which does not bind to the catalytic site but has inhibitory activity. Thus, a non-competitive inhibitor sometimes became a specific inhibitor for protein having high homogeny [30]. Although an RNA aptamer and peptide inhibitor with non-competitive and specific inhibitory activities have been selected from a pooled random sequence library [30,31], a method that allows one to choose a non-competitive inhibitor de novo is more efficient. In our study, by setting the enzyme pocket, which appears to be an allosteric binding site, as the docking area, the non-competitive peptide inhibitor was obtained as expected. The top three selected peptides that were evolved by binding index in the docking study showed much greater inhibition compared to the negative control, DDDD. The most effective peptide, SERG showed potent inhibition with a K i value of 20 µM and K D value of 60 µM.
In MOE-Dock, the docking energies were calculated as the sum of three energies, electrostatic energies, van der Waals energies and energy of the (flexible) ligand. Analysis of the individual components of the calculated binding energies (Table 4) shows that the top three peptides have lower electrostatic energies and energies of flexibility than the negative control peptide DDDD. All of the ten peptides that were expected to have high affinity also showed low electrostatic energies and energy of flexibility (data not shown). We initially expected that the peptide DDDD might show high affinity since the pocket next to the glucose-binding site in PQQGDH contains five positively charged residues. However, the docking energy of DDDD was ranked 38 th in the docking study ( Table 1). The peptide DDDD showed weak electrostatic interaction with PQQGDH due to the low flexibility that comes from the intramolecular ionic repulsion (Table 4). In fact, the peptide DDDD did not show remarkable inhibitory effect, suggesting that flexibility of a peptide inhibitor may play an important role in the interaction with the target molecule.
In future studies, the use of longer sequences that form higher-order structures will generate more specific peptide inhibitors from in silico panning. This method has strong potential to become a useful tool for structure-based noncompetitive inhibitor screening.

Conclusion
We have demonstrated the potential of in silico panning for selecting non-competitive peptide inhibitors in a more efficient manner. By choosing a target region next to catalytic site in the docking study and then evolving the peptide ligand on the basis of binding indicators, we succeeded in obtaining a non-competitive peptide inhibitor. The most effective peptide showed potent inhibition with a K i value of 20 µM and K D value of 60 µM. Our in silico panning approach should become a useful tool for screening structure-based enzyme inhibitors. This methodology has excellent potential for the screening of noncompetitive peptide ligands for allosteric binding sites of an enzyme in the early stages.

Design of the initial peptide library
Ten tetrapeptides encompassing seven different amino acid residues were identified in the first generation. The large pocket next to the catalytic site of PQQGDH, which is mainly composed of the conserved residues of Arg148, Arg406, Arg408, Arg45, and Lys28, was designated as the docking field. Arg, Lys, Asp, Glu, and Ser residues were chosen to form electrostatic and hydrogen interactions with the residues in the large pocket. The hydrophobic amino acid proline appeared frequently in the phage display peptide library examined in a previous study [26], so we also chose this residue. Gly was added to increase the flexibility of the ligand.

Calculation of docking energies of peptide ligands for the enzyme catalytic site
The Molecular Operating Environment of docking (MOE-Dock; [27,28]) was used to calculate the docking energies between peptide ligands and enzyme catalytic site. MOE-Dock is used to search for favorable binding configurations between a small ligand and a macromolecular target. The peptides were drawn by the MOE software and stored in the database.
The PDB structure of GDH (PDB entry code 1CQ1) and stored peptide were imported at the start of the docking program. Since not all X-ray crystallographic files contain hydrogen atoms, we added them to the protein using the MOE modelling suite before carrying out the docking studies. Minimizing contacts for hydrogen, the structures were subjected to an AMBER94 energy minimization protocol. The docking energy calculation was carried out within a user-specified three-dimensional docking box (3D docking box) using the simulated annealing method under the MMFF94 force field. The energy grids for docking were generated as grid-based potential fields by the MOE-Dock program, to reduce the calculation time.
We selected the large pocket next to the catalytic site of PQQGDH, which includes Arg148, Arg406, Arg408, Arg45, and Lys28, as the target area in the 3D docking box. After importing the peptide to the 3D docking box, we initiated the calculation of docking energies. Each docking energy value was calculated as the sum value of the electrostatic, Van der Waals, and flexibility energies. The interaction energy was calculated using the electrostatic and Van der Waals potential fields sampled on a grid overlaying the 3D docking box. The 3D docking box was interpolated at the atom positions by tri-linear interpolation. The Van der Waals parameters were taken from the currently active force field. The electrostatic field was calculated in the Coulombic manner using the current dielectric. MOE-Dock performed 25 independent docking runs, and wrote the resulting conformations and their energies to a molecular database file. The lowest docking energy conformation for each peptide was chosen for GA evolution. All of the docking simulations were carried out on a personal computer.

Evolution of peptide ligands using GA
After calculation of the docking energies of the initial ten peptides (described in the Methods section), the peptides were evolved by GA, to produce the next generation. First, each of the ten peptides was ranked according to its calculated docking energy, and then the top seven peptides were remained for evolution by GA. The GA program then duplicated the top three sequences and formed a pair with the remaining (4 th -7 th ) peptides. Then ten peptides were recombined into the sequence. Finally, a 10-20% mutation factor was introduced to increase sequence variety in the next generation. Each process from docking to GA evolution was repeated for four rounds.

Characterization of peptide effects on PQQGDH activity
The three peptides with the lowest docking energies, GEKD, GERD, and SERG, and the control peptide DDDD were synthesized by Qiagen (Chatsworth, CA, USA). The PDB file of PQQGDH is a monomer and does not contain Ca 2+ . In general, purified PQQGDH is a mixture of an apo-and holo-type enzyme. Therefore, we dialyzed the enzyme against 100 mM EDTA and 10 mM MOPS-NaOH buffer for 2 days, to remove the Ca 2+ ions. Each synthesized peptide was prepared at concentrations of 0.2, 1, 2, 10, and 20 µM. The prepared peptides and 0.57 nM of apo-PQQGDH were incubated in 10 mM MOPS-NaOH buffer (pH 7.0) for 10 min at room temperature. PQQGDH activity was measured using 0.6 mM PMS and 0.06 mM DCIP. Decreased in the absorbance of DCIP at 600 nm at various concentrations of glucose were monitored, and the kinetic parameters of PQQGDH were as described previously [26]. The kinetic parameters and Ki values were calculated using the KaleidaGraph software.

Surface plasmon resonance (SPR) data analysis
The GEKD and SERG peptides, which showed inhibitory effects on GDH, were synthesized with a biotin modification (Qiagen).
The binding features of the peptide ligand to PQQGDH have been studied using SPR technology based on the BIAcoreX instrument (BIAcore). The biotinylated GEKD and SERG peptides were immobilized onto the Biacore SA sensor chip using the streptavidin-biotin interaction, and 120 µg/ml of each peptide in 100 mM NaCl, 10 mM phosphate buffer (pH 7.0) was applied a flow rate of 5 µl/min for 16 min. The resonance signals of peptides GEKD and SERG reached 273 and 300 resonance units (RUs). PQQGDH for use as an analyte was prepared at concentrations of 446 µg/ml, 223 µg/ml, 89 µg/ml, and 47 µg/ml in 10 mM MOPS-NaOH buffer (pH 7.0). Analytes were injected into the immobilized peptides at a flow rate of 20 µl/ml for 4 min at 25°C. All of the SPR experiments were performed in 10 mM MOPS-NaOH buffer (pH 7.0), and 0.5% SDS was chosen for surface regeneration. The sensor-gram signal of PQQGDH with each peptide was normalized by the nothing-immobilized sensor chip. The equilibrium dissociation constant (K D ) value was calculated from the Scatchard plot.