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.

Background

Protein is the substantial basis of biological activity. The function of protein is determined by its structure which is believed to be decided by the amino acid sequence according to Anfinsen’s experiments. So the research on protein structure prediction (also called protein folding problem) is very significant and fundamental in exploring the fundamental principle to map sequence, structure, and function.

To capture the backbone of protein structure prediction, Dill and his collaborators introduced HP lattice model to simplify real world complexity in 1995 [1]. HP lattice model is an abstracted scaffold, and eventually convert the protein structure prediction problem to an optimization problem on lattice. The aim is to find the optimal structure with the lowest energy. Computationally, solving this problem is NP-hard. For this reason many researchers have been attracted to study this problem by proposing many heuristic algorithms. In recent years, for 2D HP protein folding problem, many methods have been proposed, e.g., PSO (Particle Swarm Optimization) [2], ACO (Ant Colony Algorithm) [3], ABO (Artificial Bee Colony) [4] and SOM (Self-Organizing Mapping) [5] etc.

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.

Methods

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.

At present, relative coordinates and space coordinates have been used to denote the protein conformation. For a sequence S with L amino acids, X is a string of length L−1 over the symbols {r(ight),l(eft),f(orward),d(own),u(p)} in relative coordinates, these five symbols reflect the relative location of contiguous amino acids on lattice. In space coordinates, X records the 3D coordinates of L amino acids, namely, X=(X(1),X(2)⋯X(L)) and X(l)∈N 3(l=1,2⋯L) is the coordinate of the l amino acid. In this paper, we chose the space coordinates. For example, Fig. 1 showed a conformation with 7 H-H interactions on 3D square lattice. Its conformation was denoted as X=((2,3,2),(3,3,2),(3,4,2),(3,4,3),(3,3,3),(2,3,3),(2,2,3),(3,2,3),(3,2,2),(3,1,2),(2,1,2),(2,2,2)).

Fig. 1

The structure of sequence (H 2 P 2 H 5 P 2 H) on 3D square lattice

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:

Here,

Where, E(X) is the free energy of protein conformation X, X(l) is the d 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 amino acid occupies point (i,j,k). In Eq. (9), f(l) translates the l 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 particle are represented as X ∈R and V ∈R , respectively. They are updated by the following two equations:

Table 1

The process of basic PSO algorithm

Step 1	To initialize {Xi0|i=1,2⋯m} and {Vi0|i=1,2⋯m};
Step 2	To calculate E(Xit), find Pibt and Pgbt ;
Step 3	To update Xit and Vit;
Step 4	To output P gb.

Where and are the best position of the i particle and the best position of all particles in the t 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 component of particle position becomes the component, then the component and the component (including the component) were changed subsequently. In addition, position X 1 was obtained by the sum of position X 2 and a series of shift, namely X 1=X 2+{(j ,j )}. For example, V={(2,4),(3,1)} and X=(X(1),X(2),X(3),X(4)), then

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 . If X is satisfied with the aspiration rule, it will replace the current solution X and be added to tabu list T ; Otherwise, the current solution X will be replaced by the best one X ′ (E(X ′)=min{E(X)|X∈N(X),X∉T }) and X ′ will be added to T . Generally, T 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 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.

Fig. 2

The flow of pulling strategy. Left figure is for structure before pulled. Right figure is for structure after pulled

Hybrid method

A novel hybrid method was proposed by combining modified PSO with modified Tabu search strategy algorithm, denoted as TPPSO1. Another hybrid method was taken as TPPSO2, which combined TPPSO1 with pulling strategy. Both methods employed Tabu search strategy and were applied to solve protein structure prediction problem. In TPPSO1 and TPPSO2, when P and P 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 previous particle and P would be taboo in a period of time. Differently, pulling strategy has been used in TPPSO2, so each particle could be closer to optimal position. Table 2 showed the detailed procedures of TPPSO2.

Table 2

The algorithm outline of TPPSO 2

Step 1	To initialize {Xi0|i=1,2⋯m},{Vi0|i=1,2⋯m} and T list=Ø;
Step 2	To calculate E(Xit), find Pibt and Pgbt;
Step 3	To update {Vijt|j=1,2⋯s} and {Xijt|j=1,2⋯s};
Step 4	To adjust and pull {Xijt|j=1,2⋯s} ;
Step 5	To calculate E(Xict)=min{E(Xijt)|j=1,2⋯s};
Step 6	If E(Xict)≤E(Xit) then Xit=Xict;
Step 7	To calculate E(Pgbct)=min{E(Xit)|i=1,2⋯m};
Step 8	If E(Pgbct)<E(Xgbt) then Xit=Xict,T_list=Ø;
Step 9	If E(Pgbct)=E(Xgbt) and Pgbct∉T_list then T_list=T_list+Xgbt,Xgbt=Xgbct;
Step 10	To output P gb.

Results

Numerical simulations

In order to test the feasibility of the hybrid algorithms (TPPSO1 and TPPSO2) and explore the properties of algorithms, we calculated two groups of amino acids sequences, respectively.

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 TPPSO1 (without pulling strategy) and TPPSO2 (with pulling strategy), respectively. In TPPSO1 and TPPSO2, the inertia weight ω was updated by the following formula:

Table 3

Sequences with 27 amino acids used in our study

Sequence ID	Amino acids sequence
A 1	PHPHPH 3 P 2 HPHP 11 H 2 P
A 2	PH 2 P 10 H 2 P 2 H 2 P 2 HP 2 HPH
A 3	H 4 P 5 HP 4 H 3 P 9 H
A 4	H 3 P 2 H 4 P 3 HPHP 2 H 2 P 2 HP 3 H 2
A 5	H 4 P 4 HPH 2 P 3 H 2 P 10
A 6	HP 6 HPH 3 P 2 H 2 P 3 HP 4 HPH
A 7	HP 2 HPH 2 P 3 HP 5 HPH 2 PHPHPH 2
A 8	HP 11 HPHP 8 HPH 2
A 9	P 7 H 3 P 3 HPH 2 P 3 HP 2 HP 3
A 10	P 5 H 2 PHPHPHPHP 2 H 2 PH 2 PHP 3
A 11	HP 4 H 4 P 2 HPHPH 3 PHP 2 H 2 P 2 H

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 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 TPPSO1 and TPPSO2. hELP and TPPSO2 can obtain the minimal free energy of every sequence, but EN and TPPSO1 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 TPPSO1 and TPPSO2, TPPSO1 is superior to hELP, and TPPSO2 can fold stable structures earlier than TPPSO1.

Table 4

Comparing four algorithms in eleven sequences with 27 amino acids

Sequence ID	EN	hELP	TPPSO1	TPPSO2
A 1	-9	-9(18009)	-9(1983)	-9(177)
A 2	-10	-10(9447)	-10(1304)	-10(439)
A 3	-8	-8(1420)	-8(1249)	-8(44)
A 4	-15	-15(2125)	-15(795)	-15(19)
A 5	-8	-8(2877)	-8(104)	-8(61)
A 6	-11	-12(2610)	-11(940)	-12(812)
A 7	-13	-13(3967)	-12(721)	-13(805)
A 8	-4	-4(1070)	-4(6)	-4(3)
A 9	-7	-7(363)	-7(389)	-7(14)
A 10	-11	-11(416)	-11(2784)	-11(83)
A 11	-14	-16(285)	-14(957)	-16(2672)

The number in parentheses is the iteration number before the lowest free energy values are found. TPPSO2 can find the optimal results of all sequences. TPPSO1 can’t obtain the minimal free energies for sequence A 6, A 7, and A 11 (highlighted in bold)

Especially, TPPSO2 found the lowest free energies of sequences A 6, A 7, and A 11, while TPPSO1 or EN did not. So we enumerated the structures of these sequences (See Figs. 3, 4, and 5), which were computed by TPPSO2. 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 TPPSO2 has not been employed in TPPSO1, it is understandable that pulling strategy was able to accelerate the convergence of algorithm and optimize the protein structure.

Fig. 3

This is one of structures for sequence A 6. This optimal conformation was simulated by TPPSO2 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

Fig. 4

This is one of structures for sequence A 7. This optimal conformation was simulated by TPPSO2 with 13 H-H interactions. 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

Fig. 5

This is one of structures for sequence A 11. This optimal conformation was simulated by TPPSO2 with 16 H-H interactions. Squares are for hydrophobic amino acids, and circles are for hydrophilic amino acids. It is stable and compact with minimal free energy

The average CPU time of hELP was summarized in reference [14]. We also computed the average CPU time of TPPSO1 and TPPSO2. The average CPU time of all methods were shown in Fig. 6. It is obvious that the average CPU time of every sequence of TPPSO1 is the shortest, and that of TPPSO2 is longer, because TPPSO2 added the pulling strategy. For every sequence, the average CPU time of TPPSO1 and TPPSO2 are stable and vary around 0.4 and 0.8 s respectively. However, the average CPU time of hELP is not stable. Since TPPSO1 can’t obtain the lowest free energy of sequence A 6, A 7, and A 11, we made TPPSO as the method which can fold the optimal structures of all sequences by PSO. The average CPU time of TPPSO was taken as less average CPU time of TPPSO1 and TPPSO2,which was also showed in Fig. 6. We know that the average CPU times of all sequences of TPPSO and hELP are 0.475 and 1.41 s respectively. It indicated that our method TPPSO is faster.

Fig. 6

The average CPU time of our methods and hELP. The abscissa is the number of sequence, and the ordinate is CPU time. Because TPPSO1 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 TPPSO1 and TPPSO2 are stable with respective optimal structure. CPU time of TPPSO also is stabler. But CPU time of hELP is fluctuant

Table 5 summarized the number of H-H interactions and iteration number before the lowest free energy values are found by TPPSO2 (denoted as IBE number) for eleven sequences with 27 amino acids in the Table 3. It is obvious that the more H-H interactions, the more IBE numbers. The fitting function of the number of H-H interactions and IBE number was given as follows.

Table 5

IBE number of TPPSO 2

Sequence ID	A 8	A 9	A 3	A 5	A 1	A 2	A 10	A 6	A 7	A 4	A 11
H-H a	4	7	8	8	9	10	11	12	13	15	16
IBE number b	3	14	44	61	177	439	83	812	805	19	2672

aH-H means the number of hydrophobic-hydrophobic amino acid interactions for optimal structure

biteration numbers before the lowest free energy values are found

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 TPPSO2 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.

Fig. 7

The fitting figure of IBE number and H-H interaction for sequences in table 1 by TPPSO2. The abscissa is H-H interaction of every sequence, and the ordinate is IBM number. Almost all sequences are satisfied with this fitting figure. IBE number will increase with H-H interaction for sequences with the same length. The three stars is the IBE number of test sequence with the same length 27. Their fitting function values are almost matched to computed IBE numbers

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 TPPSO2 is reasonable.

Table 6

Test sequences

Sequence ID	Amino acids sequence	H-H	IBE number	Relative error
Test 1	H 4 P 5 HP 5 H 3 P 8 H	8	51 (52.3829)	0.0271
Test 2	H 4 P 5 HP 5 H 3 P 4 HP 3 H	9	167 (177.8417)	0.0649
Test 3	(HP 2 HP)5 HP	14	956 (921.1219)	0.0365

IBE number is the iteration number of every sequence by TPPSO2 before the minimal free energy was found. The number in parentheses is IBE number calculated by fitting function

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 TPPSO2. These sequences were simulated by TPPSO1 and TPPSO2 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 TPPSO2 is shorter than one of TPPSO1, because TPPSO2 includes pulling strategy. For this reason, the structure obtained by TPPSO2 is more compact. It is illustrated in Fig. 8.

Table 7

Sequences with different lengths

Sequence ID	Amino acids sequence	Length	H-H	IBE number
B	H 4 P 2 H 7 P 3 H	17	9	2
C	HPHP 2 H 2 PHP 2 HPH 2 P 2 HPH	20	11	11
D	P 2 HP 2 H 2 P 4 H 2 P 4 H 2 P 4 H 2	25	9	139
E	P 3 H 2 P 2 H 2 P 5 H 7 P 2 H 2 P 4 H(HP 2)2	36	17	432
F	P 2 H(P 2 H 2)2 P 5 H 10 P 6(H 2 P 2)2 HP 2 H 5	48	29	976

H-H interactions were the same by TPPSO1 and TPPSO2. IBE numbers were computed by TPPSO2. These sequences were simulated by different methods, and every method only folds a part of sequences. TPPSO2 found the minimal free energy of all sequences

Fig. 8

These figures are the structures of sequence B with 20 amino acids. a is one of structures by TPPSO1 with 11 H-H interactions. b is one of structures by TPPSO2 with 11 H-H interactions. By comparing, left structure is more compact, right structure is looser

These results shows that: a) TPPSO2 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 TPPSO2, 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 site-directed mutagenesis and similar techniques. This is time-consuming and often requires the use of computational prediction methods to select the best possible combinations [16-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 TPPSO2.

Single amino acid mutation

The hybrid method TPPSO2 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. We firstly calculate the optimal H-H interactions of original sequence by TPPSO2. Then we choose one amino acid to mutate, i.e., we change it from H (P) into P (H). Then we calculate the optimal H-H interactions of mutated sequence by TPPSO2. Finally the deviation of H-H interactions between mutated sequence and original sequence was recorded.

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.

Fig. 9

These figures are H-H interactions of sequences after single amino acid mutation. The abscissa is the location of mutated amino acid, and the ordinate is the number of H-H interaction for mutated sequence. The horizonal line is the number of H-H interaction of original sequence. Two vertical lines spit sequence into three equal parts. a is mutational results of sequence B. In this figure, 100% pivotal amino acids locate at beginning or ending of sequence B. b is mutational results of sequence C. In this figure, 100% pivotal amino acids locate at beginning or ending of sequence C. c is mutational results of sequence D. In this figure, 71.4% pivotal amino acids locate at beginning or ending of sequence D. d is mutational results of sequence A 8. In this figure, 50% pivotal amino acids locate at beginning or ending of sequence A 8

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 8

The single amino acid mutation results for sequence B

D-value	-3	-2	-1	0	1	2	3
Q-value	0	1	5	7	2	2	0
R-value	0%	6%	2%	41%	12%	12%	0%

There are 9 H-H interactions by TPPSO2 for original sequence B. Every amino acid would be mutational, namely H (P) was changed into P (H). D-value is the deviation of H-H interactions between new sequence and original sequence when single amino acid was mutated. Q-value is the number of amino acids caused the deviation. R-value is the ratio of amino acids. The ratio of amino acids which caused the maximal deviation is 18% (summarization of the numbers highlighted in bold)

Table 9

The single amino acid mutation results for sequence C

D-value	-3	-2	-1	0	1	2	3
Q-value	0	5	5	4	6	0	0
R-value	0%	25%	25%	20%	30%	0%	0%

There are 11 H-H interactions by TPPSO2 for original sequence C. Every amino acid would be mutational, namely H (P) was changed into P (H). The ratio of amino acids which caused the maximal deviation is 25% (summarization of the numbers highlighted in bold)

Table 10

The single amino acid mutation results for sequence D

D-value	-3	-2	-1	0	1	2	3
Q-value	0	7	1	6	11	0	0
R-value	0%	28%	4%	24%	44%	0%	0%

There are 9 H-H interactions by TPPSO2 for original sequence D. Every amino acid would be mutational, namely H (P) was changed into P (H). The ratio of amino acids which caused the maximal deviation is 28% (summarization of the numbers highlighted in bold)

Table 11

The single amino acid mutation results for sequence A 8

D-value	-3	-2	-1	0	1	2	3
Q-value	0	5	3	7	11	1	0
R-value	0%	19%	11%	26%	41%	3%	0%

There are 8 H-H interactions by TPPSO2 for original sequence A 8. Every amino acid would be mutational, namely H (P) was changed into P (H). The ratio of amino acids which caused the maximal deviation is 22% (summarization of the numbers highlighted in bold)

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 H2 and P2 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.

Table 12

Summary of the single mutation results

Sequence ID	H0	P0	H1	P1	H2	P2
B	12	5	12	5	1	2
C	10	10	10	10	5	0
D	9	9	16	16	7	0
A 8	6	21	6	21	5	1

H0 and P0 are the number of hydrophobic and hydrophilic in the original sequence. H1 is the number of mutational hydrophobic amino acid whose H-H interactions is not more than the original one. P1 is the number of mutational hydrophilic amino acid whose H-H interactions is not less than the original one. H2 is the number of hydrophobic amino acid which caused the maximal deviation with original H-H pairs. P2 is the number of hydrophilic amino acid which caused the maximal deviation with original H-H pairs

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 deviation. It is obvious that optimal structures of 14th amino acid mutation is more compact.

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 TPPSO2. b is one of optimal structures predicted for mutated sequence B by TPPSO2, 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

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.

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 TPPSO2. b is one of optimal structures predicted for mutated sequence C by TPPSO2, 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

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 TPPSO2. 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.

Table 13

Single amino acid mutation of sequences with 27 amino acids

Sequence ID	H-H	mass number	UC number	H	P	P →H	H →P
A 8	4	9	4	6	21	19	4
A 3	8	7	2	9	18	17	8
A 5	8	8	7	9	18	12	8
A 10	11	17	7	11	16	9	11
A 4	15	13	10	14	13	3	14

H-H is the number of H-H pairs of original sequence. The same continuous amino acids were taken as a mass. Mass number is for original sequence. UC number is the number of mutational amino acid which does not change the H-H interactions. H (P) is the number of hydrophobic (hydrophilic) amino acids in the original sequence. P →H (H →P) is the number of mutational hydrophilic (hydrophobic) amino acid which affected the H-H interactions

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).

Table 14

Double amino acids mutation results for sequence B

D-vale	HH → PP (9) a	PP → HH (3) a	HP → PH (2) a	PH → HP (2) a
-2	4(H,M,T) b	0	0	0
-1	5	0	1	0
0	0	0	0	0
+1	0	1	1	1
+2	0	2(T)	0	1(T)

aThe number in parentheses is the number of adjacent amino acids in sequences. D-value is the deviation of minimal free energy caused by neighboring mutational amino acids

bThe position of mutational double amino acids. H (M,T) means that the mutational double amino acids is at head (middle, tail) of the sequence

Table 15

Double amino acids mutation results for sequence C

D-vale	HH → PP (2)	PP → HH (3)	HP → PH (7)	PH → HP (7)
-3	1(H)	0	0	0
-2	1(T)	0	0	0
-1	0	0	3	3
0	0	0	4	4
+1	0	2	0	0
+2	0	1(H)	0	0

Table 16

Double amino acids mutation results for sequence D

D-vale	HH → PP (4)	PP → HH (11)	HP → PH (4)	PH → HP (5)
-2	4(H,M,T)	0	0	0
-1	0	0	0	3
0	0	0	2	1
+1	0	5	1	1
+2	0	6(M,T)	1(T)	0

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 → PH, PH → HP. We simulated the sequence B in 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.

Table 17

Double arbitrary amino acids mutation results for sequence B

Combination	O-num	V-num	V-rate
H+H	45	45(↓)	100%
P+P	45	43(↑)	96%
H+P	50	29(↑ ↓)	58%
P+H	50	18(↑ ↓)	38%

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

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 is at the tail of sequence; c) amino acid H 3 will cause D-value with hydrophilic (P), so amino acid H 3 is very sensitive to polar, it is at the head of sequence; d) matching amino acid H 7 with arbitrary amino acid P, HH interactions are invariable, so H 7 is obtuse, but by combining H 7 with arbitrary H, it is sensitive, this amino acid is in the middle of sequence; e) H 1 and H 20 are impressible for other H, but they are stable for arbitrary P, the mutations of H 1 and H 20 with all of P lead to decrease one more HH interaction, H 1 and H 20 are at the head and tail of sequence, respectively.

Table 18

Combination D-value and pivotal amino acids results for sequence B

Combination D-value and pivotal amino acids
H+H	-4	H 1 H 3,H 18 H 20
P+P	+2	P 4 P 5,P 4 P 13,P 5 P 8,P 8 P 17,P 10 P 13,P 11 P 16,P 13 P 16,P 16 P 19
H+P	-2	H 3 P 10,H 3 P 16,H 3 P 19,H 6 P 19
P+H	-3	P 10 H 18

H+H(P+P) means that arbitrary double hydrophobic(hydrophilic) amino acids will mutate. H+P(P+H) means that hydrophobic(hydrophilic) will match with hydrophilic(hydrophobic) behind of it to mutate. D-value is the maximal deviation of H-H interactions between new sequence and original sequence when double arbitrary amino acids mutation. H H means that the i animo acid matches the j amino acid to mutate, in which two mutational amino acids are hydrophobic

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 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 sub-optimal 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 (TPPSO2) and proved its feasibility by simulating; fourthly, we derived some interesting insights for protein stability via single and double amino acid mutation perturbation.