# Folding Lennard-Jones proteins by a contact potential

## I Introduction

One of the most challenging open problems in computational molecular biology is the prediction of a protein’s structure from its amino acid sequence – the so called protein folding problem. Assuming that the folded state is thermodynamically stable  the problem can be formulated as follows: for given sequence of amino acids and interactions between them, find the conformation of minimal energy. A successful prediction is marked by the experimental validation of the structure. Such kind of prediction is still unfeasible .

In order to tackle the protein folding problem from a theoretical point of view, in practice one has to choose 1) a representation of protein structure; 2) an approximation for the energy; 3) an algorithm to minimize the energy (i.e. to find the lowest energy structure for a given amino acid sequence). A conceptually straightforward approach is to perform detailed Molecular Dynamics simulations of fully detailed all atom models. However, numerical integration of Newton’s equations is unrealistic with present computers on the time-scale of the folding process. Even more problematic is the choice of the energy function. An imperfect parameterization of the classical effective interactions between atoms could correspond to an energy landscape with minima unrelated to the desired ones.

The true energy function which dictates the folding of proteins is unknown and one must resort to simple approximations. A fairly commonly used approximation, that of contact energies, was recently proved  to be too crude to successfully fold real proteins, whose native fold is stabilized by the ”true” potential function. It is not known whether this failure is caused by the extreme complexity of the true potential, or to the oversimplification inherent in the contact energy approximation. Resolving this issue is the main purpose of this work. To do this, we pose and study in detail a very simple question:

> Can artificial “proteins”, whose constituent residues interact by a Lennard-Jones pairwise potential, be folded successfully by a pairwise contact potential ?

The aim of such a study is to uncover some of the general problems which arise when one is forced to resort to energy approximations.

The Weizmann group has developed the contact map approach to protein folding. The contact map  of a protein with $`N`$ residues is an $`N\times N`$ matrix $`𝐒`$, whose elements are defined as $`S_{ij}`$=1 if residues $`i`$ and $`j`$ are in contact, and 0 otherwise. Contact between two residues can be defined in different ways; one is to consider two amino acids in contact when their two $`C_\alpha `$ atoms are closer than some threshold (they used 8.5 Å). Another definition is based on the minimal distance between two ions that belong to the two residues . The central premise of this approach is that the contact map representation has an important computational advantage. Changing a few contacts in a map induces rather significant large-scale coherent moves of the corresponding polypeptide chain .

In order to carry out the program they had to solve two problems:

1. Finding an efficient procedure to execute a search in contact map space 
2. Developing a test to determine whether the resulting maps correspond to physically realizable conformations 

Using these techniques, they have demonstrated that a the simple contact energy approximation is unsuitable to assign the lowest energy to the native state even for one protein.  This last result highlights the fact that the bottleneck in the successful prediction of protein structure is the choice of the approximation for the energy.

The Sissa group has pursued a different strategy developing minimalist off-lattice models of proteins. By studying toy protein models they were able to reproduce the essential features of the real folding process in terms of stability and accessibility (first of all the typical “funnel” structure for the energy landscape, that has been shown  to play a critical role). They found that the case of a $`C_\alpha `$ chain with 4 species of amino acids, equipped by a suitable design procedure, is effective to represent global features of the energy landscape . In their model, the potential is the simplest generic off-lattice potential that maintains chain connectivity, provides an attractive interaction between monomers and ensures self-avoidance: a polynomial potential between beads that are nearest neighbors along the chain (to represent the $`C_\alpha `$-$`C_\alpha `$ virtual bond) and a Lennard-Jones (LJ) potential among all other beads. This off-lattice toy model of proteins has been a suitable tool to test potential extraction techniques  and, on the other hand, it useful to study of the dynamical properties can provide important insights into the real mechanisms acting in real proteins.

In this work, we investigate whether the contact map representation is as a suitable approximation to the LJ model of proteins; that is, whether there exists a suitable choice of contact energy parameters, for which the ground states of a pairwise contact potential constitute a good approximation to the lowest energy states of the LJ potential. We also discuss whether dynamics in contact map space, used in combination with Molecular Dynamics can, at least in the simple case presented here, be successfully used to perform energy minimization.

## II Protein Design with Lennard-Jones Potentials

### A Definition of the Off-Lattice Model

In this first part of the study we represent a protein by a chain of $`N`$ monomers or beads, representing the $`C_\alpha `$ atoms of the amino acids. A configuration $`𝐫`$ of the chain is defined by the coordinates $`𝐫_1,\mathrm{},𝐫_N`$ of each $`C_\alpha `$ atom in three-dimensional continuous space. We consider effective pairwise energies between amino acids. The interaction potential $`V_{i,j}`$ between the monomers $`i`$ and $`j`$ is the sum of a covalent bond term plus a LJ term:

$$V_{i,j}=\delta _{i,j1}f(r_{ij})+(1\delta _{i,j1})V_{ij}^{LJ}$$
(1)   
where

$$V_{ij}^{LJ}=\eta _{ij}\left[(\frac{\sigma _{ij}}{r_{ij}})^{12}(\frac{\sigma _{ij}}{r_{ij}})^6\right],$$
(2)   
and $`r_{ij}=|𝐫_i𝐫_j|`$ are the interparticle distances.

For the energy bond function $`f(r_{ij})`$ we choose the expression:

$$f(x)=a(xd_0)^2+b(xd_0)^4,$$
(3)   
that is the fourth-order expansion of a generic and symmetric function of the distance $`x`$. The parameters $`a`$ and $`b`$ are taken to be $`1`$ and $`100`$ respectively. We add a quartic term to the usual quadratic one  because a plain harmonic potential could induce energy localization in some specific modes, significantly increasing the time needed for equilibration.

The parameter $`d_0`$ represents the equilibrium distance of the nearest neighbors along the chain. We set $`d_0=3.8`$, as this is the experimental value for the mean distance between nearest neighbor $`C_\alpha `$ atoms along the chain in real proteins, as taken from the Protein Data Bank.

The total energy of a chain in conformation $`𝐫`$ is defined as:

$$E=\underset{i=1}{\overset{N}{}}\frac{𝐩_i^2}{2}+\underset{i=1}{\overset{N}{}}\underset{j>i}{}V_{ij}.$$
(4)   
The first term is the classical kinetic energy of the chain, where the $`𝐩_i`$’s are the canonical momenta conjugate to the $`𝐫_i`$’s.

### B First step: generation of off-lattice native-like structures

Off-lattice there are an infinite number of conformations almost all of which are not designable (i.e. there is no sequence which has the conformation as its ground state). We aim at working with dynamically accessible, designable and compact structures. We present here an outline of the procedure; further details can be found in Ref. . The simplest way to obtain such configurations is to initially deal with homopolymer chains, setting the parameters $`\eta _{ij}=\eta `$ and $`\sigma _{ij}=\sigma `$ of Eq. (2) for all the beads in the chain. By cooling such a chain below the $`\theta `$ temperature It is known  that varying the temperature an homopolymer chain presents two different phases according to the dominance of attractive or repulsive interaction energy (as it can be with the LJ potential). The transition temperature is usually called $`\theta `$ temperature. we obtain maximally compact configurations.

We fix the parameters to the constant values $`\eta =40`$ and $`\sigma =6.5`$. The parameter $`\eta `$ determines the energy scale, while $`\sigma `$ determines the interaction length between monomers. Such a value for $`\sigma `$ ensures that, in practice, two monomers significantly attract each other for distances smaller then $`8`$ Å, as we show in §IV. Such distance threshold is usually used for the inter-amino acids bond, and is determined by the requirement that the average number of $`C_\alpha C_\alpha `$ contacts for each amino acid is roughly equal to the respective numbers obtained with the all-atom definition of contacts .

We used Molecular Dynamics (MD) (integrating Newton’s laws of motion on a computer) for simulating the kinetics of the chains. We employed an efficient and precise simplectic algorithm , in which one varies the energy density $`ϵ=E/N`$, which is related to the temperature . In order to find low energy structures of the chain, we used a combination of MD and slow cooling. Starting with different, randomly selected, initial conditions we find different low energy configurations for the chain, since the minimum energy state of a homopolymer is largely degenerate. For each MD simulation the chain is equilibrated at its initial energy during 8000 integration steps. We estimate the temperature of the system by computing the average value of the kinetic energy of the chain . Then, after 8000 steps, we slowly decrease the temperature by rescaling each component of the momenta by the same factor ( $`<1`$). We use the value $`0.8`$ as a cooling factor. Then the chain is again equilibrated at the new temperature and the procedure is iterated until very low temperatures are reached. At low temperature the chain is “trapped” in the compact configuration of minimum energy.

We performed $`6`$ times this algorithm and obtained $`6`$ different native-like structures for sequences of lengths $`N=30`$, whose coordinates will be denoted by $`𝐫_\alpha ^{}`$, ($`\alpha =1,\mathrm{},6`$).

We measure the difference between two structures of equal length N, using the root mean square (RMS) distance $`D`$:

$$D(𝐫_i,𝐫_i^{})=\sqrt{\frac{1}{N}\underset{i=1}{\overset{N}{}}(𝐫_i𝐫_i^{})^2}$$
(5)   
where one structure is translated and rotated to get a minimal $`D`$. The standard procedure, described in Ref.  was used. Two structures are considered different if their distance $`D`$ is larger than $`1`$ Å. The mean distance among our 6 different native structures is about $`5`$ Å.

### C Second step: protein design

At this point we need to select, for each chosen compact structure, a sequence having its energy minimum on this structure. Furthermore, it should be able to reach the selected structure in an accessible time. In order to generate such sequences we use the design procedure of Ref. . We use $`4`$ types of amino-acids, thus we deal with $`10`$ different parameters $`\eta _{ij}`$ ($`i,j=1,\mathrm{},4`$). We set these parameters by hand, with the constraint that the inequalities $`2\eta _{\alpha \beta }<eta_{\alpha \alpha }+eta_{\beta \beta },(\alpha \beta )`$, are satisfied. The following values have been used:   
$`\eta _{1,1}=40`$, $`\eta _{1,2}=\eta _{2,1}=30`$, $`\eta _{1,3}=\eta _{3,1}=20`$, $`\eta _{1,4}=\eta _{4,1}=17`$, $`\eta _{2,2}=25`$, $`\eta _{2,3}=\eta _{3,2}=13`$, $`\eta _{2,4}=\eta _{4,2}=10`$, $`\eta _{3,3}=5`$, $`\eta _{3,4}=\eta _{4,3}=2`$, $`\eta _{4,4}=1`$.

We fix an identical amino acid composition for all the sequences and the values for each matrix element $`\eta _{ij}`$. Then we select the optimal sequence for a target structure using a Monte Carlo algorithm in sequence space. Starting with a random sequence, at fixed composition, we perform random permutations accepting or rejecting new sequences with respect to the Boltzmann factor $`P=\mathrm{exp}(\mathrm{\Delta }E/k_BT)`$, where $`\mathrm{\Delta }E`$ is the energy variation due to permutation and $`T`$ is the “temperature” parameter of the Monte Carlo optimization scheme ($`k_B`$ is set equal to 1). By slowly lowering the parameter $`T`$ we select the sequence that has the lowest energy in the target conformation. In such a way we modify the energy landscape: out of the multidegenerate minimum energy scenario of the homopolymer we select a single minimum configuration lowering its energy by designing an appropriate sequence.

This method to select the sequence cannot guarantee that the target structures are indeed global energy minima. To check the obtained sequences we slowly cool (as in §II B) each selected sequence about 50 times (each time starting from a different initial condition). The cooling simulations find, as the lowest energy configurations, the same structures as before. Hence, these target structures are probably the lowest energy structures for the selected sequences, at least among the dynamically accessible ones, and this result is confirmed by the simulation in contact map space presented below. This repeated cooling process generates also a set of alternative, higher energy, metastable configurations (typically $`23`$ for each sequence) that, if perturbed, “decay” to the corresponding global minimum configuration.

## III Derivation of a set of pairwise contact energy parameters

In the pairwise contact approximation the energy is written as

$$E^{pair}(𝐚,𝐒,𝐰)=\underset{i<j}{\overset{N}{}}𝐒_{ij}w(a_i,a_j).$$
(6)   
We denote by $`𝐚`$ the sequence of amino acids, by $`𝐒`$ the conformation (represented by its contact map ) and by $`𝐰`$ the set of energy parameters. If there is a contact between residues $`i`$ and $`j`$, the parameter $`w(a_i,a_j)`$, which represents the energy gained by bringing amino acids $`a_i`$ and $`a_j`$ in contact, is added to the energy. Two amino acids are said to be in contact if their $`C_\alpha `$ atoms are closer than an upper threshold distance $`R_U`$ whose value will be discussed below.

We now explain in detail the approximation involved in Eq. (6). Denote by $`𝒞`$ a micro-state of the system. In general the micro-state is specified by the coordinates of all the atoms of the proteins and of the water molecules of the solvent. In the present case, we are considering a Lennard-Jones model of a $`C_\alpha `$ chain, so a micro-state is completely specified by the coordinates of the N monomers comprising the chain. Since many microscopic conformations share the same contact map S, it is appropriate to define a free energy $`(𝐚,𝐒)`$ associated with this sequence and map:

$$\mathrm{Prob}(𝐒)e^{(𝐚,𝐒)}=\underset{𝒞}{}e^{E(𝒞)/k_BT}\mathrm{\Delta }(𝒞,𝐒)$$
(7)   
where $`\mathrm{\Delta }(𝒞,𝐒)`$=1 if $`𝐒`$ is consistent with $`𝒞`$ and $`\mathrm{\Delta }=0`$ otherwise; i.e. $`\mathrm{\Delta }`$ is a “projection operator” which ensures that only those configurations whose contact map is S contribute to the sum (7). For real proteins $`E(𝒞)`$ is the unknown “true” microscopic energy, here it is the Lennard-Jones energy as of Eq. (1).

Since it is impossible to evaluate this exact definition of the free energy of a map, we resort to a phenomenological approach, guessing the form of $`(𝐚,𝐒)`$ that would have been obtained had the sum (7) been carried out. $`E^{pair}(𝐚,𝐒,𝐰)`$ of eq. (6) is the simplest approximation to the true free energy. To test the extent to which this approximate form is capable of stabilizing the native map of a protein against other non-native maps, we must specify the parameters $`𝐰`$.

To derive the set $`𝐰`$ of pairwise contact energy parameters we require that the following set of conditions be satisfied for all the proteins in the database 

$$E^{pair}(𝐚,𝐒^{},𝐰)<E^{pair}(𝐚,𝐒_\mu ,𝐰).$$
(8)   
Here $`𝐒^{}`$ is the contact map of the native conformation of sequence $`𝐚`$ and $`𝐒_\mu `$ ($`\mu =1,\mathrm{}N_D`$) are contact maps taken from a set of $`N_D`$ alternative conformations. The problem is to find a set $`𝐰`$ such that for each decoy $`\mu `$ the inequality (8) is satisfied. The perceptron learning technique  is used to search for contact energy parameters $`𝐰`$ for which the set of conditions (8) are satisfied. The set $`𝐰`$ that gives the optimal solution of this problem is equivalent to the maximization of the gap between the ground state and the first excited state, as it is shown in the Appendix.

### A Generation of alternative conformations

To obtain a large ensemble of decoys we performed dynamics in contact map space collecting maps along the trajectory. The method is presented elsewhere ; here we outline only the details that are relevant for the present implementation. Our algorithm is divided in four steps.

1. We start from an existing map and perform large scale “cluster” moves. At this stage, no attempt is made to preserve physicality. The contact map which is obtained by this procedure is typically uncorrelated to the starting one.
2. The resulting map is refined by using local moves.
3. We use the reconstruction algorithm  to restore physicality by projecting the map obtained from the second step onto the physical subspace.
4. We perform a further optimization by an energy minimization in real space using a standard Metropolis crankshaft technique .

Using this algorithm, after the choice of a suitable definition of contact –as detailed in § IV B and § IV C–, we generated a set of $`N_D`$=60000 alternative conformations, 10000 for each of the 6 sequences in the database. As discussed in Ref. , the contact maps that are obtained by contact map dynamics are uncorrelated and the $`N_D`$ decoys form a representative set of low energy competitors for the native state. The important requirement is to generate a large set of uncorrelated decoys, since a good energy function must stabilize the correct structure and destabilize all the others.

This selection of candidates was performed using the set of LJ parameters $`\eta _{ij}`$ as a first guess for the pairwise contact energy parameters $`𝐰`$.

## IV Square well approximation of the LJ potential

### A Truncation of the LJ potential

In the definition of contact one has to specify the upper threshold $`R_U`$ of the distance between two beads, below which are considered to be in contact. To get an estimate of a reasonable value of $`R_U`$ we studied the truncation of the LJ potential. Using MD, we analyzed the stability of the folded conformations against a cut in the tail of the LJ potential. We used a potential equal to the Lennard-Jones function for $`r`$ smaller than $`R_C`$, and equal to zero for $`r`$ larger than $`R_C`$:

$$V_{LJ}^{}(r_{ij})=\{\begin{array}{cc}\eta _{ij}((\frac{\sigma }{r_{ij}})^{12}(\frac{\sigma }{r_{ij}})^6)\hfill & \mathrm{if}r_{ij}<R_C\hfill \\ & \\ 0\hfill & \mathrm{if}r_{ij}>R_C\hfill \end{array}$$
(9)   
For a given $`R_C`$, we determined the ground state conformations of the 6 proteins in the database. This was done by performing $`n`$ runs of Molecular Dynamics minimization (by slow cooling, as described in §II B with typically $`n=50`$). For each protein, each minimization ended in a slightly different structure $`𝐫_\alpha ^\nu (R_C)`$ ($`\alpha =1,\mathrm{},6`$, $`\nu =1,n`$) We measured the RMS distance $`D(R_C)`$, averaged over the 6 proteins in the database, between the native conformations and the conformations found as minima of the potential of Eq. (9):

$$D(R_C)=\underset{\alpha =1}{\overset{6}{}}\underset{\nu =1}{\overset{n}{}}D(𝐫_\alpha ^{},𝐫_\alpha ^\nu (R_C)).$$
(10)   
In Fig. 1 we show the results for each protein, for several choices of the cutting length $`R_C`$. In the inset the average distance $`D(R_C)`$ with its standard deviation as a function of the cutting length $`R_C`$ are shown. Remarkably, it is possible to cut the LJ potential inside the attractive well down to $`R_C`$=8 Å, still keeping the average distance below 1 Å.

### B Derivation of the region of physicality

For any value of $`R_U`$ the contact map $`𝐒`$ of any chain conformation can be easily determined. In order to construct decoys, however, we must deal with the inverse problem – to get a chain conformation starting from some contact map $`𝐒`$. In order to apply the reconstruction procedure introduced in Ref.  for a possibly non physical contact map one has to specify another length; the hard core distance $`R_L`$, defined as the minimal allowed distance between two amino acids. This is essential since without such a lower cutoff we would be able to construct chain configurations for maps with much higher numbers of contacts than the native maps and hence would identify such dense decoys as physical. Since the contact energies tend to be negative such decoys will have much lower energies than the native map and will never be able to find a solution to the inequalities eq. (12).

In Fig. 2 we show the histogram of all the pairwise distances for the six native conformations considered in this study. Apart from the peak at 3.7 Å, due to consecutive residues along the chain, the distances range from 6.1 to 18.8 Å. The vertical line indicates 6.9 Å, a distance whose role will be discussed below. There are 2436 distances between non-consecutive beads, of which 49 are below 6.9 Å and 7 below 6.5 Å.

We now consider the plane $`(R_U,R_L)`$ and set about to determine the “region of physicality” which is that part of the plane whose points allow reconstruction of the contact maps that correspond to the six native folds. To understand the problem we are facing, one should bear in mind that once we chose a value for $`R_U`$, the contact map associated with a conformation is determined. When, however, we try to reconstruct one of these maps, we must also specify a value for $`R_L`$ and if the chosen value is too high, we may find that no corresponding chain configuration exists. In this case, according to definition, the native contact maps are non-physical. The gross features of such a region can be outlined by some simple arguments. First, since we must have $`R_U>R_L`$, the region of physicality must clearly lie below the line $`R_U=R_L`$. Second, if $`R_L`$ is too high, the contact maps are non physical and there must be a (possibly horizontal)line that bounds the physical region from above. The exact location of this line must be determined numerically. For $`R_L<6.1`$ it is clear that all the native contact maps are physical, but for $`R_L>6.1`$ Å the situation is not clear. It may be able to ”catch” the few “close” contacts (those with $`d<6.9`$ Å) by locally stretching the chain. We proved numerically that this is indeed possible up to about $`R_L<6.7`$ Å, as shown in Fig. 3a and in more detail in Fig. 3b. The physical region is the shaded trapezoid open on its right side; for points in this region we are able to reconstruct chain conformations for all the six native maps (derived from their ”true” chain conformations by the choice of $`R_U`$).

### C Derivation of the region of learnability

We turn now to determine a second region in the $`(R_U,R_L)`$ plane, that of “ learnability”. At a point $`(R_U,R_L)`$ in this region it is possible to find a set of energy parameters $`𝐰`$ such that each of the six native maps have the lower energy than all their respective decoys. As was done above for the region of physicality, it is possible to predict the general shape of the region of learnability a well. First, it is limited by the same bisecting line $`R_U=R_L`$. Second, there is a vertical line at $`R_U=6.1`$ Å, at the left of which all the contact maps are empty – we are not interested in such a case. Third, there is a second vertical line, at $`R_U=19`$, at the right of which all the contact maps are filled. Also this case is not interesting. Fourth, there must be a line (approximately horizontal) below which the energy parameters are unlearnable. We expect this since for small values of $`R_L`$ amino acids are allowed to be very close to each other. Very compact conformations are possible and since the interaction parameters are on the average attractive, such compact conformations have very low energy. The exact location of this line must be determined numerically. The result is shown in Figs. 3a and 3b. For each choice of $`R_U`$ and $`R_L`$, we generated by contact map dynamics a set of $`N_D=60000`$ decoys, 10000 for each sequence in the database. Using the perceptron algorithm, we determined if such set of decoys was learnable or not. The outcome of the procedure is the triangular shaded region in Fig. 3a, whose lower edge is shown in detail in Fig. 3b.

The most interesting result is that the region of physicality and the region of learnability do not overlap. It is not possible to choose $`(R_U,R_L)`$ such that the contact maps of the native LJ conformations in the database are physical and they can be stabilized by a pairwise contact potential.

Two observation must be made. First, the points labeled as “unlearnable” are rigorous, whereas those labeled as “learnable” are tentative – it is possible that by increasing the number of decoys also those points will become unlearnable. The lower edge of the region of learnability could actually be moved upwards. The same effect should be expected by increasing the number of proteins in the database. Second, the points labeled as “physical” are rigorous, whereas those labeled as “non physical” are tentative – we could not reconstruct the corresponding contact maps, but this failure could be due to our reconstruction algorithm. Thus, the upper edge of the region of physicality could be moved upwards. However, by increasing the number of proteins in the database we expect that this upper edge would move downwards. In the limit of a large number of proteins the latter effect could possibly dominate. In conclusion, we argue that the gap between the region of learnability and the region of physicality is expected to widen when many proteins are taken into account.

## V Folding in contact map space

Since we have proved that the two regions, of learnability and physicality of the native LJ contact maps do not overlap, we can conclude that the answer to the question posed in the Introduction is negative – the contact energy approximation is unable to stabilize the native folds of LJ chains. The contact map approximation is too crude also for extremely simplified potential as the LJ potential used in this work.

The situation is, however, more subtle. We can look for contact energies that give rise to approximate solutions of the folding problem: it is impossible recover the LJ native configurations using the contact map approximation, but we want investigate how ”close” to the solution the contact map approximation can lead.

To tackle this point, we selected a working point in the learnable region. We chose $`R_U`$=8 Å, according the the minimal value of $`R_C`$ determined in Sec. IVA and $`R_L=6.9`$ Å. As shown above, for this value of $`R_L`$, native maps are non physical. In other words, for a given native map there is no chain which can realize that contact map. We found, however, that for these $`R_L`$ and $`R_U`$ it is possible to reconstruct conformations whose maps are physical and typically differing by only a few contacts from the native map. Moreover, the corresponding 3D structures are, on the average, at a RMS distance of less that 2 Å from the native LJ conformations. We decided to choose $`R_L=6.9`$ Å, as the best possible compromise between solving Eq. (8) and our aim is that some physical maps of low energy are close to the native maps.

Now we turn to answer the following important question: is the set $`𝐰`$ obtained this way useful to fold LJ proteins? The first aspect of this general question concerns the feasibility of using contact map dynamics. That is, we ask:

> Is it possible to fold the sequences in the database using contact map dynamics?

Fully successful folding in contact map space corresponds to exact recovery of the native map. Since for $`R_L=6.9`$ Å and $`R_U=8.0`$ Å the native contact maps are non physical, such an exact solution is not within our reach. Nevertheless, the contact map dynamics could select configurations slightly different from the native ones.

For each of the 6 sequences we started from some random contact map and ran the contact map dynamics procedure for 1000 steps; this took about 3 cpu hours on our HP812 computer. We now analyze in detail our results.

As already observed in Ref. , the correlation between energy and distance from the native map (the “funnel”) is not always very strong when using $`E^{pair}`$. A case of quite good correlations, obtained by contact map dynamics for sequence N<sup>o</sup> 5 (among the six used in this work), is shown in Fig. 4. Since the lowest energy contact map found in the simulation might or might not coincide with the closest in distance to the native one, in general we can obtain only a short list of candidates for the native state.

For sequence N<sup>o</sup> 5 the physical map of lowest energy is, in fact, the closest in distance to the native one. Both maps are shown in Fig. 5. The energy $`E^{pair}`$ of the native map is -35.75 and the energy of the predicted one is -33.33. The RMS distance between two corresponding conformations, shown in Fig. 6, is 2.1 Å. We emphasize again the obvious fact that since our $`𝐰`$ is a solution of (12), the native map has lower energy than all the low-energy decoys that were used for learning. Indeed, in no case have we found by contact map dynamics maps lower in energy than the native one. We summarize in Table I the results obtained for the six sequences.

The Hamming distance $`D_H`$ between two contact maps $`𝐒`$ and $`𝐒^{}`$ is defined by

$$D_H=\underset{j>i}{}\frac{|S_{ij}S_{ij}^{}|}{N_c(𝐒)N_c(𝐒^{})},$$
(11)   
which counts the number of mismatches between maps $`𝐒`$ and $`𝐒^{}`$; $`N_c(𝐒)`$ and $`N_c(𝐒^{})`$ are the numbers of contacts in the two maps.

The second and complementary question we pose is

> Are the low energy conformations generated by contact map dynamics good starting points for Molecular Dynamics minimization?

We analyze again sequence N<sup>o</sup> 5. We randomly selected 10 conformations among the 100 of lowest $`E^{pair}`$ energy that were generated by contact map dynamics and reconstructed their corresponding conformations. These were then used as starting points for a Molecular Dynamics minimization of the $`E^{LJ}`$ energy. In Fig. 7 we present the energy and distance to the native conformation for the conformations obtained by MD. As can be seen, 8 trajectories ended up very close to the native state and 2 in conformations at RMS distances of 6 and 8 Å respectively. We note that these last two initial conformations had the largest RMS distance from the native state.

We conclude that the predictions performed in contact map space using the approximate energy $`E^{pair}`$ are, on the average, suitable starting points for a Molecular Dynamics minimization of the LJ energy. Moreover, by using Molecular Dynamics it is possible to correctly rank the predictions that were obtained by contact map dynamics.

## VI Conclusions

There are four different situations for which the contact pairwise approximation has been used:

1) When the true potential is a contact pairwise one, such as for lattice models with nearest-neighbor interactions, the energy $`E^{pair}`$ of Eq. (6) is, by definition, the exact (free) energy. Models of this type were investigated in two recent studies . In both works, a database of proteins was designed using $`E^{pair}`$, with a given choice of energy parameters $`𝐰^{true}`$. The database was then used to recover a set $`𝐰`$, which ideally could be identical to $`𝐰^{true}`$. Since the native states of the proteins in the database are the lowest energy conformations of $`E^{pair}(𝐰^{true})`$, Eq. (8) can be satisfied for all the possible conformations for all the sequences in the database.

2) When the true potential is unknown, such as for real proteins, Eq. (6) is only an approximation to the true free energy (7), and there is no guarantee that such an approximation would lead to good results. The conformations stabilized by the true energy function are not necessarily stabilized against all decoys by a pairwise contact energy function. For a particular protein, crambin, it has been shown  that, in fact, such a solution does not exist and the contact maps of lowest pairwise energy are not close to the native maps.

3) The “true” potential is a contact pairwise one, and it is used to design protein sequences on the structure of existing proteins. This case was considered in Ref. , where it was shown that it is unfeasible to design sequences, using a pairwise contact potential, on structures of existing proteins. This result implies that the native structures of existing proteins are “designed” in a very peculiar way by the possibly very complex potential used by Nature. A simple pairwise contact potential can neither be tuned to stabilize the true sequences (see point 2) nor used to perform protein design on the true structures.

4) The true potential is a Lennard-Jones pairwise one and it is used to design sequences on structures typically arising from Lennard-Jones interactions. The resemblance of such structures to those of real proteins appear to be reasonable , however, as discussed in this paper, the results obtained for approximate contact potentials are remarkably different. Here we have shown that it is not possible to stabilize the exact native folds. This result is somewhat surprising – the approximation seems to be quite reasonable, since both the contact potential and the LJ potential are characterized by a single length scale. However, we also showed that in this case an approximate way to solve the folding problem can be found. We showed that a dynamics in contact map space is a suitable tool to perform energy minimization, when a correct parameterization of the energy is possible. Contact map dynamics provides, for a good choice of contact parameters, not just a unique structure corresponding to the lowest energy state, but a set of candidates for it. When this list of predictions is used as starting configurations of a MD energy minimization, which uses the original ”true” LJ form of the energy, the correct ground states are recovered.

It remains to be investigated if a hybrid method, based on “fusion” of contact map dynamics and molecular dynamics, in a way that takes advantage of the their respective advantages (that in some sense are complementary) may be a powerful tool to fold real proteins as well.

## Acknowledgments

C.C expresses her gratitude for the hospitality to the Weizmann Institute where part of this work was carried out. We thank Devarajan Thirumalai for discussions and Ron Elber for discussing with us a similar approach, based on Eq. 8 (unpublished). This research was supported by grants from the Minerva Foundation, the Germany-Israel Science Foundation (GIF) and by a grant from the Israeli Ministry of Science. During the last three months C.C. has been supported by the NSF (Grant # 96-03839) and by the La Jolla Interfaces in Science program (sponsored by the Burroughs Wellcome Fund).

## Perceptron learning

We describe here the perceptron learning technique that was used.

We first show that for any conformation the condition Eq.(8) can be trivially expressed as

$$𝐰𝐱^\mu >0$$
(12)   
To see this just note that for any map $`𝐒_\mu `$ the energy (6) is a linear function of the 10 contact energies that can appear and it can be written as

$$E^{pair}(𝐚,𝐒_\mu ,𝐰)=\underset{c=1}{\overset{10}{}}N_c(𝐚,𝐒_\mu )w_c$$
(13)   
Here the index $`c=1,\mathrm{}10`$ labels the different contacts that can occur and $`N_c(𝐚,𝐒_\mu )`$ is the total number of contacts of type $`c`$ that actually appear for the sequence $`𝐚`$ in the map $`𝐒_\mu `$. The difference between the energy of this map and the native $`𝐒^{}`$ is, therefore,

$$\mathrm{\Delta }E_\mu =\underset{c=1}{\overset{10}{}}x_c^\mu w_c=𝐰𝐱_\mu $$
(14)   
where we used the notation

$$x_c^\mu =N_c(𝐚,𝐒_\mu )N_c(𝐚,𝐒^{})$$
(15)   
and $`𝐒^{}`$ is the native map.

Each candidate map $`𝐒_\mu `$ is represented by a vector $`𝐱^\mu `$ and hence the question raised above regarding stabilization of $`𝐒^{}`$ for a sequence $`𝐚`$ becomes

> Can one find a vector $`𝐰`$ such that condition (12) holds for all $`𝐱^\mu `$?

If such a $`𝐰`$ exists, it can be found by perceptron learning.

A perceptron is the simplest neural network . It is aimed to solve the following task. Given a set $`P`$ of patterns (also called input vectors, examples) $`𝐱^\mu `$, $`\mu =1,\mathrm{},P`$, find a vector $`𝐰`$ of weights, such that the condition

$$h_\mu =𝐰𝐱^\mu >0$$
(16)   
is satisfied for every example from the set. If such a $`𝐰`$ exists for the training set, the problem is learnable; if not, it is unlearnable. We assume that the vector of “weights” $`𝐰`$ is normalized,

$$𝐰𝐰=1$$
(17)   
The vector $`𝐰`$ is “learned” in the course of a training session. The $`P`$ patterns are presented cyclically; after presentation of pattern $`\mu `$ the weights $`𝐰`$ are updated according to the following learning rule:

$$𝐰^{}=\{\begin{array}{cc}\frac{𝐰+\eta 𝐱_\mu }{|𝐰+\eta 𝐱_\mu |}\hfill & \mathrm{if}𝐰𝐱_\mu <0\hfill \\ & \\ 𝐰\hfill & \mathrm{otherwise}\hfill \end{array}$$
(18)   
This procedure is called learning since when the present $`𝐰`$ misses the correct “answer” $`h_\mu >0`$ for example $`\mu `$, all weights are modified in a manner that reduces the error. No matter what initial guess for the $`𝐰`$ one takes, a convergence theorem guarantees that if a solution $`𝐰`$ exists, it will be found in a finite number of training steps. .

If the region of parameter space whose points are solutions is large, one is interested in the optimal solution. To obtain the optimal perceptron solution, called the maximal stability perceptron, we used the Krauth and Mezard algorithm . In this algorithm the condition (16) is replaced by

$$h_\mu =𝐰𝐱^\mu >c$$
(19)   
where $`c`$ is a positive number that should be made as large as possible.

At each time step the “worst” example $`𝐱^\nu `$ is identified, namely the one such that

$$h_\nu =𝐰𝐱^\nu =\underset{\mu }{\mathrm{min}}𝐰𝐱^\mu $$
(20)   
and the example $`\nu `$ is used to update the weights according to the rule (18). The field $`h_\nu (t)`$ keeps changing at each time step $`t`$ and the procedure is iterated until it levels off to its asymptote.

### Captions to the figures

Fig. 1. Upper part: Approximation of the Lennard-Jones potential by a contact potential.   
Lower part: Distance of the conformation of lowest energy found in a MD simulation, done with a potential truncated at $`R_C`$, to the true native state of the non truncated LJ interaction. Each dot represents a different run. In the inset the average $`D(R_C)`$ is shown. Down to $`R_C=8`$ Å the distance to the true native state is below 1 Å.

Fig. 2. Distance probability distribution between amino acids in the six native conformations in this study. 

Fig. 3. Regions of physicality and of learnability in the ($`R_U`$,$`R_L`$) plane.

Fig. 4. Correlation between energy $`E^{pair}`$ and RMS distance $`D`$ to the native state and between energy and the Hamming distance $`D_H`$ (see Eq. (11)). Results refer to contact map dynamics.

Fig. 5. Contact maps of sequence N<sup>o</sup> 5: native (below diagonal) and of lowest energy found during the simulation (above diagonal).

Fig. 6. Native conformation of sequence N<sup>o</sup> 5 and its lowest energy conformation found during the simulation. The RMS distance is 2.1 Å.

Fig. 7. MD folding trajectories of sequence N<sup>o</sup> 5. Each trajectory (small connected full dots) is obtained starting from 10 conformations (large full dots) randomly chosen from the 100 conformations of lowest energy obtained by contact map dynamics. Conformations (empty large dots) collected during other MD simulation (see §II C) are also shown as reference.