Volume 11 Supplement 4

## Selected papers from the 10th International Conference on Systems Biology (ISB 2016)

- Research
- Open Access
- Published:

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

*BMC Systems Biology***volume 11**, Article number: 93 (2017)

## Abstract

### 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*
^{th} 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)).

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*)_{
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_{1}^{th}\) component of particle position becomes the \(j_{2}^{th}\) component, then the \(j_{1}^{th}\) component and the \(j_{2}^{th}\) component (including the \(j_{2}^{th}\) 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*
_{
p
},*j*
_{
q
})}. 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*
_{
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 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.

##### 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 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}.

## Results

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

#### 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 of TPPSO^{1} and TPPSO^{2} are stable and vary around 0.4 and 0.8 s respectively. However, the average CPU time of hELP is not stable. Since TPPSO^{1} 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 TPPSO^{1} and TPPSO^{2},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.

Table 5 summarized the number of H-H interactions and iteration number before the lowest free energy values are found by TPPSO^{2} (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.

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 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 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. We firstly calculate the optimal H-H interactions of original sequence by TPPSO^{2}. 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 TPPSO^{2}. 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.

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

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.

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.

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

## References

- 1
Dill KA, Bromberg S, Yue K, et al. Principles of protein folding a perspective from simple exact models. Protein Sci. 1995; 4(4):561–602.

- 2
Guo YZ, Wu Z, Wang Y, et al. Extended particle swarm optimization method for folding protein on triangular lattice. IET Sys Bio. 2016; 10(1):30–33.

- 3
Nardelli M, Tedesco L, Bechini A. Cross-lattice Behavior of General ACO Folding for Proteins in the HP Model. Proc of the 28th Annual ACM Symp on Appl Comput. 2013; 18(22):1320–1327.

- 4
Zhang Y, Wu L. Artificial Bee Colony for Two Dimensional Protein Folding. Adv Electr Eng Syst. 2012; 1(1):19–23.

- 5
Zhang XS, Wang Y, Zhan ZW, Wu LY, Chen LN. Exploring protein’s optimal HP configurations by self-organizing mapping. J Bioinf Comput Biol. 2005; 3(02):385–400.

- 6
García-Martínez JM, Garzón EM, Cecilia JM, et al. An efficient approach for solving the HP protein folding problem based on UEGO. J Math Chem. 2015; 53(3):794–806.

- 7
Lin CJ, Su SC. Protein 3D HP model folding simulation using a hybrid of genetic algorithm and particle swarm optimization. Int J Fuzzy Syst. 2011; 13:140–147.

- 8
Benítez CMV, Lopes HS. Protein structure prediction with the 3D-HP side-chain model using a master-slave parallel genetic algorithm. J Braz Comput Soc. 2010; 16:69–78.

- 9
Tsay JJ, Su SC. An effective evolutionary algorithm for protein folding on 3D FCC HP model by lattice rotation and generalized move sets. Proteome Sci. 2013; 11(1):1.

- 10
Eberhart RC, Kennedy J. A new optimizer using particle swarm theory. Proc of the sixth Int Symp on micro Mach Hum Sci. 1995; 1:39–43.

- 11
Anfinsen CB. Principles that govern the folding of protein chains. Sci. 1973; 181(4096):223–230.

- 12
Glover F. Tabu search-part I. J Comput. 1989; 1(3):190–206.

- 13
Guo YZ, Feng EM. The simulation of the three-dimensional lattice hydrophobic-polar protein folding. J Chem Phys. 2006; 125(23):234703.

- 14
Liu J, Li G, Yu J, et al. Heuristic energy landscape paving for protein folding problem in the three-dimensional HP lattice model. Comput Biol Chem. 2012; 28:17–26.

- 15
Pascal L, Stefan G, Abdullah K, Valentina C, Paul JB, Christian M, Mering C, Paola P. Cell-wide analysis of protein thermal unfolding reveals determinants of thermostability. Science. 2017; 355(6327):812.

- 16
Parthiban V, Michael MG, Schomburg D. CUPSAT: prediction of protein stability upon point mutation. Nucleic Acids Res. 2006; 34(suppl 2):W239-W242.

- 17
Cheng J, Randall A, Baldi P. Prediction of Protein Stability Changes for Single Site Mutations Using Support Vector Machines, Proteins: Str. Func Bioi. 2006; 62:1125–32.

- 18
Shortle D, Stites WE, Meeker AK. Contributions of the large hydrophobic amino acids to the stability of staphylococcal nuclease. Biochem. 1990; 29:8033–41.

- 19
Perl D, Mueller U, Heinemann U, Schmid FX. Two exposed amino acid residues confer thermostability on a cold shock protein. Nat Struct Bio. 2000; 7(5):380–3.

- 20
Miller DW, Dill KA. Ligand binding to proteins: The binding landscape model. Prot Sci. 1997; 6(10):2166–79.

- 21
Blackburne BP, Hirst JD. Evolution of functional model proteins. J Chem Phys. 2001; 115(4):1935–42.

## Acknowledgements

The open access fee was supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (XDB13000000). This work is also supported by National Natural Fund under grant number 11601288, 11422108, 61621003, and 61304178.

### Availability of data and materials

Not applicable

### About this supplement

This article has been published as part of *BMC Systems Biology* Volume 11 Supplement 4, 2017: Selected papers from the 10th International Conference on Systems Biology (ISB 2016). The full contents of the supplement are available online at https://bmcsystbiol.biomedcentral.com/articles/supplements/volume-11-supplement-4.

## Author information

### Author notes

### Affiliations

### Contributions

YG developed and implemented the methods. YW and FT participated in the development of the methods. YW conceived the protein stability experiment. All authors draft, read, and approved the final manuscript.

### Corresponding author

Correspondence to Yuzhen Guo.

## Ethics declarations

### Ethics approval and consent to participate

Not applicable

### Consent for publication

Not applicable

### Competing interests

The authors declare that they have no competing interests.

## Rights and permissions

**Open Access** This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

## About this article

#### Published

#### DOI

### Keywords

- Protein structure prediction
- HP model
- 3D lattice
- Particle swarm optimization
- Protein stability