1472-6807-2-4 1472-6807 Research article <p>A protein folding potential that places the native states of a large number of proteins near a local minimum</p> Chhajer Mukesh mchhajer@email.unc.edu Crippen M Gordon gcrippen@umich.edu

Department of Chemistry, University of North Carolina, Chapel Hill, NC 27599, U.S.A

College of Pharmacy, University of Michigan, Ann Arbor, MI 48109-1065, U.S.A

BMC Structural Biology 1472-6807 2002 2 1 4 http://www.biomedcentral.com/1472-6807/2/4 10.1186/1472-6807-2-4 12165098
5 6 2002 6 8 2002 6 8 2002 2002 Chhajer and Crippen; licensee BioMed Central Ltd. This is an Open Access article: verbatim copying and redistribution of this article are permitted in all media for any purpose, provided this notice is preserved along with the article's original URL. decoys local minimum pairwise interaction energy function solvation energy function

Abstract

Background

We present a simple method to train a potential function for the protein folding problem which, even though trained using a small number of proteins, is able to place a significantly large number of native conformations near a local minimum. The training relies on generating decoys by energy minimization of the native conformations using the current potential and using a physically meaningful objective function (derivative of energy with respect to torsion angles at the native conformation) during the quadratic programming to place the native conformation near a local minimum.

Results

We also compare the performance of three different types of energy functions and find that while the pairwise energy function is trainable, a solvation energy function by itself is untrainable if decoys are generated by minimizing the current potential starting at the native conformation. The best results are obtained when a pairwise interaction energy function is used with solvation energy function.

Conclusions

We are able to train a potential function using six proteins which places a total of 42 native conformations within ~4 Å rmsd and 71 native conformations within ~6 Å rmsd of a local minimum out of a total of 91 proteins. Furthermore, the threading test using the same 91 proteins ranks 89 native conformations to be first and the other two as second.

Background

For the development of an all-encompassing potential energy function for the protein folding problem that can simulate both the thermodynamic and kinetic processes, it is argued that one must start from first principles 1. However, given the complexity of the problem, involving thousands of atoms and an extremely large number of conformations available to these atoms, the available computational power posses a serious restriction. Alternatively, one could use the available experimental results to develop empirical potential functions using a simplified representation for the system which does not require an explicit enumeration of the entire system. Here, we present one such simplified method to train a potential energy function using a small set of proteins which places a significant number of additional native conformations near its minima.

Ideally, we would like to obtain a potential energy function which assigns correct free energy to the native and non-native states for all the proteins. However, this requires knowing the complete distribution of energy states and total number of degeneracies for both the native and denatured states, quantities which are not readily available even for a simplified representation of our system. This has led to design of potential functions with restricted objectives.

Knowledge-based potentials are obtained from a survey of protein crystal structures. In this case, there is no training of parameters. The residue-residue interaction can be defined as contact vs. no-contact type234 or, alternatively, can also be defined as distance dependent567. The resulting potentials depend upon the training set of proteins, reference state, and the functional form8. Thomas and Dill910 have pointed out that these potential functions do not assign the correct value to the energy parameters, though they do rank them correctly.

On the other hand, the potential parameters are explicitly trained in the optimization-based methods. In these cases, the goal is to obtain a potential function which either maximizes the difference between the native and average non-native state energy (Z-score11), maximizes the ratio of folding to glass transition temperatures 12 or requires that all the non-natives be of energy higher than the native 13. Here, the native state is represented by a single conformation. While the first two approaches provide reasonable results, they do not guarantee complete success in terms of fold recognition even for the training set proteins. A somewhat better criterion is to require that all the non-native conformations have energy higher than the native. In this case, if the energy function is chosen such that the total energy of the conformation is a linear function of its adjustable parameters, then the problem can be stated as a linear programming problem with a suitable objective function1314151617181920.

There are five factors that govern the training of a potential function using optimization-based methods: (i) representation of amino acids, (ii) choice of the interaction energy function, (iii) training set of proteins, (iv) generation of the alternative conformations or decoys and (v) the objective function.

The representation of the amino acids depends upon the amount of atomic-level detail one wishes to incorporate in the calculations, starting from a single point per residue representation91011 to an all-heavy atom representation2021. As more and more details are included, the complexity of the calculations increases; however, a coarse-grained model can not be expected to produce very refined structure predictions. While the lattice models have usually stayed with single point representation, the continuous state models have used different variations from a single point per residue to all-atom representation.

The choice of the interaction energy function to some extent depends upon the physical property of the amino acid which is deemed important, as well as on the type of amino acid representation used. The most commonly used pairwise inter-residue interactions can be treated as contact/no contact 141522, discrete distance ranges 161723, and continuously varying functions of distance 1318192425. The energy function should be flexible enough without causing overfitting. While having more parameters does provide the flexibility, having too many of them is not always helpful, as seen by Park et. al. 25 where a 80,000 parameter potential performed worse than a simple contact type function. Such discontinuous functions could also lead to problems if used for kinetic calculations or for local optimization. Furthermore, it has been shown that distance dependent energy functions perform better than contact type in a continuous conformation space 24.

Furthermore, excluded-volume effects play an important role in the performance of a potential function. As Thomas & Dill 10 have pointed out, one of the reasons for failure of the potential functions is the absence of any excluded-volume effects in the interaction energy models. Similarly, Park & Levitt 24 show that a van der Waals type energy function gives better results than contact type or solvation potential alone, and a combination of surface energy function with a van der Waals type energy function performs even better.

Another crucial requirement for the development of potential functions using optimization-based methods is the generation of realistic and challenging decoys. A great deal depends upon the set of decoys used to define the non-native state. It has been demonstrated by various investigators 141926 that using a set of decoys obtained mostly by threading does not enforce a very stringent training criterion, and the potential functions so obtained fail to even place the training set native conformations near a local minimum, let alone the global minimum. Similarly, one can obtain low energy conformations by manipulation of conformations in a contact map representation, though this may lead to geometrically unrealizable conformations 14. Other methods which have been used to generate low energy decoys include Monte-Carlo simulations 27 and molecular dynamics simulations 2829, both of which require significant computation time, discretization of conformational variables 24, and inclusion of rigid pieces of structures 30.

A necessary, though not sufficient, condition for a well behaved potential function is that the native conformation be at or near a minimum of the potential function. As pointed out earlier 141926, the potential functions derived with mostly threaded conformations fail to satisfy this condition, even when a large set of conformations is used. The reason for such a failure is that the training conformations are fixed and not allowed to relax given a potential function. This results in most of the conformations being of significantly higher energy and only a small set of conformations end up providing an active constraint set. However, if the inactive conformations are allowed to relax, they would provide a much more stringent set of constraints and will improve the performance of the potential function significantly. In this work, we place a milder condition on the training by only asking that the native be near a local minimum. No condition is placed for the global minimum. We use energy minimization to generate new decoys which are physically realizable low energy decoys and provide better training. Furthermore, a physically meaningful objective function is used during the quadratic programming. Here we train a potential function using only a few single domain, monomeric proteins which do not require any hetero groups or ligands to stabilize their folded state. These chains have less than 10% sequence alignment (using MOE 31) and significantly different crystal structures (average rmsd of ~11 Å). Once trained, the potential function is able to place 42 proteins, even including some multi-chain proteins and two CASP3 proteins, to within ~4 Å of the experimental native conformations. Here we present a method that trains a reasonably good potential function using only a few proteins.

Results and Discussion

As we have shown in an earlier work 26, the generation of good decoys is crucial for the training of the potential function. The potential functions trained mostly by the threaded conformations, even using a large set of training proteins, do not put the native conformations of even the training set proteins at a local minimum or near it, a necessary, though not a sufficient, condition for the stability of the native conformation. This is a serious shortcoming of the potential function, since without ensuring the local minimum condition, there is no hope of ensuring that the native is the global minimum conformation, much less that the native is thermodynamically stable. Here we use energy minimization starting from the native conformation to generate new decoys which depend upon the current potential function. However, since this process is much more time consuming than using a library of conformations, we use a small set of proteins to train our potential function, though the training is much more rigorous. We also use small proteins since computation time increases roughly as a square of the chain length. The proteins we use have very little sequence identity and have very dissimilar crystal structures. This allows for various residues to be in different environments. For the first set of three proteins (Set A: {layi, 1bk2, 1ubi}; see Table 1), we perform calculations using all the three energy functions to find the one which works best and then use that energy function to get a better potential function using a slightly expanded set of proteins. All the proteins we use in the training set are single chain, single domain proteins which do not require any hetero groups or any other ligands to form a stable folded state in the aqueous solution. The potential function so trained is then used on the test set to evaluate its performance. The test set contains some close homologs but mostly quite distinct proteins compared to the training set in terms of their size, sequence, and crystal structures.

<p>Table 1</p>

Effect of Solvation with three training proteins

without solvation (PIE3)

with solvation (PSE3)



protein

ρ

rmsd (Å)

ρ

rmsd (Å)


1ayi

0.165

2.03

0.150

1.87

1bk2

0.212

2.04

0.118

1.13

1ubi

0.128

1.47

0.168

1.94

1ame

0.386

3.99

0.606

6.04

1ayd

0.492

6.39

0.440

5.60

1bm8

0.651

8.17

0.390

4.79

1enh

0.363

3.55

0.324

3.27

1gzi

0.355

3.66

0.456

4.43

1igd

0.928

11.15

0.220

2.50

1mjc

0.750

8.45

0.313

3.29

1pgb

0.450

4.73

0.278

2.82

1ptf

0.820

9.85

0.450

5.07

1pwt

0.296

3.03

0.187

1.90

1szt

0.867

11.50

0.181

2.92

1vcc

0.625

7.67

0.423

4.91

There are two additional constants, G and ρcut (described in Parameter Adjustment section), which need to be fixed before the calculation can be performed. In our previous calculations 26, we had fixed G = 0.3 and ρcut = 0.10. The value of G defines the minimum energy separation between the native and an alternate conformation. We would like this value to be large so that the native is stable. The value of G depends upon the kind of energy function being used and in this work, we vary the value between 0.3 and 60.0. Similarly, the value of ρcut defines the basin of the native conformation and this value is varied between 0.05 and 0.25. To obtain the best potential function, i.e., a potential function that places the maximum number of native conformations of test set proteins close to its minima, we vary these two parameters and repeat the training process.

PIE function

The PIE function, shown in Figure 1, has 3 parameters per interaction type, hence a total of 900 parameters. The potential function was so trained that when the native conformations of the training set proteins are energy minimized, the resulting local minimum conformation is within 2.5 Å rmsd. For all choices of 0.3 ≤ G ≤ 60.0 and 0.05 ≤ ρcut ≤ 0.25, the potential was trainable, i.e., the potential function had enough flexibility that we did not run in to a situation where the QP could not come up with a solution, though at times solutions were suboptimal. The performance of these potential functions on the test set proteins varied greatly. In general high G and low ρcut values resulted in a very rough energy surfaces which led to large energy changes for small conformational variations. On the other hand, low G and high ρcut values resulted in a potential function that is too flat and allows easy conformational changes. For these proteins, we found that a value of G = 10.0 and ρcut = 0.12 gave us the best results. This potential function (PIE3) is able to place four of the 12 test protein native conformations to within ~4 Å rmsd of a local minimum (see Table 1). The other eight are between 4–12 Å rmsd.

<p>Figure 1</p>

Terms of pairwise interaction energy function, eq. (2).

Terms of pairwise interaction energy function, eq. (2). Heavy solid line is f0, thin solid line is f1, dashed is f2, and dotted is f3.

SE function

The SE function is very simple and has only 21 adjustable parameters, one for each residue and one for the overcrowding which prevents a complete collapse of the chain as well as prevents any unrealistic chain overlaps. Once again, we used the same three proteins and used various combinations of G and ρcut to obtain solutions. However, this time, with the SE function by itself, we failed to obtain a potential function even for a single protein, let alone all three proteins together. As soon as the conformations were generated by energy minimization process, we obtained a set of inequalities which resulted in a null feasible solution region for QP. This shows that the SE function by itself does not have enough flexibility to give a trained potential function when decoys are generated by more rigorous methods. Yet, this function was successfully trained by threading using experimental crystal structures 18.

PSE function

As with the PIE, the same three proteins are used for the training purposes. Since this is a sum of the other two energy functions, there are 920 parameters, the one parameter for overcrowding in the SE function having been removed since the f0 part of the PIE function performs a similar function. Once again, we obtain various sets of trained parameters depending upon the value of G and ρcut, and once again we find that G = 10.0 and ρcut = 0.12 give the best results for the test proteins. In this case (PSE3), six of these proteins remain within ~4 Å rmsd and all 12 are within ~6 Å. Considering that the test set contains both close homologs and significantly different proteins from the training proteins in terms of their size, sequence, and structure, this is a very encouraging result since all 12 test proteins remain within ~6 Å compared to ~12 Å for PIE. It shows that using a PSE function, one can train a potential function using only a few proteins which would be applicable to many other proteins. The comparison between the results of the PIE and PSE3 functions is shown in Table 1.

The training for PSE3 was performed using three proteins (Set A = {1ayi, 1bk2, 1ubi}), none of which contains Cys. However, there were many other interactions which remained untrained in the final version of PSE3. To get all 300 pairwise interactions to train at least moderately, we had to either select larger proteins which contain all residues in sufficient number or enlarge the training set even if all the residues are not present in any one of the protein. Since computation time increases as square of the chain length, we chose to enlarge the training set using smaller proteins. A set of four proteins (Set B: {1enh, 1bk2, 2era, 1ubi}) contains all the 20 amino acids collectively. However, in the final trained potential function (PSE4), eight of the interactions still remain completely untrained. This set was then further expanded by including two additional proteins (Set C: {1enh, 1bk2, 2era, 1ubi, 1ail, 1dsl}). The final trained potential function using Set C (PSE6) has all the pairwise interactions at least mildly trained. Set C contains two proteins (1enh and 1ail) with mostly α-helical structure, three (1bk2, 2era and 1dsl) with mostly β-sheet structure and one (1ubi) with α/β structure. However, there is very little sequence identity between any of these proteins, and their crystal structures are also very different (less than 10% sequence alignment and rmsd of 8–14 Å in their crystal structures using MOE). The corresponding solvation parameters for various residues are shown in Table 2.

<p>Table 2</p>

Solvation parameters

Residue

PSE4

PSE6


Gly

0.39958884300

-0.17727756470

Ala

0.86811242830

-0.36213127280

Val

0.62981201440

-0.36741883040

Leu

0.36371810330

0.82156988920

He

0.89053343600

0.09748013457

Cys

1.39536632100

0.09762673560

Met

-0.42954544170

0.15720350550

Phe

-0.01532194602

-0.22049718860

Pro

1.60142719700

0.03948009827

Tyr

-0.37437233800

0.71270611040

His

0.45098491050

-0.44487856400

Trp

1.01909618300

0.26052565020

Ser

-0.19148377760

-0.44395269350

Thr

0.47362962090

0.11467398110

Lys

-0.14119705180

-0.15725262960

Arg

-0.30176394520

-0.30074012620

Asp

0.42807283590

-0.58936068250

Asn

0.30127878790

0.52122075820

Glu

-0.19447595190

0.69933440120

Gln

-0.58651775250

0.01869831326

The best results for Set B (PSE4) are obtained using G = 9.9 and ρcut = 0.12. The potential function PSE4 places 14 of the 15 test set native conformations to within ~4 Å of a local minimum (see Table 3) which is a significant improvement over PSE3. The PSE4 was further applied to another 72 proteins (see Table 4) for a total of 87 test proteins. Of these 87 test proteins, 31 native conformations were placed within ~4 Å rmsd of a local minimum (including three CASP3 proteins), and 59 are within ~6 Å rmsd. Thus, PSE4 places a total of 35 native conformations (including the training set proteins) of very different sequences and crystal structures to within 4 Å rmsd and 63 to within 6 Å rmsd of a local minimum out of a total of 91 proteins. The set of 91 proteins contains a variety of proteins, some of which are multi-chain (e.g., 1d3b) or require a large group for stabilization (e.g., 1cc5).

<p>Table 3</p>

Comparison of PSE4 and PSE6

protein

fold type

chain length

PSE4

PSE6



rmsd (Å)

ρ

rmsd(Å)

ρ


1ail

α

70

2.24

0.178

1bk2

β

57

1.65

0.173

2.45

0.251

1dsl

β

87

2.10

0.181

1enh

α

54

1.63

0.160

2.19

0.216

1ubi

α/β

76

1.03

0.090

2.41

0.211

2era

β

62

1.22

0.106

2.36

0.212

1a19.A

α/β

90

3.03

0.256

3.20

0.264

1ame

α/β

67

4.02

0.400

3.76

0.375

1ayd

α/β

101

3.53

0.288

3.44

0.276

1ayi

α

87

3.96

0.323

5.42

0.438

1bm8

α/β

99

2.90

0.237

4.97

0.409

1dsl

β

87

3.98

0.362

1gzi

α/β

65

3.21

0.327

3.31

0.334

1igd

α/β

61

2.66

0.234

3.33

0.290

1mjc

β

69

3.50

0.334

2.90

0.278

1ops

α/β

64

2.38

0.238

3.46

0.347

1pgb

α/β

56

3.01

0.296

2.88

0.275

1ptf

α/β

88

4.58

0.421

3.86

0.340

1pwt

β

61

2.17

0.211

3.39

0.323

1szt

α

68

2.99

0.187

7.91

0.533

1vcc

α/β

77

3.67

0.307

3.36

0.287

1yhb

β

87

7.44

0.523

<p>Table 4</p>

Comparison of PSE4 and PSE6 (additional tests)

protein

chain length

PSE4

PSE6



rmsd (Å)

ρ

rmsd (Å)

ρ


1acp

77

4.88

0.436

6.22

0.577

1ail

70

5.91

0.475

1amm

174

6.86

0.464

3.72

0.238

1bdo

79

3.29

0.289

3.30

0.288

1beo

98

4.83

0.381

3.77

0.304

1ble

161

4.85

0.338

4.16

0.280

1bor

56

6.95

0.750

5.42

0.619

1c94.A

37

2.48

0.149

5.10

0.318

1cc5

83

7.04

0.657

13.89

1.273

1cd8

114

4.18

0.321

4.88

0.363

1cfe

135

6.50

0.474

7.87

0.543

1chd

198

4.11

0.279

4.96

0.333

1cm

46

3.91

0.396

4.30

0.484

1tf

68

2.52

0.248

3.20

0.302

1cyo

88

3.83

0.320

5.94

0.492

1d3b.A

72

2.93

0.268

2.63

0.240

1d3b.B

81

3.65

0.275

3.60

0.261

1div

149

11.11

0.473

11.06

0.489

1ecd

136

5.38

0.375

5.02

0.367

1exg

110

3.61

0.274

4.03

0.303

1f3g

150

4.88

0.357

3.83

0.279

1fbr

93

9.01

0.604

6.61

0.428

1fdx

54

9.29

0.869

5.72

0.703

1fkf

107

5.59

0.450

3.56

0.279

1fxd

58

4.01

0.428

3.23

0.341

1hfh

120

9.37

0.527

4.12

0.227

1hoe

74

3.12

0.285

6.47

0.560

1ife

131

4.53

0.341

3.91

0.278

1jpc

108

5.55

0.399

6.01

0.427

1knb

186

5.18

0.320

4.64

0.291

1kuh

132

8.71

0.653

6.68

0.477

1lba

146

8.15

0.598

3.90

0.288

1lzl

130

5.18

0.388

4.33

0.324

1mai

119

4.99

0.380

3.73

0.280

1paz

120

3.87

0.305

3.94

0.317

1pcy

99

6.70

0.537

5.29

0.447

1pdo

129

4.66

0.342

10.08

0.691

1pkp

145

4.19

0.287

7.10

0.486

1poa

118

8.65

0.681

6.70

0.517

1poc

134

13.60

0.982

8.66

0.587

1r69

63

5.19

0.517

3.21

0.311

1ra9

159

7.24

0.497

3.34

0.220

1rie

127

5.17

0.393

4.98

0.364

1rsy

135

6.47

0.439

4.98

0.328

1skz

104

6.60

0.458

7.04

0.512

1snc

135

3.76

0.275

6.96

0.492

1vhh

157

5.12

0.361

6.81

0.465

1whi

122

3.49

0.264

3.75

0.280

1xnb

185

3.22

0.213

3.47

0.229

1ycc

108

6.50

0.515

5.35

0.429

1yua

122

8.53

0.626

7.49

0.489

2abd

86

8.38

0.752

5.21

0.432

2 end

137

8.20

0.534

5.71

0.356

2fx2

147

4.18

0.309

4.87

0.350

2fxb

81

4.08

0.382

3.40

0.322

2hbg

147

5.72

0.400

3.99

0.278

21hb

29

9.24

0.703

3.42

0.299

2mcm

112

3.77

0.303

4.87

0.360

2rhe

114

3.47

0.263

4.25

0.315

2m2

155

7.15

0.507

8.33

0.542

2sns

141

6.61

0.459

7.74

0.522

2stv

184

6.31

0.366

6.91

0.421

2tgi

112

7.88

0.463

4.05

0.238

3b5c

85

5.62

0.471

5.23

0.443

3chy

128

3.30

0.261

3.46

0.268

451c

82

5.65

0.487

5.74

0.514

4fd1

106

6.15

0.510

3.27

0.286

4fxn

138

4.78

0.357

5.43

0.415

4icb

76

4.67

0.431

3.85

0.340

5rxn

54

4.93

0.509

4.10

0.421

7rsa

124

8.70

0.605

5.74

0.404

The best result for Set C (potential function PSE6 with G = 13.1 and ρcut = 0.15) places 11 out of 14 test set proteins to within ~4 Å rmsd (see Table 3). In all, out of 91 proteins including the training set, 71 are within ~6 Å and 42 are within ~4 Å (including two CASP3 proteins). The fact that PSE6 is able to place 71 out of 91 native conformations within ~6 Å rmsd of a local minimum for proteins of very different sequences, crystal structures, and length points to the robustness of the potential function, and, therefore, the usefulness of the approach.

Threading test

To further test the validity and usefulness of the final potential function PSE6, we performed an ungapped threading test using all the 91 proteins. We first energy minimized each of the native conformations and obtained the conformation of the local minimum using potential PSE6. These local minimum conformations were then used as decoys during ungapped threading tests for each protein. The threaded conformations were not energy minimized. In the threading test, instead of the usual single point representation for each residue, we use a five point representation for each residue, since the potential PSE6 has been trained using this representation. Out of 91 proteins, 89 native conformations are ranked first and the other two (1dsl, 1snc) are ranked second. This test is conducted using a total of 180,878 decoys. Therefore, we see that a potential function trained using the procedure outlined in this work also performs well when subjected to the threading test, whereas the potentials designed using only threaded conformations do not perform well when decoys are generated using energy minimization.

One of the drawbacks of this procedure is that the training process is not sequential and, therefore, quite time consuming, i.e., addition of inequalities for a new protein after the training for one protein is complete does not necessarily give better results. Any change in protein set, G, ρcut or any other strategy change requires starting all over again, since the set of conformations generated for the training are based on the current set of parameters, i.e., the process is memory dependent. However, given that the potential energy surface is of very high dimension and highly irregular, it is to be expected.

By only considering the energy minimization of the experimental native conformation, we reduce the amount of time required to train the potential function. Each potential function could be trained within a week running on four Sun workstations (CPU speeds: 135 MHz-450 MHz). However, we had to perform many training runs to get the best possible potential function for each set, and at this point it guarantees only the local minimum condition for the native state, not the global one.

Conclusions

A necessary, though not sufficient, condition for the stability of the native state of a protein is that the experimental native conformation be at or near a minimum of the potential function. We show that by using only a small set of proteins, potential functions can be trained that put a large number of native states near a local minimum of the potential function. While this does not guarantee the stability of the native state, this does ensure that the native conformation is not sitting at highly unstable points on the potential surface, a condition which is encountered in other potential functions. Our best potential function (PSE6) obtained using training Set C (containing six proteins) is able to place a total of 42 native conformations to within ~4 Å rmsd and 71 native conformations to within ~6 Å rmsd of a local minimum of the potential function out of the 91 proteins. We also find that while the pairwise interaction energy function is trainable using the more rigorously generated decoys, the solvation potential alone is not. The solvation potential has only 21 adjustable parameters and can not provide enough flexibility to satisfy all the constraints in our case, though it was trainable using threaded conformations from PDB. The best results are obtained using the pairwise energy function in combination with the solvation energy function. We further test the final potential function PSE6 with threaded conformations for 91 proteins and find that 89 of the native conformations are ranked first and the remaining two are ranked second. The training process is further improved by the use of a physically meaningful quantity for the objective function in the QP optimization which helps in identifying a better solution.

Methods

Generation of native structure

For the purposes of training, we select only single domain, monomeric proteins which do not require any hetero atoms or ligands for their stabilization in the folded state. Furthermore, the atom positions are available for the main chain heavy atoms and at least the corresponding Cβ atom of the residue. The crystal structures (obtained from PDB32) for the protein chains are fitted to a standard geometry continuous state model (very similar to the one used by Dill et. al.33). In the fitted model, each side chain is represented by a single interacting site, located at the Cβ, while keeping all the main chain heavy atoms, i.e., each residue is represented by five interacting sites (united atom types). Standard values for bond lengths and bond angles are used, and all peptide bonds are kept in the trans conformation. Thus, only torsion angles (φ, ψ) are allowed to vary between -180° and 180°. The details of the fitting procedure are described in ref. 20 and will not be repeated here. The fitted model is within 0.5 Å rmsd 34 of the PDB structure and is used as the native structure in our calculations. Thus, the native state is a single energy, non-degenerate state in our calculations.

Generation of the non-native decoys

The conformations of the non-native decoys are generated using the following procedures, (i) Parameter-Independent Decoys: Starting from the native conformation, all pairs of torsional angles are perturbed by -30° ≤ δ ≤ 30° at a time. For example, for a chain of 50 residues, there are 100 torsion angles and, therefore, 4950 pairs of torsional angles. This would give us 4950 different decoys, but see the Solution Procedure section for the precise protocol used. These conformations do not depend upon the current potential function, (ii) Parameter-Dependent Decoys: The native conformation is energy minimized with respect to torsion angles using the current potential function.

Energy functions

Since we generate conformations by energy minimization with respect to the torsion angles which are allowed to change continuously, we also need an energy function that is a continuous and differentiable function of angles. We use three different types of energy functions: the pairwise interaction energy (PIE) function, the solvation energy (SE) function, and the pairwise interaction with solvation energy (PSE) function. In our case, all these functions are continuous functions of separation between atoms and, therefore, the conformations can be energy minimized with respect to torsion angles to generate decoys once we use the standard geometric representation for the conformations.

Pairwise interaction energy (PIE) function

The pairwise interaction between a pair of atoms is represented by 19

where pk,titj are adjustable parameters, ti is the type of united atom i, and xij is the distance between united atoms i and j. The function fk(xij) is defined over a limited distance interval as follows:

where b = 4.0, a0 = 0.0, a1 = 4.0, a2 = 6.0 and a3 = 8.0. The first term in each interaction, f0(xij), provides some steric repulsion for atoms sitting too close to each other and its coefficient is always fixed at 10 which prevents a complete collapse of the chain. The other three functions have adjustable coefficients which are allowed to vary between -20 and 20. The functional form of fk(xij)'s is shown in Figure 1. Each interaction energy is represented by four parameters. The maximum interaction distance is 12 Å. Equation (1) has following properties:

etitj(0) = 10

etitj(xij) = 0 for xij ≥ 12 Å

The interaction energy function etitj (xij) is continuous and differentiable for all xij ≥ 0, and is a linear function of adjustable parameters.

Furthermore, we consider 24 atom types: 20 side chain atom types representing 20 different amino acids and four main chain heavy atoms (Cα, C, N, and 0). This gives us a total of 24 × 25/2 = 300 different types of interactions and a total of 1200 parameters. In the present study, 300 of these parameters have a fixed value of 10, and 900 are adjustable. The total pairwise energy is

where the summation is over all pairs of united atom types separated by at least three covalent bonds. Since the total pairwise energy of any conformation is just the sum of individual pairwise interaction energies, the total pairwise energy is also a linear function of the adjustable parameters.

Solvation energy (SE) function

The solvation energy represents the change in energy due to the burial of various residues. This is an average contribution for each residue based upon how many other residues are surrounding it and given by

where pk(i) is the energy contribution for a residue of type k at sequence position i along the chain surrounded by other residues, and Ps is the penalty for overcrowding. Here, the first summation is over all residue pairs separated by at least two residues along the chain backbone, and second summation is over all those residues for which the total number of contacts exceeds the maximum allowable number. The number of excess contacts Si are

The function Ctitj (xij) determines the burial of each residue where xij is the separation between two side chain Cβ atoms situated at positions i and j. These are sigmoidal functions (see Figure 2). Both Ctitj (xij) and Ck(i), max were obtained by Dombkowski and Crippen18 from the survey of PDB crystal structures, and these functions are predetermined in our training procedure. Only the pk(i)'s and PS are determined here. This energy function has therefore only 21 adjustable parameters.

<p>Figure 2</p>

Functional form of Ctitj (xij), eq. (4).

Functional form of Ctitj (xij), eq. (4).

Pairwise interaction with solvation energy (PSE) function

The PSE function includes both the pairwise interaction energy and the solvation energy without the overcrowding part since the f0(xij) part of the PIE function ensures that the chain does not collapse to a point. The total energy in this case is given by

and it has a total of 920 adjustable parameters.

Parameter adjustment using quadratic programming

Given the above energy functions, the total energy of any conformation is a linear function of the adjustable parameters, irrespective of the energy function being used. This property allows us to use quadratic programming (QP) to adjust the values of these parameters. We choose the following objective function for the optimization purposes:

where αn,k is the sum of the coefficients of the Np parameters pk when the derivative of the energy of a conformation is evaluated at the native conformation of the n-th protein. While the first term minimizes the gradient of the energy with respect to dihedral angles at the native conformation for all training set proteins, the second term keeps the values of parameters small to reduce the roughness of the energy surface. Having a continuous energy function is essential for using such a physically meaningful objective function. By minimizing the magnitude of the gradient of energy of the native conformation, we obtain a set of parameters such that the native conformation approximates a stationary point of the potential function. The total number of parameters NP depends upon the type of energy function is being used, namely 900 for the PIE function, 21 for the SE function, and 920 for the PSE function. The constraints are of the form

ΔE = Enon - Enat >g     (8)

where

where G and ρcut are adjustable constants. If ρ < ρcut, we consider the decoy to be conformationally identical to the native and no inequality is generated even if the energy of the decoy is lower than that of the native. The constants G and ρcut define the minimum energy separation between the native and decoys, and the radius of the basin of the native, respectively. They may be changed from one training run to another; however, not during the same training run. Here, ρ is given by 35

where Rgj is the radius of gyration of the jth conformation and D1,2 is the customary root mean square deviation in Cα coordinates after optimal superposition of the two conformations, taken to be the native and an alternate. By allowing g to vary linearly with ρ near the native structure, we avoid sudden jumps in the potential energy surface and keep the variation in the energy function smooth near the native. Finally, we need to provide bounds on our parameters, and we have chosen pk = 10 for all pk's which multiply f0(xij)'s and -20 ≥ pk ≤ 20 for the rest of them. These values are chosen so that we are able to find a feasible solution while keeping most of the parameters away from the extreme values. Due to the choice of , we use quadratic programming (QP) instead of linear programming to solve the set of inequalities. 36

Solution procedure

We obtain the initial set of parameters by minimizing using QP. The decoys are generated using the methods described earlier. First, the decoys are generated by perturbing pairs of torsion angles by an amount δ. For a chain of n residues, there are 2n(2n-1)/2 such pairs. Even for a chain with 50 residues, this would produce 4950 conformations. Initially, since the parameter set is completely untrained, a large number of inequalities are generated for a training set of three to six proteins ranging in length from 50–90 residues. To keep the number of inequalities to a reasonable value, we first perturb pairs of angles separated by a certain number of covalent bonds n and its multiples. For example, in the first pass, we change the pair of angles separated by n = 9 covalent bonds or multiples of it, i.e., 9, 18, 27, etc. In this case, starting from the fifth torsion angle from the N-terminus, we perturb pairs of angles located at (5+9k, 5+9l) positions along the chain where kl = 0, 1, 2, etc, collect all the inequalities and run them through the QP to obtain a new set of parameters. Next, we take pairs of angles separated by 4 backbone bonds and add new inequalities to the previous set and obtain a new set of parameters by QP. This process is repeated for pairs of angles separated by two, one and zero backbone bonds. Next, we repeat this process for all pairs of angles by perturbing them by δ ± 2°, δ ± 4° and δ ± 6° where 12° ≤ |δ| ≤ 20°. This is the parameter independent generation of decoys.

The next step is parameter dependent generation of decoys. Here, we use energy minimization of native conformation (using BFGS method 37) to generate new decoys. During each cycle, we begin from the native conformation and energy minimize the conformations with respect to the torsion angles. After every 100 minimization steps, the conformational distance of the current conformation from the native is calculated in terms of ρ. If the ρ value is greater than the cutoff value ρcut, we say that the current conformation is sufficiently different from the native one, i.e., it does not belong to the basin of native, and since the energy of this conformation is lower than the native, a new constraint is generated. After complete minimization of the native conformations, if any new inequalities are generated, they are added to the set of inequalities and a new set of parameters is obtained by QP. If no new inequalities are generated, then the training is complete and a trained potential has been obtained. To keep the number of inequalities within a manageable limit, we periodically discard inequalities corresponding to nonnatives having energy significantly higher than the native state, i.e., having high slack value greater than S for the current parameter set. Initially, S is set at 500 and slowly reduced to 200.

Authors' Contributions

M.C. designed the fitting protocols, wrote the computer programs, chose training and test data, and carried out the calculations. G.M.C. contributed overall planning, advice, and discussions throughout. Both have read and approved the manuscript.

Acknowledgements

This work was supported by a grant from NIH (Grant #1-R01-GM59097-01). The authors would like to thank Dr. Y. Z. Ohkubo for allowing us to use some of his computer code. We would also like to thank the PDB contributors.

<p>Effective energy functions for protein structure prediction.</p> Lazaridis T Karplus M Curr Opin Struct Biol 2000 10 139 145 10.1016/S0959-440X(00)00063-4 10753811 <p>Estimation of effective contact energies from protein crystal structures: Quasi-chemical approximation.</p> Miyazawa S Jernigan R Macromolecules 1985 18 534 552 <p>Residue-residue potentials with a favorable contact pair term and an unfavorable high packing density term, for simulation and threading.</p> Miyazawa S Jemigan R J Mol Biol 1996 256 623 644 10.1006/jmbi.1996.0114 8604144 <p>Exploring conformational space with a simple lattice model for protein structure.</p> Hinds D Levitt M J Mol Biol 1994 243 668 682 7966290 <p>Derivation and testing of pair potentials for protein folding. When is the quasichemical approximation correct?</p> Skolnick J Jaroszewski L Kolinski A Godzik A Protein Sci 1997 6 676 688 9070450 <p>Boltzmann's principle, knowledge-based mean fields and protein folding. An approach to the computational determination of protein structures.</p> Sippl M J Comput-Aided Mol Design 1993 7 473 501 <p>Factors influencing the ability of knowledge-based potentials to identify native sequence-structure matches.</p> Kocher J-P Rooman M Wodak S J Mol Biol 1994 235 1598 1613 10.1006/jmbi.1994.1109 8107094 <p>Are proteins ideal mixtures of amino acids? Analysis of energy parameter sets.</p> Godzik A Kolnski A Skolnick J Protein Sci 1995 4 2107 2117 8535247 <p>An interactive method for extracting energy-like quantities from protein structure.</p> Thomas P Dill K Proc Natl Acad Sci USA 1993 93 11628 11633 10.1073/pnas.93.21.11628 <p>Statistical potentials extracted from protein structures: How accurate are they?</p> Thomas P Dill K J Mol Biol 1996 257 457 469 10.1006/jmbi.1996.0175 8609636 <p>How to derive a protein folding potential? A new approach to an old problem.</p> Mirny L Shakhnovich E J Mol Biol 1996 264 1164 1179 10.1006/jmbi.1996.0704 9000638 <p>Protein tertiary structure recognition using optimized Hamiltonians with local interactions.</p> Goldstein R Luthey-Shulten Z Wolynes P Proc Natl Acad Sci USA 1992 89 9029 50058 1409599 <p>Contact potential that recognizes the correct folding of globular proteins.</p> Maiorov V Crippen G J Mol Biol 1992 227 876 888 1404392 <p>Pairwise contact potentials are unsuitable for protein folding.</p> Vendruscolo M Domany E J Chem Phys 1998 109 11101 11108 10.1063/1.477748 <p>Can a pairwise contact potential stabilize native protein folds against decoys obtained by threading?</p> Vendruscolo M Najmanovich R Domany E Proteins 2000 38 134 148 10.1002/(SICI)1097-0134(20000201)38:2<134::AID-PROT3>3.0.CO;2-A 10656261 <p>Distance-dependent, pair potential for protein folding: Results from linear optimization.</p> Tobi D Elber R Proteins 2000 41 40 46 10.1002/1097-0134(20001001)41:1<40::AID-PROT70>3.3.CO;2-L 10944392 <p>On the design and analysis of protein folding potentials.</p> Tobi D Shafran G Linial N Elber R Proteins 2000 40 71 85 10.1002/(SICI)1097-0134(20000701)40:1<71::AID-PROT90>3.0.CO;2-3 10813832 <p>Disulfide recognition in an optimized threading potential.</p> Dombkowski A Crippen G Protein Eng. 2000 13 679 689 10.1093/protein/13.10.679 11112506 <p>Potential energy function for continuous state model of globular proteins.</p> Ohkubo Y Crippen G J Comput Biol 2000 7 363 379 10.1089/106652700750050835 11108468 <p>Constructing smooth potential functions for protein folding.</p> Crippen G J Mol Graph Mod 2001 19 87 93 10.1016/S1093-3263(00)00132-7 <p>An all-atom distance dependent conditional probability discriminatory function for protein structure prediction.</p> Samudrala R Moult J J Mol Biol 1998 275 895 916 10.1006/jmbi.1997.1479 9480776 <p>Folding Lennard-Jones proteins by a contact potential.</p> Clementi C Vendruscolo M Maritan A Domany E Proteins 1999 37 544 553 10.1002/(SICI)1097-0134(19991201)37:4<544::AID-PROT5>3.3.CO;2-Z 10651270 <p>Towards an energy function for the contact map representation of proteins.</p> Park K Vendruscolo M Domany E Proteins 2000 40 237 248 10.1002/(SICI)1097-0134(20000801)40:2<237::AID-PROT60>3.3.CO;2-G 10842339 <p>Energy functions that discriminate X-ray and near native folds from well-constructed decoys.</p> Park B Levitt M J Mol Biol 1996 258 367 392 10.1006/jmbi.1996.0256 8627632 <p>Factors affecting the ability of energy functions to discriminate correct from incorrect folds.</p> Park B Huang E Levitt M J Mol Biol 1997 266 831 846 10.1006/jmbi.1996.0809 9102472 <p>A protein folding potential that has a stable native state and correct free energy of unfolding,</p> Chhajer M Crippen G 2002 <p>Learning effective amino acid interactions through interactive stochastic techniques.</p> Micheletti C Seno F Banavar J Maritan A Proteins 2001 42 422 431 10.1002/1097-0134(20010215)42:3<422::AID-PROT120>3.0.CO;2-2 11151013 <p>Discriminating compact nonnative structures from the native structure of globular proteins.</p> Wang Y Zhang H Li W Scott R Proc Natl Acad Sci USA 1995 92 709 13 42689 7846040 <p>Using a hydrophobic contact potential to evaluate native and near-native folds generated by molecular dynamics simulations.</p> Huang E Subbiah S Tsai J Levitt M J Mol Biol 1996 257 716 725 10.1006/jmbi.1996.0196 8648635 <p>Statistical mechanics of protein folding by exhaustive enumeration.</p> Crippen G Ohkubo Y Proteins 1998 32 425 437 10.1002/(SICI)1097-0134(19980901)32:4<425::AID-PROT3>3.3.CO;2-Z 9726414 <p>MOE, Version 1999.09.</p> Chemical Computing Group Inc 1999 http://www.chem-comp.com <p>The Protein Data Bank.</p> Berman H Westbrook J Feng Z Gilliland G Bhat T Weissig H Shindyalov I Bourne P Nucleic Acids Res 2000 28 235 242 102472 10592235 10.1093/nar/28.1.235 <p>Protein structure and energy landscape dependence on sequence using a continuous energy function.</p> Dill K Phillips A Rosen J J Comput Biol 1997 4 227 239 9278057 <p>A discussion of the solution of the best rotation to relate two sets of vectors.</p> Kabsch W Acta Cryst 1978 A34 827 828 <p>Size-independent comparison of protein 3-dimensional structures.</p> Maiorov V Crippen G Proteins 1995 22 273 283 7479700 <p>CPLEX, Version 6.6. ILOG, Inc</p> 2000 http://www.ilog.com <p>Remark on "Algorithm 500: Minimization of unconstrained multi-variate function [E4]".</p> Shanno D Phua K ACM Trans Math Software 1980 6 618 622 10.1145/355921.355933