the prediction of a protein structure from its primary sequence is one of the most interesting problems in computational biology ( 1 ) .
native proteins usually fold much too fast ( by at least tens of magnitude ) to involve an exhaustive search ( 2 ) .
it is a classical puzzle of the protein folding that biological proteins could not have originated from random sequences . despite a tremendous amount of efforts and
progresses over many decades , the problem remains essentially unsolved . in 1963 , anfinsen and his colleagues made a remarkable discovery that the amino acid sequence of a protein was fully sufficient to specify the molecule s ultimate 3d shape and biological activity ( 3 ) . for most single domain proteins ,
the information coded in the amino acid sequence is sufficient to determine the three - dimensional folded structure , which is the minimum free energy structure .
based on this theory , the protein is described by the complete list of the atoms in a molecule , with connectivities , bond lengths , angles , and force constants between all pairs of atoms ( 4 ) .
this all - atom model involves complex energy force and needs astronomical computational time . to understand the folding mechanism , it is useful to study simplified models such as the hydrophobic - hydrophilic model introduced by dill 5 .
, 6 .. lattice protein folding models have been playing important roles in theoretical studies of protein folding ( 7 ) . in these models ,
a protein is represented by a self - avoiding chain of beads placed on a discrete lattice , with two types of beads used to mimic hydrophobic and polar ( hp ) .
the advantage of hp lattice models is that they are simple enough to be amenable to thorough theoretical study , which can provide fruitful insights to feed back to or test against realistic models and experiments .
the biological foundation of this model is the believed theory that the first - order driving force of protein folding 8 .
, 9 . is due to a  hydrophobic collapse  in which those residues that prefer to be shielded from water ( hydrophobic residue ) are driven to the core of the protein , while those that interact more favorably with water ( polar residues ) remain on the outside of the protein .
previous papers researched the problem of protein folding on the cubic lattice model whose goal was to find the fold with the maximum number of contacts between non - covalently linked hydrophobic amino acids . yet , a significant drawback of the cubic lattice is the  parity problem  . in this paper , we present a triangular lattice model that overcomes the shortcomings of the cubic lattice model .
based on this model , we enumerated all possible compact structures and hp sequences . for a 333 cubic lattice model of size n=27 ,
the total number of compact structures is 103,346 and the number of all possible hp sequences is 2 . in the enumeration study , if the energy of every sequence folded into every compact structure is valued , the total number of evaluation will be 2103,3461.3910 . for the triangular models of size 3 + 4 + 5 + 4 + 3 and size 4 + 5 + 6 + 5 + 4 , the computation is 220,4861.0710 and 21,474,7822.4810
the algorithm decreased the computation by computing the objective energy of tree non - leaf nodes .
the parallel experiments proved that the fast tree search algorithm yielded an exponential speed - up factor of (1.486log2mlog2maxop ) , in which m is the number of different compact construct strings , maxop is the optimal string bound per leaf .
the methods currently used for the tertiary structure modeling are based on cubic lattice models to enumerate the minimization of the energy as a function of the topological contacts .
as previously stated , the quadrate lattice model exists a defect referred as the  parity problem  , in which only the residue in the even position of the primary sequence and another one in the odd site can form the topological contact .
the non - bonded neighbor ca nt be made between two even residues or two odd residues .
we generated and enumerated all the compact self - avoiding walks on the triangular lattice model of size 4 + 5 + 6 + 5 + 4 ( figure 1 ) .
the second is the two - mirror symmetry , of which one is the axial reflection and the other is the diagonal reflection .
so we got the differently directed 5,903,128/4=1,475,782 self - avoiding walks with this model . among the 1,475,782 paths
, we obtained reversal symmetries of 738,189 pairs and head - tail symmetries of 596 pairs ( figure 1 ) , so the number of different compact structures in this model was 738,189 + 596=738,785 .
it is energetically favorable for hydrophobic amino acids to occupy core sites , where there is low exposure to water . in this model
, we denoted a sequence of amino acids by i , and took only two types of amino acids , hydrophobic and polar .
the energy of a sequence folded into a structure was taken to be the sum of the contributions from each amino acid upon burial away from water : e=i=1nisi , where si is a structure - dependent number characterizing the degree of burial of the i amino acid in the chain .
larger si corresponds to a smaller surface area accessible to the solvent . for a structure on a 3d lattice ,
there are four different kinds of sites : center , face , edge , and corner .
therefore , in principle , there could be four different values of si . on a 2d lattice model
, we took only two values for si , and defined a string si for each structure with si = 1 if the i site is a core and si = 0 if it is a surface ( figure 1 ) .
out of the 1,475,782 compact structures , the number of distinct structure strings was 219,093 , among which there were 25,825 lattice conformations , and each represented exactly one structure .
each structure string has exactly ten 1 s and fourteen 0 s . our goal was to find a target structure string { sj } , which had unique and minimum energy corresponding to the target sequence string { i}.
it is obvious that the target structure string must possess a certain similarity with the target sequence string .
according to the observation , we organized the structure string into a binary tree and clustered similar structure strings into the same tree node ( figure 2 ) .
the algorithm decreased the computation by computing the objective energy of non - leaf nodes to locate the target structure string .
the distinct structure strings was 219,093 in the model of 4 + 5 + 6 + 5 + 4 .
there were thus 109,656 distinct strings that we kept in the calculation excluding the reversal structure strings .
each node of the tree represented a subset of these strings and maintained the following three kinds of information .
first , a structure string would have the value 1 at the i position if and only if all the strings corresponding to this node have 1 s at the i position .
second , a structure string would have the value 1 at the i position if and only if there is a table entry in this node that has a 1 or 0 at the i position .
third , each string in a node would have 1 s at some undecided positions . for each string , missing ones are the sum of these 1 s . by construction , each string has exactly ten 1 s , so the number of missing ones is equal to 10 minus the sum of known ones .
that is , missing ones are single integers no greater than 10 for each node .
we splited each node at the position that made two child nodes as tightly clustered as possible , and measured this clustering for each child as the entropies for each site , with the tightest clustering corresponding to the minimum entropy .
specifically , for each node we regarded the site of minimum entropy s of its set of structure strings as a branch point ( 10):s = min{(pilogpi+qilogqi)},in which pi is the probability of the i position being 1 , and qi is the probability of the i position being 0 .
those nodes partition the strings in the parent according to the value of the given position i : one child has the entire parent strings where i=1 and the other child has all of the strings where i=0 .
each leaf node at the end of the tree contains a small list of structure strings  in the following experiments we defined the parameter as max .
given a sequence string { i } and a node of the tree , we hoped to obtain the bounds of all structure strings represented by the node .
clearly , for all strings in the node the upper bound can be expressed by : bupper = min{i*(u+k),i*k+m}. given such a tree , where the root corresponds to all of the structure strings , the question here is how to search the tree .
we computed the upper bound of the sequence string according to previous formula , and called this the objective value .
if there are any leaf nodes that achieve this objective , go to next sequence . if not , repeat with the reversed version of the sequence strings . again
, if we achieve the objective value , go to next sequence . if not , decrease the objective value by 1 and try again .
given an objective value , we searched the tree as follows , starting at the root node . if the upper bound on the node indicates that objective value is unachievable , backdate and search other nodes .
if one string or none has been found that satisfies the objective , backdate and search other nodes , else if two strings have been found that achieve this objective , exit and transfer the next sequence string . if the current search node is not a leaf node , try each of the child nodes .
we first tried the child that matches { i}. the following experiments proved that only this step would yield a speed - up factor 1.782 in the model of size 4 + 5 + 6 + 5 + 4 .
using the fast search tree , we were able to completely enumerate all possible hp sequences and compact structures in the 4 + 5 + 6 + 5 + 4 triangular lattice model .
for every sequence , we rapidly computed the energy value of all the compact structures , found the structure with the minimal energy value , and recorded the minimal energy value and energy value of the first excited state ( the second minimal energy value ) .
the overall computation is highly parallelizable because each sequence can be done independently . in order to implement the calculation of ground states for all sequences in parallel , it is useful to divide the sequences into groups and use these groups as the unit of parallelism .
we performed our computations with all 2 + 2 hp sequences excluding reversed strings , which produced 129 groups .
these 129 groups were executed on the legend group deepcomp1800 -p4 xeon 2  ghz -myrinet/ 512 large array multiple processors , which is a collection of 24 computers , each containing two intel pentium pro microprocessors of 2  ghz .
every machine ran the linux operating system and had at least 512  mb of memory .
we divided the space of sequences into 129 groups fi ( 0  i < 129 ) .
num  denotes the sign of group with the initial value of 0 ; parameter  count  denotes the number of performed group with the initial value of 0 .
the operation of main processor p0 includes : ( 1 ) take out m groups f0 , f1 ,  , fm1 ,
send the groups to processor p1 , p2 ,  , pm , respectively ; num = m ; count=0 ; ( 2 ) perform the data of fnum group , count added by 1 and num added by 1 ; ( 3 ) receive the result of processors pi ( 1  i  m ) ; count added by 1 ; ( 4 ) if num < n1 ( n denotes the number of groups ) , send the data of fnum group to processor pi , else notify processor pi to exit ; num added by 1 ; ( 5 ) if num < n1 , transfer to step 2 ; if count equals to n , transfer to step 6 ; ( 6 ) collect the results of all the other processors pi ( 1  i  m ) and exit .
the operation of other processors pi ( 1  i  m ) includes receiving the group sent by the main processor p0 , performing the fast tree search , and sending the result to the main processor p0 .
each leaf node at the end of the tree contains a small list of structure strings .
we considered the maximum structure string number ( max ) per leaf as a variant in our experiment . the creating tree time ( tt ) , search time ( tc ) , and total time ( total = tt+tc ) through the experiment
the computation time showed in figure 3 is the cpu time and its unit is second .
it is obvious that search time and total time are the lowest when max is 8 .
the reason is that when creating the tree , every node needs to search all the structure strings and the number of nodes in the tree increases exponentially as max decreases , so the creating tree time increases exponentially . at the same time ,
futile search is eliminated by computing the objective value and aim structure string is located rapidly ; hence the search time decreases exponentially .
this is because when max is small , the leaf node contains less structure strings , but information in its parent node is enough to describe that ( figure 4 ) , so continuous bisection will increase the computation of objective energy value , which makes search time increase .
when max equals to the total number of structure strings , there is only a node in the tree and the situation equals to enumerate all the structure strings .
the experiments proved that the fast tree search algorithm yielded an exponential speed - up factor of (1.486log2mlog2maxop ) in the 4 + 5 + 6 + 5 + 4 triangular lattice model , in which m is the total number of distinct structure strings , and maxop is optimal max .
whenever a ground state structure string is found for a sequence , the reversed sequence necessarily has the reversed structure string as a ground state .
we excluded the reversed structure string and reversed sequence string in order to simplify and statistically analyze the process .
there were totally 109,656 such structure strings unrelated by rotational , reflection , or reverse labeling symmetries .
for a given sequence , the ground state structure is found by calculating the energy of all compact structures .
we completely enumerated all the ground states of all 2 + 2 possible sequences , and found that only 181,375 sequences have unique ground states , and then we calculated the designability of each compact structure .
there are structures that can be designed by an enormous number of sequences , and there are poor structures that can only be designed by a few sequences .
the number of structures decreases monotonically as the number of sequences increases ( figure 5 ) .
the result offers the evidence that the highly designable structures can tolerate more mutations during evolution because these structures can be designed by more different sequences . under different designability situations
, we obtained that the ratio of sequences with the ground state energy 10 , 9 , 8 , and 7 among 181,375 non - degeneracy sequences was 62.64% , 33.50% , 3.83% , and 0.32% , respectively .
the statistical relationship between the sequences with ground state energy 10 , 9 , 8 and the designability is shown in fig .
from the figure we can find that the sequences with the ground state energy 10 account for 95% in the structure of the highest designability , while those with the ground state energy 9 and 8 account for 5% and 0 , respectively ; the sequences with the ground state energy 10 account for 52% in the structure of the lowest designability , while those with the ground state energy 9 and 8 account for 37% and 11% , respectively .
this indicates that the compact structures of highly designability are , on average , thermodynamically more stable than other structures .
considering a structure string to be a chain of 0 s and 1 s linked by n1 links of three types , 0 - 0 , 1 - 0 or 0 - 1 , and 1 - 1 , with n00 , n10 or n01 , and n11 being the numbers of such links , respectively , we analyzed the relationship between designability and n01 ( figure 8) . when n01=10 , the structures of high designability occur more frequently ( figure 8a ) , and the relatively high frequency of low designability structures are due to the large number of structure strings with designability  1  ( ns=1 ) .
the tendency is approximately the same as n01=12 ( figure 8a ) , except lower frequency .
when n01= 8 , 9 , 11 ( figure 8b ) , the structures of high designability occur less frequently while the structures of low designability occur more frequently .
the tendencies of n01= 6 , 7 , 13 ( figure 8c ) are approximately the same as each other and have no much effect on the designability .
the occurrence of such kind of phenomenon is because that the structures of high designability have regular secondary structures . in the top twelve structures of high designability ( figure 9 )
, the number of n01 should be in a suit scope to form a regular second structure ; it is appropriate to form a regular helix structure when n01=10 .
there is no regular sheet in the structures of high designability because of the energy computation simplification of the lattice model .
we suppose that each dot in the 24d hypercube represents the sequence and structure in the lattice model . in the hp model ,
the energy of a sequence folded into a particular structure is the hamming distance between their binary strings .
hence , the number of sequences that fold uniquely to a particular structure  the designability of the structure  is the set of vertices lying closer to that structure than to any others .
it happens that in the hypercube the smallest hamming distance between two structures is approximately proportional to the difference in their respective n10 numbers .
this is evident in figure 10 , where the smallest hamming distance is plotted against the difference in n01 for all the pairs among the 25,825 binary structures on a 4 + 5 + 6 + 5 + 4 lattice , and is consistent with results given by li et al .
( 11 ) in which x(p ) ( the degree of clustering of hydrophobic residues ) is analogous to n01 . to see whether what we have observed so far has anything to do with real proteins
, we compared five sequences , p15 , each being a concatenation of a set of real proteins or ( 4 + 5 + 6 + 5 + 4 ) lattice binary peptides .
p1 , the representative non - redundant 350 proteins culled from protein data bank ( pdb ; www.rcsb.org/pdb/ ) ; p2 , the sections in p1 that fold into helices ; p3 , the sections in p1 that fold into sheets ; p4 , the 2,919 peptides mapped to the highest designabilities ; p5 , the 2,077 peptides mapped to the lowest designabilities .
firstly we defined the frequency distribution function as : fi(l)(m)=fi(l)fi(l)z , in which fi(l)(m ) denotes the frequency of the m binary word of length l occurring in sequence pi and there are 2 words of length l. fi(l ) denotes the average frequency of fi(l)(m):fi(l)=mfi(l)(m)2l.z=(m(fi(l)(m)fi(l))2)1/2 denotes the normalized frequency distribution function .
the pairwise overlaps were defined as : oij(l)=m=12lfi(l)(m)fj(l)(m),where i = 2 , 3 ; j = 4 , 5 ; l = 4 ~ 14 . the relationship between pairwise overlaps of different sequences
it is seen that p4 ( p5 ) is positively ( negatively ) correlated with p3 . for all values of l , the strongest correlation occurs between the model sequence of high designability ( p4 ) and the real protein sequence rich in sheets ( p3 ) .
the sequence of low designability is strongly correlated with the sequence rich in helix ( p2 ) , and the strong correlation occurs between the two model sequences of high ( p4 ) and low ( p5 ) designabilities .
the strong correlation between p2 and p5 is to some extent an artifact of the lattice model .
since we only considered the simple structure string , fractional compact structures will be washed out due to path degeneracy .
in the opinion of algorithm , we created the cluster tree by using the relationship between structure strings .
the algorithm decreased the computation by computing the objective energy of non - leaf nodes .
the parallel experiments proved that the fast tree search algorithm yielded an exponential speed - up in models of size 4 + 5 + 6 + 5 + 4 . in this paper , we have presented the two - dimensional triangular lattice model to study the designability of protein folding . through enumerating all the possible compact structures and hp sequences , we found that different compact structures have rather different designability .
the compact structures of high designability exhibit lower ground state energy , showing that these structures are , on average , thermodynamically more stable than other ones .
the research offers strong evidence that compact structures of high designability are more regular and geometrically symmetric . in the opinion of structure ,
the triangular lattice model has no parity problem and its surface is more similar to the natural proteins , therefore it is possible to get more information from the research of highly designable compact structures on triangular lattice models .
 using a triangular lattice model to study the designability of protein folding , we overcame the parity problem of previous cubic lattice model and enumerated all the sequences and compact structures on a simple two - dimensional triangular lattice model of size 4 + 5 + 6 + 5 + 4 . 
 we used two types of amino acids , hydrophobic and polar , to make up the sequences , and achieved 223 + 212 different sequences excluding the reverse symmetry sequences . 
 the total string number of distinct compact structures was 219,093 , excluding reflection symmetry in the self - avoiding path of length 24 triangular lattice model . 
 based on this model , we applied a fast search algorithm by constructing a cluster tree . 
 the algorithm decreased the computation by computing the objective energy of non - leaf nodes . 
 the parallel experiments proved that the fast tree search algorithm yielded an exponential speed - up in the model of size 4 + 5 + 6 + 5 + 4 . 
 designability analysis was performed to understand the search result . 