PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of jbiolphysspringer.comThis journalToc AlertsSubmit OnlineOpen Choice
 
J Biol Phys. 2009 August; 35(3): 245–253.
Published online 2009 June 11. doi:  10.1007/s10867-009-9152-1
PMCID: PMC2710459

Structure optimization of the two-dimensional off-lattice hydrophobic–hydrophilic model

Abstract

A two-dimensional off-lattice protein model with two species of monomers, hydrophobic and hydrophilic, was studied. Low-energy configurations in the model were optimized using the improved energy landscape paving (ELP+) method. In ELP+, the energy landscape paving (ELP) was first applied to search for the low-energy states. After the ELP led to the basins of the local energy minima, the additional degree-of-freedom of bond length was introduced, and the gradient descent method was then used to search for lower energy states near the local minima. Numerical results show that the proposed methods are quite effective for finding the ground states of proteins. A comparison between ELP+ and other methods is made.

Keywords: Protein structure prediction, Two-dimensional off-lattice model, Energy landscape paving, Local search

Introduction

Predicting the native structure of proteins from their amino acid sequences is one of the most challenging problems in bioinformatics. Despite many decades of research we still lack a detailed understanding of the relation between chemical composition and structure (and consequently function) of proteins. The theoretical analysis of the protein structure prediction problem, or the protein folding problem, faces two major difficulties. One is the determination of the potential energy function. The effective energy function can generally distinguish the native states from non-native states of protein molecules. The other is that detailed protein models are often characterized by a rough energy landscape with a huge number of local minima separated by high-energy barriers.

Since solving such a problem is too difficult for realistic protein models, in recent years the theoretical community has introduced and examined some highly simplified, but still nontrivial, models which make the following assumption: instead of considering all 20 different kinds of amino acids in real proteins, the models comprise only two prototypes of residues: hydrophobic (H) and the hydrophilic (or polar, P) monomers. These models include a large family of HP lattice models [13] and HP off-lattice models [4, 5]. Even though these models are highly simplified, to solve the corresponding folding problem remains NP-hard.

In this article, we introduced a so-called AB model by Stillinger et al. [6, 7], where the hydrophobic monomers are now labeled by A and the hydrophilic or polar ones by B. Based on the minimum free-energy theory [8], many authors have developed various techniques to search for the low-energy states of AB proteins in the two-dimensional (2D) case. The methods of optimizing low-energy states include conventional Metropolis-type Monte Carlo (MMC) procedures [7], the memory Tabu search (MTS) method [9], the off-lattice pruned–enriched Rosenbluth method (PERM) and, after subsequent conjugate gradient minimization (PERM+) [10], the annealing contour Monte Carlo (ACMC) method [11], as well as conformational space annealing (CSA) [12] and statistical temperature annealing (STA) methods [13]. For four amino acids chains with lengths 13 ≤ n  55 studied in Refs. [6, 7, 913], in this article we introduced another global optimization method, the energy landscape paving (ELP) method [14, 15], to optimize low-energy configurations. After the ELP led to the basins of the local energy minima, the additional degree-of-freedom of bond length was introduced, and the local search based on gradient descent was then used to search lower energy states near the local minima. Numerical results show that the proposed methods are very promising for finding the ground states of proteins.

Two-dimensional AB model

For an n-mer chain living in two dimensions, the distances between consecutive monomers along the chain are fixed to unit length, while nonconsecutive monomers interact through a modified Lennard–Jones potential. There is an angle θi between bond vectors ui and ui+1, where θi (1 ≤ i  n  2) denotes the bond angle at monomer i + 1 and satisfies –πθi < π. When the rotating direction of the angle is counterclockwise, its value is positive. The shape of the n-mer is either specified by the n  2 bond angles θ1,..., θn  2 or by the n  1 bond vectors u1, u2, ..., un−1. The definitions of the bond angle θi and the bond vector ui are shown in Fig. 1. The energy function [6, 7, 913, 16, 17] consists of two types of contributions: bond angle and Lennard–Jones. It can be written as

equation M1

Here, rij is the distance between monomers i and j (with i < j ). Each ζi is either A or B, and C(ζi, ζj) is +1, +1/2, and −1/2, respectively, for AA, BB, and AB (or BA) pairs, giving strong attraction between AA pairs, weak attraction between BB pairs, and weak repulsion between A and B.

Fig. 1
Definitions of bond vector ui and bond angle θi

Obviously, the AB model with no explicit representation of side chains or hydrogen bonds is able to provide only a coarse-grained approximation to the complexities of real proteins. However, its off-lattice constructions and the uses of both Lennard–Jones and bond-angle contributions still enable it to capture some of the basic features of protein structures.

The protein structure prediction problem for the AB model can be formally defined as follows: given a monomer chain s = ζ1ζ2ζ3...ζn, we try to find an energy-minimizing configuration of s in the 2D plane, i.e., find X* [set membership] G (s) so that E (X*) = min {E (X)|X [set membership] G(s)}, where G(s) is the set of all the valid configurations of s.

Theory and methods

Energy landscape paving method

In order to locate low-energy states of the sequences under consideration, we used the so-called energy landscape paving (ELP) minimizer [14, 15], which combines ideas from energy landscape deformation [18] and Tabu search [19]. The ELP minimization is a Monte Carlo (MC) optimization method, but with a modified energy expression designed to steer the search away from regions that have already been explored. This means that if a state X with energy E(X) is hit, the energy E is increased by a “penalty” and replaced by energy equation M2. In this expression, the “penalty” term f (H(q, t) is a function of the histogram H(q, t) in a prechosen “order parameter” q. The histogram is updated at each MC step, thus becoming a function of “time” t. The statistical weight for a state X is defined as

equation M3

where kBT is the thermal energy at the (low) temperature T, and kB is the Boltzmann constant. In the ELP simulations of the 2D AB model, we used the potential energy itself as an order parameter and chose equation M4 as the replacement of the energy E. The temperature was set to T = 5 K.

In a regular low-temperature simulation, the probability to escape a local minimum depends on the height of the surrounding energy barriers. Within ELP, the sampling weight of a local minimum state decreases with the time the system stays in that minimum, and consequently the probability to escape the minimum increases until the local minimum is no longer favored. The system will then explore higher energies until it falls into a new local minimum. In practice, the accumulated histogram function H(E, t) from all previously visited energies at the MC step t helps the simulation escape local entrapment and surpass high-energy barriers more easily.

Configuration update mechanism

In ELP, each MC step must update the current configuration. We used alternately the following gradient descent update mechanism and a rotary update mechanism to update configuration.

Gradient descent update mechanism

Consider the current configuration X = (X1, X2,....,Xn), where Xi = (xi, yi) denotes the position vector of the ith amino acid (i = 1, 2, ..., n). Since the length of the bonds is fixed ([mid ]ui[mid ] = 1, i = 1, 2, ..., n  1), the ith amino acid lies on a circle with radius of unit length around the (i  1)th amino acid. For updating the configuration X, we produced randomly a positive integer j(2 ≤ j  n, where n denotes chain length) and changed temporarily the position of the amino acid j to that of the amino acid j * in the anti-gradient direction of the energy function E(X) at the amino acid j, that is, Xj * = equation M5, where Xj and Xj * denote the position vectors of the amino acids j and j *, respectively, ε = 0.2 is the iterative step length factor, and equation M6 is the iterative search direction. Link up the two points j   1 and j *, and the linked line intersects the circle at a point equation M7. Move amino acid j to the position of equation M8 and amino acids j + 1, j + 2, ..., n are wholly translated as a rigid body with fixed trend of the partial chain (see Fig. 2). We thus obtained a new configuration N. Obviously, all the bond vectors in the configuration N are still unit length.

Fig. 2
Gradient descent update mechanism

Produce a randomly positive integer k (2  k  n). Based on the configuration N, we move the position of the amino acid k to new position k and amino acids k + 1, k + 2, ..., n are wholly translated as a rigid body according to the above-mentioned similar method we obtain a new configuration X . Thus, based on the configuration X, we obtain the new configuration X after applying two gradient methods.

In gradient descent update mechanism, the search has a tendency to fall into traps of local minima. ELP must spend a great deal of time to jump out of the traps. To enhance the efficiency of ELP and be able to sample also very narrow and deep valleys of the energy landscape, we executed alternately the gradient descent update mechanism and the following rotary update mechanism. After the simulation executed the gradient descent update mechanism 5,000 times, the rotary update mechanism was implemented 50 times.

Rotary update mechanism

The procedure of the rotary update mechanism is displayed in Fig. 3. To update the configuration X, we selected randomly a positive integer i (2 ≤ i  n). Since the ith amino acid lies on the circumference with radius unity around the (i  1)th amino acid, we selected randomly point i′ on circumference with maximum opening angle 2θmax (the dark area in Fig. 3). In our simulations of the AB model, we used a very small opening angle, cos(θmax) = 0.9, in order to be able to sample very narrow and deep energy basins. The ith amino acid was moved to the point i′ and the (i + 1)th, (i + 2)th, ..., nth amino acids were wholly translated as a rigid body with fixed trend of the partial chain. Thus an updated configuration X was obtained.

Fig. 3
Rotary update mechanism

Local search based on gradient descent

In ELP, the motivation for introducing the penalty term H(E, t) is to escape local minima. But, there is a technical flaw: the collection of the histogram function H(E, t) is performed by dividing the energy space into finite-size bins. Consider an attempted move in ELP which yields a new lower energy minimum that has never been visited before and happens to fall into the same bin containing other energies previously visited in the earlier steps. Undesirably, due to the penalty term H(E, t), the likelihood of accepting this new energy minimum becomes small. As a result, the ELP minimization may miss this new lower energy state near local minima. For this, we proposed the following local search method.

We introduced the additional degrees-of-freedom of bond lengths so that all monomers can freely move near the local minima derived by the ELP minimization. Suppose that there exists a spring with the original length 1 between the centers of the ith and (i + 1)th monomers (i = 1, 2, ..., n  1). According to Hooke’s law, the elastic potential energy between two adjacent monomers is proportional to the square of the length of the spring transformation. So the elastic potential energy of the system can be indicated as follows:

equation M9

where k is a physical coefficient characterizing the rigidities of all the springs. Thus, the sum of the potential energy is U = E + E . Obviously, E is the “penalty” term in the potential energy U. When the physical coefficient k is great enough, ri,i + 1 →1 (i = 1,2,...,n  1) and the optimal solution of the unconstrained optimization problem min U(X) is also the optimal solution of the optimization problem min E(X) with the constraint

equation M10
1

For the unconstrained optimization problem min U(X), we employed the gradient descent method to find low-energy configurations for a given monomer chain. In the beginning of the computation, we let k [set membership] [1,10]. Subsequently, we modified the physical coefficient k through multiplying by the factor 1.3 per 100,000 iterative steps. During the beginning of the calculation, the physical coefficient of all the springs was small, such that all monomers could move freely and easily attain low-energy states. Subsequently, as the calculation was carried out, the physical coefficient k increased gradually so as to increase the “penalty” and satisfy constraint (1) gradually, and finally the interaction of the springs disappeared. Obviously, when the physical coefficient rose to a large number (e.g., ≥1010), that is, when the springs turned rigid, constraint (1) was satisfied approximately, and a global energy minimum configuration was found, or a new trap of local minimum occurred, which had lower energy than the local minima derived by the ELP minimization.

Simulations

To search low-energy configurations for a given monomer chain, the ELP minimization was first executed. In our ELP simulations, the initial configuration was obtained through selecting the lowest-energy one from 10 randomly produced configurations. After the ELP was run up to 8 × 107 MC sweeps, the lowest-energy configuration was picked out from all energy states previously visited. By introducing the additional degree-of-freedom of bond lengths, the local search based on gradient descent was then used to search lower energy states near the lowest-energy state found by the ELP minimization. The gradient descent minimization was repeated up to five times. During the beginning of each round of gradient minimization, the physical coefficient was restored to initial value so that all monomers could freely move and easily attain low-energy states. After five rounds of gradient minimization, if the simulations did not find states with lower energy than the target value, we chose randomly another initial configuration for a new round of ELP simulation and subsequent local search based on gradient descent.

Results

In this article, we restricted ourselves to the 2D off-lattice AB model with Fibonacci sequences. The Fibonacci sequences are defined recursively by S0 = A, S1 = B, Si + 1 = Si  1 * Si. Here “*” is the concatenation operator. The first few sequences are S2 = AB, S3 = BAB, S4 = ABBAB, and so forth. They have the lengths given by ni + 1 = ni  1 + ni, i.e., given by the Fibonacci numbers. Following Refs. [6, 7, 913], in this article we considered the sequences with lengths n = 13, 21, 34 and 55 listed in Table 1.

Table 1
Test sequences and the lowest energies obtained by ELP+ in the 2D AB model, in comparison with those by MMC [7], PERM+ [10], ACMC [11], CSA [12], and STA [13], respectively

For each sequence, the ELP minimization was run 50 times independently on a PIV 2.0 GHz computer. After each round of ELP simulation, the additional degree-of-freedom of bond length was introduced, and the local search based on gradient descent was then used to search lower energy states near the lowest-energy state found by the ELP minimization. The calculated results showed that, along with the increment in the calculation steps in the AB model, all monomers ultimately tended to be stable and constraint (1) was satisfied approximately. The error margin was smaller than 10  10, i.e.,equation M17.

The lowest energies found by the ELP minimization and subsequent local searches based on gradient descent (ELP+) in 50 runs are listed in Table 1. We compare these results with minimum energies by conventional Metropolis-type Monte Carlo (MMC) procedures [7], the off-lattice pruned–enriched Rosenbluth method (PERM) and subsequent conjugate gradient minimization (PERM+) [10], the annealing contour Monte Carlo (ACMC) method [11], as well as conformational space annealing (CSA) [12] and statistical temperature annealing (STA) [13]. Table 1 shows that the results obtained by ELP+ are better than those of MMC and PERM+ for all four sequences. It can also be seen from Table 1 that ELP+ found slightly lower energy states for the 55-mer sequence than ACMC, CSA, and STA, while the putative ground-state energies found by ELP+ for 13-mer and 21-mer sequences are identical to those by ACMC, CSA, and STA. But ACMC, CSA and STA found lower energy states than ELP+ for the 34-mer sequence.

Figure 4 shows the lowest-energy configurations obtained by ELP+ for the four Fibonacci sequences listed in Table 1, where solid dots and empty circles indicate hydrophobic and hydrophilic monomers, respectively. The configurations for 13-mer and 21-mer sequences are identical to those by PERM+, ACMC, CSA, and STA, as can be expected from their almost identical ground-state energies in Table 1. This implies that all five methods find ground-state configurations reasonably well for these short chains. For 34-mer and 55-mer sequences, on the other hand, the lowest-energy configurations by ELP+ are different from those by PERM+, ACMC, CSA, and STA, as the energy difference between them is obvious.

Fig. 4
The lowest-energy configurations for the four sequences listed in Table 1, found by applying the ELP+ algorithm to the 2D AB model

Summary and anticipation

In this article, we introduced a global optimization method, energy landscape paving (ELP). We used ELP and subsequent local search (ELP+) to optimize low-energy configurations of a 2D off-lattice AB protein model. The numerical results showed that the proposed methods are very promising for finding the ground states of proteins. Of course, the chosen model still has some factors that are not very realistic. But this does not affect the fact that the AB model is suitable for benchmarking for the protein structure prediction problem.

Although ELP+ has powerful global optimization performance, it is still easy to get into traps of local minima as a result of random sampling. In future work, we hope to improve the sampling method in ELP, and apply it to all-atom models with realistic potentials to design various kinds of more high performance algorithms for the protein structure prediction problem.

Acknowledgements

This work was supported by Foundation of Nanjing University of Information Science and Technology and Qing Lan Project.

References

1. Shakhnovich, E.I.: Proteins with selected sequences fold into unique native conformation. Phys. Rev. Lett. 72(24), 3907–3911 (1994). doi:10.1103/PhysRevLett.72.3907 [PubMed]
2. Camacho, C.J., Thirumalai, D.: Kinetics and thermodynamics of folding in model proteins. Proc. Natl. Acad. Sci. U. S. A. 90(13), 6369–6372 (1993). doi:10.1073/pnas.90.13.6369 [PubMed]
3. Dill, K.A., Bromberg, S., Yue, K., Fiebig, K.M., Yee, D.P., Thomas, P.D., Chan, H.S.: Principles of protein folding: a perspective from simple exact models. Protein Sci. 4(4), 561–602 (1995) [PubMed]
4. Honeycutt, J.D., Thirumalai, D.: The nature of folded states of globular proteins. Biopolymers 32(6), 695–709 (1992). doi:10.1002/bip.360320610 [PubMed]
5. Fukugita, M., Lancaster, D., Mitchard, M.G.: Kinematics and thermodynamics of a folding heteropolymer. Proc. Natl. Acad. Sci. U.S.A. 90(13), 6365–6368 (1993). doi:10.1073/pnas.90.13.6365 [PubMed]
6. Stillinger, F.H., Head-Gordon, T., Hirshfeld, C.L.: Toy model for protein folding. Phys. Rev. E 48(2), 1469–1477 (1993). doi:10.1103/PhysRevE.48.1469 [PubMed]
7. Stillinger, F.H., Head-Gordon, T.: Collective aspects of protein folding illustrated by a toy model. Phys. Rev. E 52(32), 2872–2877 (1995). doi:10.1103/PhysRevE.52.2872 [PubMed]
8. Anfinsen, C.B.: Principles that govern the folding of protein chains. Science 181(96), 223–230 (1973). doi:10.1126/science.181.4096.223 [PubMed]
9. Yue, X.H., Tang, H.W., Guo, C.H.: A Tabu search and its application in 2D HP off-lattice model. Comput. Appl. Chem. 22(12). 1101–1105 (2005)
10. Hsu, H.P., Mehra, V., Grassberger, P.: Structure optimization in an off-lattice protein model. Phys. Rev. E 68(3), 037703 (2003). doi:10.1103/PhysRevE.68.037703 [PubMed]
11. Liang, F.: Annealing contour Monte Carlo algorithm for structure optimization in an off-lattice protein model. J. Chem. Phys. 120(14), 6756–6763 (2004). doi:10.1063/1.1665529 [PubMed]
12. Kim, S., Lee, S.B., Lee, J.: Structure optimization by conformational space annealing in an off-lattice protein model. Phys. Rev. E. 72(1), 011916 (2005). doi:10.1103/PhysRevE.72.011916 [PubMed]
13. Kim, J., Straub, J.E.: Structure optimization and folding mechanism of off-lattice protein models using statistical temperature molecular dynamics simulation: statistical temperature annealing. Phys. Rev. E 76(1), 011913 (2007). doi:10.1103/PhysRevE.76.011913 [PubMed]
14. Hansmann, U.H.E., Wille, L.T.: Global optimization by energy landscape paving. Phys. Rev. Lett. 88(6), 068105 (2002). doi:10.1103/PhysRevLett.88.068105 [PubMed]
15. Schug, A., Wenzel, W., Hansmann, U.H.E.: Energy landscape paving simulations of the trp-cage protein. J. Chem. Phys. 122(19), 194711 (2005). doi:10.1063/1.1899149 [PubMed]
16. Liu, J.F., Huang, W.Q.: Studies of finding low energy configurations in off-lattice protein models. J. Theor. Comput. Chem. 5(3), 587–594 (2006). doi:10.1142/S0219633606002453
17. Liu, J.F., Huang, W.Q.: Quasi-physical algorithm of an off-lattice model for protein folding problem. J. Comput. Sci. Technol. 22(4), 574–597 (2007). doi:10.1007/s11390-007-9067-x
18. Besold, G., Risbo, J., Mouritsen, O.G.: Efficient Monte Carlo sampling by direct flattening of free energy barriers. Comput. Mater. Sci. 15(3), 311–340 (1999). doi:10.1016/S0927-0256(99)00023-3
19. Cvijovic, D., Klinowski, J.: Taboo search: an approach to the multiple minima problem. Science 267(3), 664–666 (1995). doi:10.1126/science.267.5198.664 [PubMed]

Articles from Journal of Biological Physics are provided here courtesy of Springer Science+Business Media B.V.