Hybrid method to solve HP model on 3D lattice and to probe protein stability upon amino acid mutations

Background Predicting protein structure from amino acid sequence is a prominent problem in computational biology. The long range interactions (or non-local interactions) are known as the main source of complexity for protein folding and dynamics and play the dominant role in the compact architecture. Some simple but exact model, such as HP model, captures the pain point for this difficult problem and has important implications to understand the mapping between protein sequence and structure. Results In this paper, we formulate the biological problem into optimization model to study the hydrophobic-hydrophilic model on 3D square lattice. This is a combinatorial optimization problem and known as NP-hard. Particle swarm optimization is utilized as the heuristic framework to solve the hard problem. To avoid premature in computation, we incorporated the Tabu search strategy. In addition, a pulling strategy was designed to accelerate the convergence of algorithm based on the characteristic of native protein structure. Together a novel hybrid method combining particle swarm optimization, Tabu strategy, and pulling strategy can fold the amino acid sequences on 3D square lattice efficiently. Promising results are reported in several examples by comparing with existing methods. This allows us to use this tool to study the protein stability upon amino acid mutation on 3D lattice. In particular, we evaluate the effect of single amino acid mutation and double amino acids mutation via 3D HP lattice model and some useful insights are derived. Conclusion We propose a novel hybrid method to combine several heuristic strategies to study HP model on 3D lattice. The results indicate that our hybrid method can predict protein structure more accurately and efficiently. Furthermore, it serves as a useful tools to probe the protein stability on 3D lattice and provides some biological insights.

One issue for 2D lattice model is that it's too simplified to constrain the amino acid sequence on a 2D plane. One step forward is to fold the sequence on 3D lattice and make it a better and native approximation. So far, several algorithms have been applied for 3D HP protein structure prediction problem, such as UEGO (Universal Evolutionary Global Optimization) [6], GA (Genetic Algorithms) [7], TS (Tabu Search) [8], EA (Evolutionary Algorithm) [9] and so on. Each method has its advantage to capture some special structure in the problem. In this paper, we aim to propose a hybrid method and improve the efficiency to solve the 3D HP protein structure prediction problem.
PSO was introduced by Kennedy and Eberhart [10]. It is a swarm intelligence optimization algorithm which imitates the foraging behaviors of birds and fish. As a simple meta-heuristic, it has been used to solve optimization problem with nonlinear, non-differentiable, and multi-modal function. Originally, this algorithm was designed for solving continuous optimization problem. Here, we started from the basic PSO framework and firstly extend the algorithm to the combinatorial optimization, into which we formally formulate the HP model on 3D lattice. In addition, we improved PSO as follows: a) redefined velocity for discrete model; b) employed modified Tabu search strategy to avoid premature convergence; c) designed pulling strategy to speed up convergence.
We showed that our hybrid algorithm can predict structures of amino acid sequences with different length efficiently. With this useful tool, we simulated the effects after single amino acid mutation and double amino acids mutation, respectively. Some biological insights are obtained.
The remainder of this paper is organized as follows. Firstly, a mathematical model was established for 3D HP problem. Secondly, we explained the PSO algorithm and proposed modified Tabu search method and pulling strategy. Thirdly, the performance of our algorithm was validated. Fourthly, the amino acid mutation result was obtained and analyzed. Finally, conclusions were presented.

Combinatorial optimization formulation for 3D HP lattice model
In HP model, every amino acid sequence is abstracted as an alphabetic string with H (hydrophobic amino acid) and P (hydrophilic amino acid). The protein conformation is a self-avoiding path on a 2D lattice. It is assumed that the main driving forces of the formation of the tertiary structure are the interactions among hydrophobic amino acids which are adjacent on lattice but not adjacent in the sequence, denoted as H-H interactions. The free energy of a protein conformation (X) is expressed by the number of H-H interactions. Based on Anfinsen's assumption [11], the configuration tends to form a core in the spatial structure shield from the surrounding solvent by hydrophilic amino acids with the minimal free energy. So the more H-H interactions, the lower the free energy. We assumed that the free energy equals to the minus number of H-H interactions. HP lattice model has been used for solving protein structure prediction problem on 2D and 3D lattices widely. In this paper, we focused on the 3D HP square lattice model.
Based on the abstraction and minimum energy principle, we established the optimization model (OM) for protein structure prediction problem on 3D square lattice as following: Fig. 1 The structure of sequence (H 2 P 2 H 5 P 2 H) on 3D square lattice min E(X) (1) Here, x i,j,k (l) = 1 if the X(l) = (i, j, k) 0 else (8) f (l) = 1 if the l th amino acid is H 0 if the l th amino acid is P Where, E(X) is the free energy of protein conformation X, X(l) d is the d th component of X(l), M(X) is the number of H-H interactions in conformation X, r expresses the number of adjacent hydrophobic pairs in amino acid sequence and · is Hamming distance. Equations (2), (3) and (4) constrain that every amino acid occupies only one lattice point, each lattice point cannot be used more than once and adjacent amino acids in the chain occupy the adjacent points on the lattice. Equation (8) presents whether the l th amino acid occupies point (i, j, k). In Eq. (9), f (l) translates the l th H (or P) of the amino acid sequence into 1 (or 0).
Solving the simplified HP model is NP-complete even on two dimensional lattice. Then we have to seek help from heuristic algorithms. Particle swarm optimization, one of the stochastic algorithm, serves as a powerful approximation method.

Hybrid algorithm The basic PSO algorithm
Particle swarm optimization (PSO) is a heuristic framework that optimizes an objective function by iteratively improve a candidate solution. The motivation is to have a population of candidate particles, and move these particles around in the search-space according to simple mathematical formulae over the particle's position and velocity. Each particle's movement is influenced by its local best known position, but is also guided toward the best known positions in the search-space, which are updated as better positions are found by other particles. Finally it is expected to move the swarm toward the best solution. The advantage of PSO is that it makes no assumptions about the problem and can search very large spaces of candidate solutions.
In basic PSO algorithm (See Table 1), m particles search the optimal position simultaneously with dynamic velocity. Particle velocity is affected by iteration, own cognition, and social cognition of particle. Particularly, each particle can remember not only its own flight experience, but also the trajectories of all particles. In n dimensional search space, the position and velocity of the i th particle are represented as X i ∈ R n and V i ∈ R n , respectively. They are updated by the following two equations: Where P t ib and P t gb are the best position of the i th particle and the best position of all particles in the t th iteration, respectively. Inertia weight (ω), self confidence (c 1 ) and swarm confidence (c 2 ) are input parameters, r 1 , r 2 are two separately generated uniformly distributed random numbers in the range [0,1].

The modified PSO algorithm
Definitions To solve the optimization model, we redefined position and velocity of PSO on 3D lattice. Particle position was orderly expressed by protein conformation (X). Velocity of particle was defined as a series of shift (j 1 , j 2 ), which means that the j th 1 component of particle position becomes the j th 2 component, then the j th 1 component and the j th 2 component (including the j th 2 component) were changed subsequently. In addition, position X 1 was obtained by the sum of position X 2 and a series of shift, namely Table 1 The process of basic PSO algorithm Step 1 To initialize {X 0 i |i = 1, 2 · · · m} and {V 0 i |i = 1, 2 · · · m}; Step 2 To calculate E(X t i ) , find P t ib and P t gb ;

Step 3
To update X t i and V t i ;

Step 4
To output P gb . Clearly, X + V is a new position. Nevertheless, the new position may not satisfy the constraints in the OM model. An adjustment strategy is needed to ensure the new position was valid.
Modified Tabu search strategy Premature convergence is one of the major difficulty to solve OM model by PSO algorithm. To further improve the modified PSO, we adopted the idea of Tabu search which was proposed by Glover [12]. This method was briefly described as follows.
Tabu search is a meta-heuristic method that maintains only one solution in the iteratively searching process. Given an initial solution X, the idea is to calculate and compare its neighboring solutions N(X). The best solution is chosen as candidate solution X c . If X c is satisfied with the aspiration rule, it will replace the current solution X and be added to tabu list T list ; Otherwise, the current solution X will be replaced by the best one X (E(X ) = min{E(X)|X ∈ N(X), X / ∈ T list }) and X will be added to T list . Generally, T list is a first-in first-out (fifo) memory with limited length. So particles would not search the solutions which have been found for a while, simultaneously, the better solutions would not always be taboo.
Neighbourhood of solution and aspiration rule are the key components of Tabu search. In our 3D HP problem, feasible solution is a 3D self-avoiding path. It was not easy to figure out its neighboring solutions from a given solution. According to Eqs. (10) and (11), we got similar solutions by changing r 1 , r 2 at the same iteration for the same particle in PSO, then these solutions constituted  Table 2 The algorithm outline of TPPSO 2
a neighbourhood. When candidate solution was better than the current solution, we would ignore whether the candidate solution was taboo or not.
Pulling strategy The convergence rate of modified PSO with Tabu search strategy is not fast enough and the conformations obtained by this modified PSO may be too loose. The following strategy was designed in order to improve the algorithm.
In native protein structure, hydrophobic amino acids concentrate inside of conformation and they were surrounded by hydrophilic amino acids. If hydrophilic amino acids were pulled out of the central of protein structure, the structure will be more compact and more stable. Without changing structure's legitimacy, this strategy was defined as pulling strategy. In order to make pulled structure to satisfy the self-avoiding constraints, only one amino acid could be pulled to its vacant diagonal position once. Figure 2 showed the move and result of one pulling.
The number in parentheses is the iteration number before the lowest free energy values are found. TPPSO 2 can find the optimal results of all sequences. TPPSO 1 can't obtain the minimal free energies for sequence A 6 , A 7 , and A 11 (highlighted in bold) Hybrid method A novel hybrid method was proposed by combining modified PSO with modified Tabu search strategy algorithm, denoted as TPPSO 1 . Another hybrid method was taken as TPPSO 2 , which combined TPPSO 1 with pulling strategy. Both methods employed Tabu search strategy and were applied to solve protein structure prediction problem. In TPPSO 1 and TPPSO 2 , when P ib and P gb were found, s alternative particles would be produced by Eqs. (10) and (11) for each particle.
We selected different r 1 and r 2 for finding alternative particles. These alternative particles might not satisfy the constraints, therefore they should be adjusted. Then the best alternative particle would replace the Fig. 3 This is one of structures for sequence A 6 . This optimal conformation was simulated by TPPSO 2 with 12 H-H interactions. Squares are for hydrophobic amino acids, and circles are for hydrophilic amino acids. In this structure all hydrophobic amino acids are surrounded in center. It is stable with minimal free energy Squares are for hydrophobic amino acids, and circles are for hydrophilic amino acids. In this structure almost all hydrophobic amino acids are surrounded in center. It is stable with minimal free energy previous particle and P gb would be taboo in a period of time. Differently, pulling strategy has been used in TPPSO 2 , so each particle could be closer to optimal position. Table 2 showed the detailed procedures of TPPSO 2 .

Numerical simulations
In order to test the feasibility of the hybrid algorithms (TPPSO 1 and TPPSO 2 ) and explore the properties of algorithms, we calculated two groups of amino acids sequences, respectively.  The average CPU time of our methods and hELP. The abscissa is the number of sequence, and the ordinate is CPU time. Because TPPSO 1 can't obtain the minimal free energy of sequence A 6 , A 7 and A 11 , we chose smaller CUP time with minimal free energy for all sequences, denoted as TPPSO. In the figure, the CPU time of TPPSO 1 and TPPSO 2 are stable with respective optimal structure. CPU time of TPPSO also is stabler. But CPU time of hELP is fluctuant

Simulation of sequences with 27 amino acids
We selected 11 sequences with 27 amino acids (See Table 3) which were also computed by EN [13] and hELP [14]. These sequences were used to test the performances of TPPSO 1 (without pulling strategy) and TPPSO 2 (with pulling strategy), respectively. In TPPSO 1 and TPPSO 2 , the inertia weight ω was updated by the following formula: The Time is the circular times and Maxtime is the maximum number of iterations which is 3000 in our implementation. For each particle, we chose c 1 = c 2 = 1, r 11 = rand(0.9, 1), r 12 = rand(0.82, 0.92), r 13 = rand(0.74, 0.84), r 21 = rand(0.9, 1), r 22 = rand(0.85, 0.95), r 23 = rand(0.8, 0.9) to produce three similar but not identical alternative particles. In this test, T list only contained ten particles.
According to Table 4, we knows that all the sequences in Table 3 were simulated by EN, hELP, and our method TPPSO 1 and TPPSO 2 . hELP and TPPSO 2 can obtain the minimal free energy of every sequence, but EN and TPPSO 1 can't find the minimal free energies of sequence A 6 , A 7 , and A 11 which are bigger. It illustrated our method can successfully predict the protein structure on 3D square lattice. The number in parentheses is the iteration number when the lowest free energy values are found. By comparing the results of hELP with TPPSO 1 and TPPSO 2 , TPPSO 1 is superior to hELP, and TPPSO 2 can fold stable structures earlier than TPPSO 1 .
Especially, TPPSO 2 found the lowest free energies of sequences A 6 , A 7 , and A 11 , while TPPSO 1 or EN did not. So we enumerated the structures of these sequences (See Figs. 3, 4, and 5), which were computed by TPPSO 2 . It can be seen that these conformations were more native, furthermore, the hydrophobic amino acids were concentrated and surrounded by hydrophilic amino acids. Because pulling strategy of TPPSO 2 has not been employed in TPPSO 1 , it is understandable that pulling strategy was able to accelerate the convergence of algorithm and optimize the protein structure.
The average CPU time of hELP was summarized in reference [14]. We also computed the average CPU time of TPPSO 1 and TPPSO 2 . The average CPU time of all methods were shown in Fig. 6. It is obvious that the average CPU time of every sequence of TPPSO 1 is the shortest, and that of TPPSO 2 is longer, because TPPSO 2 added the pulling strategy. For every sequence, the average CPU time  (13) where x is the number of H-H pairs, and y is the IBE number. The figure of fitting function was exhibited in Fig. 7. Except for sequences A 4 and A 10 , the IBE number of others are all close to the fitting function. The IBE number of sequence A 4 and A 10 are not satisfied with the fitting function, because in these sequences H amino acids and P amino acids are very dispersive, but in other sequences H segments or P segments are longer. We believed that IBE number of TPPSO 2 is mainly affected by the number of H-H interactions for sequences with the same length. It tends to be larger with more H-H interactions. Moreover, the length of H or P segments will affect the IBE number.
In order to further verify the above conclusion, we simulated three test sequences with the same length (See Table 6). It is obvious that H segments and P segments of these test sequences are longer. IBE numbers for test sequences are close to the fitting curve (See Fig. 7) and all relative errors were showed in Table 6. It means that our inference about IBE number of TPPSO 2 is reasonable.

Simulation of sequences with different length
We also computed several sequences with different length which have not been solved by EN and hELP. Moreover, Table 7 recorded the H-H interactions and IBE number of TPPSO 2 . These sequences were simulated by TPPSO 1 and TPPSO 2 respectively. The results of two methods are the same (See Table 7). They have the same H-H interactions. But we knows that the CPU time of TPPSO 2 is shorter than one of TPPSO 1 , because TPPSO 2 includes pulling strategy. For this reason, the structure obtained by TPPSO 2 is more compact. It is illustrated in Fig. 8. These results shows that: a) TPPSO 2 is able to solve sequences with different length and the obtained characteristic of protein structure is significant. b) pulling strategy improved the performance. c) Tabu search strategy  avoided prematurity effectively. d) For TPPSO 2 , the longer the sequence, the more the IBE number.

Probing protein stability upon amino acid mutation
Protein stability determines whether a protein will be in its native folded conformation or a denatured state. The folded, biologically active conformation of a protein is believed more stable than the unfolded, inactive conformations [15]. Thus, making proteins more stable is important in medicine and basic research. Amino acid mutations are widely used in protein design and analysis techniques to increase or decrease stability. These mutations are carried out experimentally using sitedirected mutagenesis and similar techniques. This is timeconsuming and often requires the use of computational prediction methods to select the best possible combinations [16][17][18][19]. With the efficient hybrid method at hand, we aim to probe the protein stability on 3D lattice. Particularly, we will simulate how single-site or double amino acid mutation affects protein stability. i.e., predicting the protein stability changes upon amino acid mutations with TPPSO 2 .

Single amino acid mutation
The hybrid method TPPSO 2 has been tested to solve protein structure prediction problem. Now, we focused on single amino acid mutation, whether and which amino acid affects the stability of protein structure. The experiments is designed as follows.   Sequences with different length In order to probe the stability of amino acid mutation, we chose four sequences with different lengths. These sequences were mentioned in the above section. They are sequence B, C, D and A 8 . Figure 9 recorded the H-H interactions of every single amino acid mutational sequence. From the results, we found that some mutational amino acids will result in a bigger deviation. We call those pivotal amino acids. The ratios of pivotal amino acids are 100%, 100%, 71.4% and 50% respectively for the four sequences. Also we deduced that those pivotal amino acids tend to locate at the beginning or end of sequence.
Tables 8, 9, 10 and 11 recorded the characters of mutated sequences including the deviation between mutated sequence and original sequence, the quantity of amino acid which mutated to cause the deviation and the ratio of every deviation. According to the results, we found that single amino acid mutation has the maximal and minimal deviation 2 and -2. The results also indicated that the ratio of maximal deviation is around 22%. Table 12 summarized the effect of hydrophobic (hydrophilic) amino acid mutation on H-H interactions. We know that hydrophobic amino acid mutation would not make the H-H interactions increase and hydrophilic amino acid mutation would not lead the H-H interactions decrease. It means that hydrophilic amino acid mutation will result in more compact structure, while hydrophobic amino acid mutation will result in the looser structure. By comparing the results of H 2 and P 2 in the Table 12, we suppose that hydrophobic amino acid is more impressible than hydrophilic amino acid to reflect stability of protein structure for sequence with different lengths.
The structures of sequence B before and after mutation are showed in Fig. 10. Since the 14th amino acid was changed from P to H, the number of H-H interaction increases and the deviation is 2, which is the maximal  H 0 and P 0 are the number of hydrophobic and hydrophilic in the original sequence. H 1 is the number of mutational hydrophobic amino acid whose H-H interactions is not more than the original one. P 1 is the number of mutational hydrophilic amino acid whose H-H interactions is not less than the original one. H 2 is the number of hydrophobic amino acid which caused the maximal deviation with original H-H pairs. P 2 is the number of hydrophilic amino acid which caused the maximal deviation with original H-H pairs deviation. It is obvious that optimal structures of 14th amino acid mutation is more compact. The structures of sequence D before and after mutation are showed in Fig. 11. Since the 6th amino acid was changed from H to P, the number of H-H interaction decreases and the deviation is 2 which is the maximal deviation. It is obvious that optimal structure of original sequence is more compact.

Sequences with the same length
We selected five sequences from Table 3 to test what kind of protein structures are more stable upon single amino acid mutation by TPPSO 2 . We changed every amino acid of these sequences, then recalculated and recorded the H-H interactions of every mutated sequence. Table 13 showed that: 1) For the sequences with the same length, UC number is likely larger with more mass number. The exception might be caused by the fact that UC number reflected the stability of protein structure. Namely, sequence with larger UC number was susceptible to single amino acid mutation. 2) According to H-H interactions and the number of hydrophobic amino acids, the sequence with more hydrophobic amino acid usually had more H-H interactions. 3) From P→H and H→P, most of mutational hydrophobic amino acids made the H-H interactions changed. Relatively, only a part of hydrophilic amino acids affect the number of H-H interactions. We conclude that hydrophobic amino acid is more impressible than hydrophilic amino acid to reflect stability of protein structure for sequence with the same length.
All these results illustrated that: a) the more hydrophobic amino acids, the more H-H interactions; b) sequence with more H-H interactions tends to be more stable when single amino acid is mutated; c) hydrophobic amino acid mutation tends to alter the protein structure largely.
According to the above observations, we summarize that the sequence with more hydrophobic amino acids will be less susceptible to single amino acid mutation.

Double neighbouring amino acids mutation
Amino acid does not work alone and multiple amino acids coordinate to maintain stability and perform function. Our in-silicon simulation allows us to go beyond single amino acid mutation and explore the combinatorial effect of amino acid mutation. In this section, we explore the effect of double neighbouring amino acids mutation (two adjacent amino acids are mutated) in protein folding. Double neighbouring amino acids mutations were classified as HH → PP, PP → HH, HP → PH, PH → HP.
We simulated three sequences (B, C, D) with different length when adjacent amino acids mutated. The maximal deviation and locations of pivotal amino acids were conserved (See Tables 14, 15, 16).

Fig. 10
These figures are the structures of sequence B before mutating and after mutating respectively. a is the optimal structure with 9 H-H interactions predicted for original sequence B by TPPSO 2 . b is one of optimal structures predicted for mutated sequence B by TPPSO 2 , in which the 14th amino acid was mutated. The mutated amino acid is denoted by arrow in figures. Since the 14th amino acid was changed from P to H, the number of H-H interaction increases and the deviation is 2. It is obvious that right structure is more compact Fig. 11 These figures are the structures of sequence D before mutating and after mutating respectively. a is the optimal structure with 9 H-H interactions predicted for original sequence B by TPPSO 2 . b is one of optimal structures predicted for mutated sequence C by TPPSO 2 , in which the 6th amino acid was mutated. The mutated amino acid is denoted by arrow in figures. Since the 6th amino acid was changed from H to P, the number of H-H interaction decreases and the deviation is 2. It is obvious that left structure is more compact Tables 14, 15 and 16 recorded the variation of H-H interactions and the position of pivotal double amino acids . According to these tables, we concluded that: a) If double amino acids mutation was HH → PP or PP → HH, the H-H interactions must be changed. But PH→HP and HP→PH maybe have variation. b)HH → PP and PP → HH must make the H-H interactions decrease and increase, respectively. c) The effect of double adjacent amino acids mutation which belongs to HP → PH or PH → HP was finite. d) The position of pivotal double adjacent amino acids mutation tend to locate be at the head or tail of sequence.

Double arbitrary amino acids mutation
We continued to explore the combinatorial effect of amino acid mutation. In this section, we check the effect of double amino acids mutations with arbitrary distance in protein folding. The amino acid mutations were classified ed as HH → PP, PP → HH, HP  Table 7 with 20 amino acids. There are 10 hydrophobic amino acids and 10 hydrophilic amino acids in sequence B. We folded the conformations of all of mutation sequences.
Form Table 17, we knew that combinations of H+H and P+P are more sensitive and easier to affect the stability of protein structure. We simulated every mutational sequence. The results showed that a) all HH→PP mutations will decrease HH interactions, namely HH interactions won't increase; b) 43 PP→HH mutations will increase HH interactions, 2 PP→HH mutations won't change HH interactions, in other words, HH interactions won't decrease; c) HP→PH or PH→HP mutations will not influence HH interactions.
The simulation results in Table 18 indicated that a) closer H and H (or P and P) can result in D-value; b) amino acid H 18 respectively matches amino acid H 20 and amino acid P 10 to obtain maximal deviation of structure, so amino acid H 18 is the most sensitive amino acid, which   According to the above observations, we summarized that a) double arbitrary amino acids mutation will be more sensitive to affect protein stability; b) double amino acids mutation with the same hydrophilic or hydrophobic property is more unstable than double amino acids mutation with different property; c) most of sensitive combinations are at the head or tail of sequence.

Discussion
As many research results indicate, HP model is very useful for modelling protein properties though it is simple and has many disadvantages. It captures the main difficulty of the real world problem. HP model has been applied in investigation of ligand binding to proteins [20]. The distinct influences of function, folding, and structure on the evolution of HP model are studied, by exhaustive enumeration of conformation and sequence space on a two dimensional lattice, which costs four week's computation [21]. These research all show that our effort to fold the HP chain by a hybrid method on 3D lattice is necessary and important.
Also we propose to use HP model to probe the protein stability. HP model serves as a very efficient tool here. The  H+H(P+P) means that arbitrary double hydrophobic(hydrophilic) amino acids will be mutated. H+P(P+H) means that hydrophobic(hydrophilic) will match with hydrophilic(hydrophobic) behind of it to mutate. O-num is original combination number in sequence. V-num is the number of combinations with which minimal free energy were altered after mutating. V-rate is the rate of V-num. The arrows in parentheses indicate increase or decrease of the free energy simplification of 20 amino acids to H, and P types dramatically reduce the possible mutation pattern. Especially we can easily perform the double mutation only considering four combinations. Those insights from the HP model can serve as novel hypothesis to guide experiments. We also need to point out that the protein stability results and conclusions are heavily depending on the optimal solution of 3D HP model. We demonstrate the results in some small scale problems. When we want to generalize the study, we need to further improve the hybrid algorithm. In our study, the computational experiments show that the new hybrid algorithm is efficient for short sequences. When the input space is bigger, there will be some suboptimal solutions and more difficult to find the minimal energy configurations. It's really a challenge for large scale HP model. The conformation space grows rapidly as the chain length increases. A possible method is to introduce divide-and-conquer strategy. We can also consider to combine with other algorithms or start from a good initial point from biological view. It will be our future work in devising such an algorithm for large protein.

Conclusion
In this paper, we studied protein structure prediction problem on 3D square lattice. We summarize the findings of this work as follows. Firstly, we formulated the protein structure prediction problem on 3D lattice into a combinatorial optimization problem; secondly, basic PSO algorithm has been enhanced to deal with discrete optimization problem; thirdly, we proposed a novel hybrid method (TPPSO 2 ) and proved its feasibility by simulating; fourthly, we derived some interesting insights for protein stability via single and double amino acid mutation perturbation.