# Simple Models of the Protein Folding Problem

## 1 Introduction

The word “protein” originates from the Greek word proteios which means “of the first rank” . Indeed, proteins are building blocks and functional units of all biological systems. They play crucial roles in virtually all biological processes. Their diverse functions include enzymatic catalysis, transport and storage, coordinated motion, mechanical support, signal transduction, control and regulation, and immune response. A protein consists of a chain of amino acids whose sequence is determined by the information in DNA/RNA. There are 20 natural amino acids nature uses to make up proteins. These differ in size and other physical and chemical properties. The most important difference however, as far as the determination of the structure is concerned, is their hydrophobicity, i.e. how much they dislike water. An open protein chain, under normal physiological conditions, will fold into a three-dimensional configuration to perform its function. This folded functional state of the protein is called the native state. For single domain globular proteins which are our focus here, the length of the chain is of the order of 100 amino acids (from $`30`$ to $`400`$). Proteins of longer chains usually form multidomains each of which can usually fold independently. In Fig. 1 is shown a globular protein flavodoxin whose function is to transport electrons. Like most water soluble single domain globular proteins, it is very compact with a roughly rounded shape. The folded geometry of the chain can be best viewed in the cartoonish ribbon diagram (Fig. 1c) of the backbone configuration (Fig. 1b). One can see that the geometry of this protein structure is several $`\alpha `$-helices sandwiching an $`\beta `$-sheet. The folded geometries, often referred to as folds, of proteins usually look far more regular than random, typically possessing secondary structures (e.g., $`\alpha `$-helices and $`\beta `$-sheets) and sometime even having tertiary symmetries. (One can recognize an approximate mirror symmetry in Fig. 1c.) One of the main goals of the protein folding problem is to predict the three-dimensional folded structure for a given sequence of amino acids.

The protein folding problem is the kind of biological problem that has an immediate appeal to physicists. A protein can fold (to its native state) and unfold (to a flexible open chain) reversibly by changing the temperature, pH, or the concentration of some denaturant in solution. While the study of denaturation of proteins can be traced back at least 70 years when Wu  pointed out that denaturation was in fact the unfolding of the protein from “the regular arrangement of a rigid structure to the irregular, diffuse arrangement of the flexible open chain”. A turning point was the work of Anfinsen on the so-called “thermodynamic hypothesis” in the late 50s and early 60s. Anfinsen  and later many others demonstrated that for single domain proteins 1) the information coded in the amino acid sequence of a protein completely determines its folded structure, and 2) the native state is the global minimum of the free energy. These conclusions should be somewhat surprising to physicists. For the configurational “(free) energy landscape” of a heteropolymer of the size of a protein is typically “rough”, in the sense that there are typically many metastable states some of which have energies very close to the global minimum. How could a protein always fold into its unique native state with the lowest energy? The answer is evolution. Indeed, random sequences of amino acids are usually “glassy” and usually can not fold uniquely. But natural proteins are not random sequences. They are a small family of sequences, selected by nature via evolution, that have a distinct global minimum well separated from other metastable states (Fig. 2). One might ask: what are the unique and yet common properties of this special ensemble of proteinlike sequences? In other words, can one distinguish them from other sequences without the arguably impossible task of constructing the entire energy landscape? The answer lies in the heart of the question we introduce in the next paragraph and is the focus of this discussion.

There are about 100,000 different proteins in the human body. The number is much larger if we consider all natural proteins in the biological world. Protein structures are classified into different folds. Proteins of the same fold have the same major secondary structures in the same arrangement with the same topological connections , with some small variations typically in the loop region. So in some sense, folds are distinct templates of protein structures. Proteins with a close evolutionary relation often have high sequence similarity and share a common fold. What is intriguing is that common folds occur even for proteins with different evolutionary origins and biological functions. The number of folds is therefore much lower than the number of proteins. Shown in Fig. 3 is the cumulative number of solved protein domains  along with the cumulative number of folds as a function of the year. It is increasingly less likely that a newly solved protein structure would take a new fold. It is estimated that the total number of folds for all natural proteins is only about 1000 . Some of the frequently observed folds, or “superfolds” , are shown in Fig. 4. Among apparent features of these folds are secondary structures, regularities, and symmetries. Therefore, as in the case of sequences, protein structures or folds are also a very special class. One might ask: Is there anything special about natural protein folds–are they merely an arbitrary outcome of evolution or is there some fundamental reason behind their selection ? Is the selection of protein structures coupled with the selection of protein sequences? We will now address these questions via a thorough study of simple models.

## 2 Simple Models and the Designability

The dominant driving force for protein folding is the so-called hydrophobic force . The 20 amino acids differ in their hydrophobicity and can be very roughly classified into two groups: hydrophobic and polar . Hydrophobic amino acids have greasy side chains made of hydrocarbons and like to stick together in water to minimize their contact with water. Polar amino acids have polar groups (with oxygen or nitrigen) in their side chains and do not mind so much to be in contact with water. The simplest model of protein folding is the so-called “HP lattice model” , whose structures are defined on a lattice and whose sequences take only two “amino acids”: H (hydrophobic) and P (polar) (see Fig. 5). The energy for a sequence folded into a structure is simply given by the short range contact interactions

$$H=\underset{i<j}{}e_{\nu _i\nu _j}\mathrm{\Delta }(𝐫_i𝐫_j),$$
(1)   
where $`\mathrm{\Delta }(𝐫_i𝐫_j)=1`$ if $`𝐫_i`$ and $`𝐫_j`$ are adjoining lattice sites but $`i`$ and $`j`$ are not adjacent in position along the sequence, and $`\mathrm{\Delta }(𝐫_i𝐫_j)=0`$ otherwise. Depending on the types of monomers in contact, the interaction energy $`e_{\nu _i\nu _j}`$ will be $`e_{\mathrm{HH}}`$, $`e_{\mathrm{HP}}`$, or $`e_{\mathrm{PP}}`$, corresponding to H-H, H-P, or P-P contacts, respectively (see Fig. 5) . We  choose these interaction parameters to satisfy the following physical constraints: 1) compact shapes have lower energies than any non-compact shapes; 2) H monomers are buried as much as possible, expressed by the relation $`e_{\mathrm{PP}}>e_{\mathrm{HP}}>e_{\mathrm{HH}}`$, which lowers the energy of configurations in which Hs are hidden from water; 3) different types of monomers tends to segregate, expressed by $`2e_{\mathrm{HP}}>e_{\mathrm{PP}}+e_{\mathrm{HH}}`$. Conditions 2) and 3) were derived from the analysis  of the real protein data contained in the Miyazawa-Jernigan matrix  of inter-residue contact energies between different types of amino acids. Since we consider only the compact structures all of which have the same total number of contacts, we can freely shift and rescale the interaction energies, leaving only one free parameter. Throughout this section, we choose $`e_{\mathrm{HH}}=2.3`$, $`e_{\mathrm{HP}}=1`$ and $`e_{\mathrm{PP}}=0`$ which satisfy conditions 2) and 3) above. The results are insensitive to the value of $`e_{\mathrm{HH}}`$ as long as both these conditions are satisfied. (The analysis in Ref.  on the interaction potential of amino acids arrived at a form

$$e_{\mu \nu }=h_\mu +h_\nu +c(\mu ,\nu ),$$
(2)   
where $`h_\mu `$ is the hydrophobicity of the amino acid $`\mu `$ and $`c`$ is a small mixing term. The additive term, i.e. the hydrophobic force, dominates the potential. The choice of $`e_{\mathrm{HH}}=2.3`$ in our study can be viewed as a result of a hydrophobic part -2 plus a small mixing part -0.3. Several authors have investigated the effect of the mixing contribution as a small perturbation to the additive potential .)

We have studied the model (1) on a three dimensional cubic lattice and on a two dimensional square lattice . For the three dimensional case, we analyze a chain composed of 27 monomers. We consider all the structures which form a compact $`3\times 3\times 3`$ cube. There are a total of 51,704 such structures unrelated by rotational, reflection, or reverse labeling symmetries . For a given sequence, the ground state structure is found by calculating the energies of all compact structures. We completely enumerate the ground states of all $`2^{27}`$ possible sequences. We find that only $`4.75\%`$ of the sequences have unique ground states and thus are potential proteinlike sequences. We then calculate the designability of each compact structure. Specifically, we count the number of sequences $`N_S`$ that have a given compact structure $`S`$ as their unique ground state. We find that compact structures differ drastically in terms of their designability, $`N_S`$. There are structures that can be designed by an enormous number of sequences, and there are “poor” structures which can only be designed by a few or even no sequences. For example, the top structure can be designed by $`3,794`$ different sequences ($`N_S=3,794`$), while there are $`4,256`$ structures for which $`N_S=0`$. The number of structures having a given $`N_S`$ decreases monotonically (with small fluctuations) as $`N_S`$ increases (Fig. 6a). There is a long tail to the distribution. Structures contributing to the tail of the distribution have $`N_S>>\overline{N_S}=61.7`$, where $`\overline{N_S}`$ is the average number. We call these structures “highly designable” structures. The distribution is very different from the Poisson distribution (also shown in Fig. 6a) that would result if the compact structures were statistically equivalent. For a Poisson distribution with a mean $`\overline{N_S}=61.7`$, the probability of finding even one structure with $`N_S>120`$ is $`1.76\times 10^6`$.

The highly designable structures are, on average, thermodynamically more stable than other structures. The stability of a structure can be characterized by the average energy gap $`\overline{\delta _S}`$, averaged over the $`N_S`$ sequences that design the structure. For a given sequence, the energy gap $`\delta _S`$ is defined as the minimum energy difference between the ground state energy and the energy of a different compact structure. We find that there is a marked correlation between $`N_S`$ and $`\overline{\delta _S}`$ (Fig. 6b). Highly designable structures have average gaps much larger than those of structures with small $`N_S`$, and there is a sudden jump in $`\overline{\delta _S}`$ for structures with $`N_S^c1,400`$. This jump is a result of two different kinds of excitations a ground state could have. One is to break an H-H bond and a P-P bond to form two H-P bonds, which has an (mixing) energy cost of $`2E_{\mathrm{HP}}E_{\mathrm{HH}}E_{\mathrm{PP}}=0.3`$. The other is to change the position of an H-mer from relatively buried to relatively exposed, so the number of H-water bonds (the lattice sites outside the $`3\times 3\times 3`$ cube are occupied by water molecules) is increased. This kind of excitations has an energy $`1`$. The jump in Fig. 6b indicates that the lowest excitations are of the first kind for $`N_S<N_S^c`$, but are a mixture of the first and the second kind for $`N_S>N_S^c`$.

A striking feature of the highly designable structures is that they exhibit certain geometrical regularities that are absent from random structures and are reminiscent of the secondary structures in natural proteins. In Fig. 7 is shown the most designable structure along with a typical random structure. We examined the compact structures with the 10 largest $`N_S`$ values and found that all have parallel running lines folded in a regular manner.

We have also studied the model on a 2D lattice. We take sequences of length 36 and fold them into compact $`6\times 6`$ structures on the square lattice. There are $`28,728`$ such structures unrelated by symmetries including the reverse-labeling symmetry. In this case, we did not enumerate all $`2^{36}`$ sequences but randomly sampled them to the extend where the histogram for $`N_S`$’s reached a reliable distribution. Similar to the 3D case, the $`N_S`$’s have a very broad distribution (Fig. 8a). In this case the tail decays more like an exponential. The average gap also correlates positively with $`N_S`$ (Fig. 8b). Again similar to the 3D case, we observe that the highly designable structures in 2D also exhibit secondary structures. In the 2D $`6\times 6`$ case, as the surface-to-interior ratio approaches that of real proteins, the highly designable structures often have bundles of pleats and long strands, reminiscent of $`\alpha `$ helices and $`\beta `$ strands in real proteins; in addition, some of the highly designable structures have tertiary symmetries (Fig. 9).

To ensure that the above results are not an artifact of the HP model, we have studied model (1) with 20 amino acids . In this case the interaction energies $`e_{\nu _i\nu _j}`$, where $`\nu _i`$ can now be any one of the 20 amino acids, are taken from the Miyazawa-Jernigan matrix –an empirical potential between amino acids. For the 3D $`3\times 3\times 3`$ system and the 2D $`6\times 6`$ system, the total numbers of sequences are $`20^{27}`$ and $`20^{36}`$, respectively, which are impossible to enumerate. So we randomly sampled the sequence space. Similar to the case of the HP model, $`N_S`$’s have a broad distribution in both 3D and 2D cases. Furthermore, the $`N_S`$’s correlate well with the ones obtained from the HP model (see Fig. 10). Thus the highly designable structures in the HP model are also highly designable in the 20-letter model . With 20 amino acids, there are few sequences that will have exactly degenerate ground states. For example, in the case of $`3\times 3\times 3`$ about $`96.7\%`$ of the sequences have unique ground states. However, many of these ground states are almost degenerate, in the sense that there are compact structures other than the ground state with energies very close to the ground state energy. If we require that for a ground state to be truly unique there should be no other states of energies within $`g_c`$ from the ground state energy, then the percentage of the sequences that have unique ground states is reduced to about $`30\%`$ and $`8\%`$ for $`g_c=0.4k_B`$T and $`g_c=0.8k_B`$T, respectively.

## 3 A Geometrical Interpretation

A number of questions arise: Among the large number of structures, why are some structures highly designable? Why does designability also guarantee thermodynamic stability? Why do highly designable structures have geometrical regularities and even symmetries? In this section we address these questions by using a geometrical formulation of the protein folding problem .

As we have mentioned before, the dominant driving force for protein folding is the hydrophobicity, i.e. the tendency for hydrophobic amino acids to hide away from water. To model only the hydrophobic force in protein folding, one can assign parameters $`h_\nu `$ to characterize the hydrophobicities of each of the 20 amino acids . Each sequence of amino acids then has an associated vector $`𝐡=(h_{\nu _1},h_{\nu _2},\mathrm{},h_{\nu _i},\mathrm{},h_{\nu _N})`$, where $`\nu _i`$ specifies the amino acid at position $`i`$ of the sequence. The energy of a sequence folded into a particular structure is taken to be the sum of the contributions from each amino acid upon burial away from water:

$$H=\underset{i=1}{\overset{N}{}}s_ih_{\nu _i},$$
(3)   
where $`s_i`$ is a structure-dependent number characterizing the degree of burial of the $`i`$-th amino acid in the chain. Eq. (3) is essentially a solvation model  at the residue level and can also be obtained by taking the mixing term of Eq. (2) to zero .

To simplify the discussion, let us consider only compact structures and let $`s_i`$ take only two values: 0 and 1, depending on whether the amino acid is on the surface or in the core of the structure, respectively. Therefore, each compact structure can be represented by a string $`\{s_i\}`$ of 0s and 1s: $`s_i=0`$ if the $`i`$-th amino acid is on the surface and $`s_i=1`$ if it is in the core (see Fig. 11a for an example on a lattice). Let us make further simplification by using only two amino acids: $`\nu _i=\mathrm{H}`$ or P, and let $`h_\mathrm{H}=1`$ and $`h_\mathrm{P}=0`$. Thus a sequence $`\{\nu _i\}`$ is also mapped into a string $`\{\sigma _i\}`$ of 0s and 1s: $`\sigma _i=1`$ if $`\nu _i=\mathrm{H}`$ and $`\sigma _i=0`$ if $`\nu _i=\mathrm{P}`$. Let us call this model the PH (Purely Hydrophobic) model. Assuming every compact structure of a given size has the same numbers of surface and core sites and noting that the term $`_i\sigma _{i}^{}{}_{}{}^{2}`$ is a constant for a fixed sequence of amino acids and does not play any role in determining the relative energies of structures folded by the sequence, Eq. (3) is then equivalent to :

$$H=\underset{i=1}{\overset{N}{}}(\sigma _is_i)^2.$$
(4)   
Therefore, the energy for a sequence $`\stackrel{}{\sigma }=\{\sigma _i\}`$ folded onto a structure $`\stackrel{}{s}=\{s_i\}`$ is simply the distance squared (or the Hamming distance in the case where both $`\{\sigma _i\}`$ and $`\{s_i\}`$ are strings of 0s and 1s) between the two vectors $`\stackrel{}{\sigma }`$ and $`\stackrel{}{s}`$.

We can now formulate the designability question geometrically. We have two ensembles or spaces: one being all the sequences $`\{\stackrel{}{\sigma }\}`$ and the other all the structures $`\{\stackrel{}{s}\}`$. Both are represented by $`N`$-dimensional points or vectors where $`N`$ is the length of the chain. The points of all the sequences are trivially distributed in the $`N`$-dimensional space. In the case of the PH model, the points representing sequences are all the vertices of an $`N`$-dimensional hypercube (all possible $`2^N`$ strings of 0s and 1s of length $`N`$). On the other hand, the points representing all the structures $`\{\stackrel{}{s}\}`$ have a very different distribution in the $`N`$-dimensional space. The $`\stackrel{}{s}`$’s are constrainted and correlated. For example, in the case of the PH model where $`s_i=0`$ or 1, not every string of 0s and 1s actually represents a structure. In fact, only a very small fraction of the $`2^N`$ strings of 0s and 1s correspond to structures. If we consider only compact structures where $`_is_i=n_c`$ with $`n_c`$ the number of core sites, then the structure vectors $`\{\stackrel{}{s}\}`$ cover only a small fraction of the vertices of a hyperplane in the $`N`$-dimensional hypercube.

Now imagine putting all the sequences $`\{\stackrel{}{\sigma }\}`$ and all the structures $`\{\stackrel{}{s}\}`$ together in the $`N`$-dimensional space (see Fig. 12 for a schematic illustration). (In a more general case it would be simplest to picture if one normalizes $`\{\stackrel{}{h}\}`$ and $`\{\stackrel{}{s}\}`$ so that $`0h_i,s_i1`$.) From Eq. (4), it is evident that a sequence will have a structure as its unique ground state if and only if the sequence is closer (measured by the distance defined by Eq. (4)) to the structure than to any other structures. Therefore, the set of all sequences $`\{\stackrel{}{\sigma }(\stackrel{}{s})\}`$ that uniquely design a structure $`\stackrel{}{s}`$ can be found by the following geometrical construction: Draw bisector planes between $`\stackrel{}{s}`$ and all of its neighboring structures in the $`N`$-dimensional space (see Fig. 12). The volume enclosed by these planes is called the Voronoi polytope around $`\stackrel{}{s}`$. $`\{\stackrel{}{\sigma }(\stackrel{}{s})\}`$ then consists of all sequences within the Voronoi polytope. Hence, the designabilities of structures are directly related to the distribution of the structures in the $`N`$-dimensional space. A structure closely surrounded by many neighbors will have a small Voronoi polytope and hence a low designability; while a structure far away from others will have a large Voronoi polytope and hence a high designability. Furthermore, the thermodynamic stability of a folded structure is directly related to the size of its Voronoi polytope. For a sequence $`\stackrel{}{\sigma }`$, the energy gap between the ground state and an excited state is the difference of the squared distances between $`\stackrel{}{\sigma }`$ and the two states (Eq. (4)). A larger Voronoi polytope implies, on average, a larger gap as excited states can only lie outside of the Voronoi polytope of the ground state. Thus, this geometrical representation of the problem naturally explains the positive correlation between the thermodynamic stability and the designability.

As a concrete example, we have studied a 2D PH model of $`6\times 6`$ . For each compact structure, we devide the 36 sites into 16 core sites and 20 surface sites (see Fig. 11a). Among the $`28,728`$ compact structures unrelated by symmetries, there are $`119`$ that are reverse-labeling symmetric. (For a reverse-labeling symmetric structure, $`s_i=s_{N+1i}`$.) So the total number of structures a sequence can fold onto is $`(28,728119)\times 2+119=57,337`$, which map into $`30,408`$ distinct strings. There are cases in which two or more structures map into the same string. We call these structures degenerate structures and a degenerate structure can not be the unique ground state for any sequence in the PH model. Out of the $`28,728`$ structures, there are $`9,141`$ nondegenerate structures (or $`18,213`$ out of $`57,337`$). A histogram for the designability of nondegenerate structures is obtained by sampling the sequence space using 19,492,200 randomly chosen sequences and is shown in Fig. 11b. The set of highly designable structures are essentially the same as those obtained from the HP model discussed in the previous section. To further probe how structure vectors are distributed in the $`N`$-dimensional space, we measure the number of structures, $`n_\stackrel{}{s}(d)`$, at a Hamming distance $`d`$ from a given structure $`\stackrel{}{s}`$. Note that all the $`57,337`$ structures are distributed on the vertices of the hyperplane defined by $`_is_i=16`$. There are a total of $`C_{36}^{16}=7,307,872,110`$ vertices in the hyperplane. If the structure vectors were distributed uniformly on these vertices, $`n_\stackrel{}{s}(d)`$ would be the same for all structures and would be: $`n^0(d)=\rho N(d)`$, where $`\rho =57,337/7,307,872,110`$ is the average density of structures on the hyperplane and $`N(d)=C_{16}^{d/2}C_{20}^{d/2}`$ is the number of vertices at distance $`d`$ from a given vertex. In Fig. 13a, $`n_\stackrel{}{s}(d)`$ is plotted for three different structures with low, intermediate, and high designabilities, respectively, along with $`n^0(d)`$. We see that a highly designable structure typically has fewer neighbors than a less designable structure, not only at the smallest $`d`$s but out to $`d`$s of order 10-12. Also, $`n_\stackrel{}{s}(d)`$ is considerably larger than $`n^0(d)`$ for small $`d`$ for structures with low designability. These results indicate that the structures are very nonuniformly distributed and are clustered–there are highly populated regions and lowly populated regions. A quantitative measure of the clustering environment around a structure is the second moment of $`n_\stackrel{}{s}(d)`$,

$$\gamma ^2(\stackrel{}{s})=<d^2><d>^2=4\underset{ij}{}s_is_jc_{ij},$$
(5)   
where

$$c_{ij}=<s_is_j><s_i><s_j>$$
(6)   
and $`<>`$ denotes average over all structures. In Fig. 14a we plot the designability $`N_S`$ of a structure vs. its $`\gamma `$. We see that while larger $`N_S`$ implies smaller $`\gamma `$ it is not true vice versa. This is because that $`N_S`$ is very sensitive to the local environment at small $`d`$s while $`\gamma `$ is more a global measure.

What are the geometrical characteristics of the structures in the highly populated regions and lowly populated regions, respectively? This is something we are very interested in but know very little about. Naively, the structures in the highly populated regions are typical random structures which can be easily transformed from one to another by small local changes. On the other hand, structures in lowly populated regions are “atypical” structures which tend to be more regular and “rigid”. They have fewer neighbors so it is harder to transform them to other structures with only small rearrangements. One geometrical feature of highly designable structures is that they have more surface-to-core transitions along the backbone, i.e. there are more transitions between 0s and 1s in the structure string for a highly designable structure than average . We found a good correlation between the number of surface-core transitions in a structure string $`\stackrel{}{s}`$, $`T(\stackrel{}{s})`$, and $`\gamma (\stackrel{}{s})`$ (Fig. 13b). Thus, a necessary condition for a structure to be highly designable is to have a small $`\gamma `$ or a large $`T`$.

A great advantage of the PH model is that it is simple enough to test some ideas immediately. Two quantities often used to characterize structures are the energy spectra $`𝒩(E,\stackrel{}{s})`$  and $`𝒩(E,\stackrel{}{s},C)`$ . The first one is the energy spectrum of a given structure, $`\stackrel{}{s}`$ over all sequences, $`\{\stackrel{}{\sigma }\}`$:

$$𝒩(E,\stackrel{}{s})=\underset{\{\stackrel{}{\sigma }\}}{}\delta [H(\stackrel{}{\sigma },\stackrel{}{s})E].$$
(7)   
The second one is over all sequences of a fixed composition $`C`$ (e.g. fixed numbers of H-mers and P-mers in the case of two-letter code), $`\{\stackrel{}{\sigma }\}_C`$:

$$𝒩(E,\stackrel{}{s},C)=\underset{\{\stackrel{}{\sigma }\}_C}{}\delta [H(\stackrel{}{\sigma },\stackrel{}{s})E].$$
(8)   
It is easy to see that if two structure strings $`\{s_i\}`$ and $`\{s_i^{}\}`$ are related by permutation, i.e. $`s_i=s_{k_i}^{}`$, for $`i=1,2,\mathrm{},N`$, where $`k_1,k_2,\mathrm{},k_N`$ is a permutation of $`1,2,\mathrm{},N`$, then $`𝒩(E,\stackrel{}{s})=𝒩(E,\stackrel{}{s^{}})`$ and $`𝒩(E,\stackrel{}{s},C)=𝒩(E,\stackrel{}{s^{}},C)`$. Thus all maximally compact structures have the same energy spectra Eqs. (7) and (8). Therefore, structures differ in designability, not because they have different energy spectra Eqs. (7) and (8) , but because they have different neighborhood in the structure space.

## 4 Folding Dynamics and Thermodynamic Stability

Will highly designable structures also fold relatively fast? This question is addressed in detail in Ref.  (see also Ref. ). A quantity often used to measure how much a sequence is “proteinlike” is the $`Z`$ score,

$$Z=\frac{\mathrm{\Delta }}{\mathrm{\Gamma }},$$
(9)   
where $`\mathrm{\Delta }`$ is the average energy difference between the ground state and all other states and $`\mathrm{\Gamma }`$ is the standard deviation of the energy spectrum. $`Z`$ score was first introduced in the inverse folding problem  and later used in protein design . It has been shown in the context of the Random Energy Model  that $`Z`$ score is related to $`T_f/T_g`$ where $`T_f`$ is the folding temperature and $`T_g`$ the glass transition temperature . We have found a good and negative correlation between the folding time and the $`Z`$ score of the compact structure energy spectrum . In the context of the PH model (3), for a sequence $`\stackrel{}{h}`$ and its ground state $`\stackrel{}{s}`$,

$`\mathrm{\Delta }`$ $`=`$ $`{\displaystyle \underset{i}{}}h_i(s_i<s_i>),`$ (10)   
$`\mathrm{\Gamma }`$ $`=`$ $`\sqrt{{\displaystyle \underset{ij}{}}h_ih_jc_{i,j}},`$ (11)   
where $`c_{ij}`$ is given by Eq. (6). So in principle for every structure $`\stackrel{}{s}`$ one can maximize the $`Z`$ score with respect to $`\stackrel{}{h}`$ to get the “best” or “ideal” sequence for $`\stackrel{}{s}`$ that gives the highest $`Z`$ score, $`Z_S`$. It is however much easier to obtain a lower bound $`Z_{S}^{}{}_{}{}^{}`$ for $`Z_S`$ by letting $`\stackrel{}{h}=\stackrel{}{s}`$: $`Z_{S}^{}{}_{}{}^{}=\mathrm{\Delta }^{}/\mathrm{\Gamma }^{}`$ with

$`\mathrm{\Delta }^{}`$ $`=`$ $`{\displaystyle \underset{i}{}}(s_i^2s_i<s_i>),`$ (12)   
$`\mathrm{\Gamma }^{}`$ $`=`$ $`\gamma /2,`$ (13)   
where $`\gamma `$ is given by Eq. (5). In Fig. 14b, the $`\mathrm{\Delta }^{}`$ for all the $`6\times 6`$ compact structures are plotted against $`N_S`$ for the PH model. There is little if any correlation between $`N_S`$ and $`\mathrm{\Delta }^{}`$ for the $`6\times 6`$ PH model. Thus, correlations between $`N_S`$ and $`Z^{}`$ in this model come mainly from the one between $`N_S`$ and $`\mathrm{\Gamma }^{}=\gamma /2`$ (Fig. 14a). So a large $`Z^{}`$ is a necessary but not sufficient condition for a structure to have a large $`N_S`$.

## 5 Summary

We have demonstrated with simple models that structures are very different in terms of their designability and that high designability leads to thermodynamic stability, “proteinlike” structural motifs, and foldability. Highly designable structures emerge because of an asymmetry between the sequence and the structure ensembles. Our results are rather robust and have been demonstrated recently in larger lattice models  and in off-lattice models . A broad distribution of designability has also been found in RNA secondary structures . However, the set of all sequences designing a good structure, instead of forming a compact “Voronoi polytope” like in proteins, forms a “neutral network” percolating the entire space . It would be interesting to study the similarities and differences of the two systems. Finally, our picture indicates that the properties of the proteinlike sequences are intimately coupled to that of the proteinlike (i.e. the highly designable) structures; the picture unifies various aspects of the two special ensembles. It also suggests that understanding the emergence and properties of the highly designable structures is a key to the protein folding problem.

This work was done in collaboration with Hao Li, Ned Wingreen, Régis Mélin, Robert Helling, and Jonathan Miller. I am grateful to Jeannie Chen for her critical reading of the manuscript.