Protein folding simulations with genetic algorithms and a detailed molecular description¹

Abstract

We have explored the application of genetic algorithms (GA) to the determination of protein structure from sequence, using a full atom representation. A free energy function with point charge electrostatics and an area based solvation model is used. The method is found to be superior to previously investigated Monte Carlo algorithms. For selected fragments, up to 14 residues long, the lowest free energy structures produced by the GA are similar in conformation to the corresponding experimental structures in most cases. There are three main conclusions from these results. First, the genetic algorithm is an effective method for searching amongst the compact conformations of a polypeptide chain. Second, the free energy function is generally able to select native-like conformations. However, some deficiencies are identified, and further development is proposed. Third, the selection of native-like conformations for some protein fragments establishes that in these cases the conformation observed in the full protein structure is largely context independent. The implications for the nature of protein folding pathways are discussed.

Introduction

In order to determine protein structure from amino acid sequence, two central problems must be overcome: some form of free energy function must be developed that is able to distinguish between the functional conformation and all the rest, and a conformational search method must be devized that can find that functional conformation using available computing resources.

The most realistic method to date for studying protein motion is Cartesian space molecular dynamics, using an explicit description of all the protein and solvent molecule atoms Allen and Tildesley 1989, Brooks et al 1991. Force fields at this level are highly evolved Aqvist et al 1985, Weiner et al 1984, Brooks et al 1983, Hagler 1985, Lee et al 1995, Kang et al 1996, and a number of studies have shown that there are significant free energy minima near to the functional conformation Brunne et al 1995, Kitson et al 1993. It is not known whether, given enough computer time, simulations started from arbitrary points would converge to the functional conformation. The very large computing requirements have so far prevented simulation of the behavior of proteins or even peptides for more than a few nanoseconds, whereas a period of microseconds to milliseconds would be needed to reproduce observed in vitro folding behavior of peptides (Williams et al., 1996) and small proteins (e.g. Kuszewski et al., 1994). The principal limitation on the simulation scale is that the energy surface with a full atom description is very fine grained. That is, the energy changes rapidly as a function of atomic position. This problem becomes more acute as a molecular system becomes more compact.

To overcome these difficulties, two parallel strategies have been pursued. First, the description of the system has been simplified, to smooth the energy surface. Simplification ranges from the use of implicit solvation models (Lazaridis et al., 1995), slightly (Brooks & Head-Gordon, 1991), moderately (Srinivasan & Rose, 1995) or drastically (Sun et al., 1995) reduced side-chain descriptions, approximate main-chain descriptions (Brooks & Head-Gordon, 1991), and lattice models (Kolinski & Skolnick, 1994), extending to one lattice point per two residues (Park & Levitt, 1996). It is clear that there must be a price to pay for these simplifications, but so far there have been few studies Mounge et al 1995, DeBolt and Skolnick 1996, Huang et al 1996 to establish to what degree the different simplifications affect the ability to identify the functional conformation reliably. We have taken a relatively conservative approach. An implicit solvent description is used, together with all oxygen, nitrogen, carbon and sulfur atoms, as well as polar hydrogen atoms. An empirical force field, including a Coulomb law description of electrostatics with partial atomic charges is used. Advantage is taken of potential of mean force methods (Avbelj & Moult, 1995b) to parameterize this and other terms against experimental protein structures. The ability of this force field to identify conformations close to the functional one has been investigated (Braxenthaler et al., 1996).

The second strategy to over come the search problem is to take larger steps in the conformational space. Covalent geometry (i.e. bond lengths and bond angles, and planarity of conjugated atoms sets) is approximately constant for subgroups of protein atoms, so that conformational freedom may be expressed in terms of the values of the dihedral angles around single bonds. Although molecular dynamics in this space allows a larger time step, greater computational cost partly offsets the gains (Rice & Brünger, 1994). The restrictive move size and the cost of calculating derivatives in traditional molecular dynamics is overcome by using a Monte Carlo (MC) or genetic algorithm (GA) procedure. Much larger changes in conformation can be obtained through the MC procedure. Large changes in one or more dihedral angles are introduced, and the resulting conformation accepted if the evaluated free energy decreases, or if the increase is not too large. A number of MC studies of peptide folding have been made Abagyan and Totrov 1994, Avbelj and Moult 1995a. The disadvantage is that, because the space of a protein molecule is densely packed with atoms, many moves will be rejected as energetically unacceptable.

We have previously investigated the effectiveness of a torsion space MC procedure at finding the experimental conformation of fragments of protein molecules (Avbelj & Moult, 1995a). The method was partly successful, in that in a number of cases low free energy, native-like conformations were generated. However, competing lower free energy but less native-like conformations were often encountered, and the lowest free energy structures were not as low as the corresponding minimized experimental structures. Three possible explanations for these limitations are: (1) Since these are fragments of proteins, the preferred conformation may not necessarily be that found in the full context. (2) The free energy function may have a large number of false minima. (3) The search may not have converged to the lowest free energy conformation possible. The fact that the experimental structures consistently had slightly lower values supported the hypothesis that this last factor played a significant role. We therefore sought a more effective way of searching in torsion space in the compact states of a protein molecule.

Genetic Algorithms (GA) offer one way of searching more effectively in crowded spaces. They were first introduced by Holland (1975) to simulate the processes of natural selection at the genetic level, but rapidly gained acceptance as general search methods (Goldberg, 1989). Possible solutions to a search problem are represented by genes, and a fitness function is used to evaluate the fitness of each gene. A population of genes evolves by three mechanisms. (1) Point mutations, in which the value at a specific place on a single gene is changed. (2) Crossovers, in which a portion of one gene replaces the equivalent portion of another one. (3) Survival of the fittest genes from one generation to the next. Tests on two and three-dimensional lattice models of proteins Unger and Moult 1993a, Unger and Moult 1993b have demonstrated that GAs are more effective than MC at finding the global energy minimum. It is still less clear how well the method performs on more detailed descriptions of a protein chain.

Dandekar and Argos 1992, Dandekar and Argos 1994 used genetic algorithms to fold C^α backbone models of proteins in real space off lattice. A simplified bit-string encoding of φ/ψ space, allowing each φ/ψ pair to adapt one of seven possible angle combinations (Rooman et al., 1991) was used. Each side-chain was represented by a sphere of 1.9 Å and the global energy or fitness function consisted of hydrogen bonding, secondary structure preference and a hydrophobic scatter term, each of which were scaled according to heuristic constants. The method generally relied on the correct pre-assignment of secondary structure for success.

Sun et al. (1995) used a method similar to that of Dandekar and Argos 1992, Dandekar and Argos 1994 in that all secondary structure was explicitly introduced. A full backbone representation was used, with one virtual atom per side-chain. A simulation consisted of mutations of five degree torsion changes and crossovers at non-secondary structure positions. This procedure was able to find conformations which had a low root-mean-square (RMS) deviation to the experimental structure in five of the seven test cases.

Sun (1993) has presented the most detailed GA for protein folding to date. A full atom polypeptide backbone together with one atom at the centroid of each side-chain was used. The GA was driven by a dictionary of di,, tri-, tetra-, and penta-peptides. A peptide segment was chosen from the database of peptides by sequence. The search method was applied to two small proteins melletin and pancreatic polypeptide. The conformation of a small cyclical peptide was also reproduced. All the structures were generated with an RMS deviation of 1.7 Å within the final population in the GA simulation. These results were impressive, although the proteins simulated were present in the dictionary which was used to derive the GA. The search method is said to be approximately 100 to 200 times more efficient than MC simulated annealing protocols.

A more detailed review of the use of GAs for protein structure prediction may be found in Pedersen & Moult (1996).

The search method, presented here, is tested on a set of 28 fragments of protein structures, up to 14 residues long. The fragments were selected on the basis of experimental data and energetic criteria indicating a preference to adopt a native-like structure independent of the presence of the rest of the protein. They were identified from a survey of the current literature on the folding of proteins and protein fragments. These fragments have the advantage of being relatively small, and therefore present a more tractable search problem than complete proteins. A second advantage is that generation of a native-like conformation in an objective search provides added information about the context independence of a fragment, and its possible role in early folding.

The structure of the rest of the paper is as follows. In the first section, the force field is briefly introduced. Next, details of the GA implementation are described, together with a study of optimization of the GA parameters. Then the results of determining the conformation of a standard set of protein fragments are presented, followed by an analysis of some of the factors contributing to problem cases. Finally, we discuss the significance of the results for the further development of search methods and force fields, and explore the implications of the results for the mechanisms of protein folding.