# Simulations of Protein Folding

## 1 PROTEINS

A protein is a linear chain of amino acids. The proteins of natural living organisms are composed of 20 different types of amino acids. A typical protein is a polymer of 300 amino acids, of which there are $`20^{300}=2\times 10^{390}`$ different possibilities. The human body uses about 80,000 different proteins for most of its functionality, including structure, communication, transport, and catalysis.

The order of the amino acids in the proteins of an organism is specified by the order of the base pairs in the deoxyribonucleic acid, DNA, of its genome. Human DNA consists of $`10^9`$ base pairs with a total length of 3m. Since three base pairs specify an amino acid, the code for the 80,000 human proteins requires only $`3\times 300\times 80,000=7\times 10^7`$ base pairs or 7% of the genome.

### 1.1 Amino Acids

The twenty amino acids differ only in their side chains. The key atom in an amino acid is a carbon atom called the $`\alpha `$-carbon, C<sub>α</sub>. Four atoms are attached to the C<sub>α</sub> by single covalent bonds: a hydrogen atom H, a carbonyl-carbon atom C, a nitrogen atom N, and the first atom of the side chain R of the amino acid. The carbonyl carbon C is connected to an oxygen atom by two covalent bonds and to a hydroxyl group OH by another covalent bond; the nitrogen atom N is attached to two hydrogen atoms, forming an amine group NH<sub>2</sub>. The backbone of an amino acid is the triplet N, C<sub>α</sub>, C.

Of the 20 amino acids found in biological systems, 19 are left handed. If one looks at the C<sub>α</sub> from the H, then the order of the structures CO, R, and N is clockwise (CORN). The one exception is glycine in which the entire side chain is a single hydrogen atom; glycine is not chiral.

### 1.2 Globular Proteins

There are three classes of proteins: fibrous, membrane, and globular. Fibrous proteins are the building materials of bodies; collagen is used in tendon and bone, $`\alpha `$-keratin in hair and skin. Membrane proteins sit in the membranes of cells through which they pass molecules and messages. Globular proteins catalyze chemical reactions; enzymes are globular proteins.

Under normal physiological conditions, saline water near pH=7 at 20–40 <sup>o</sup>C, proteins assume their native forms. Globular proteins fold into compact structures. The biological activity of a globular protein is largely determined by its unique shape, which in turn is determined by its primary structure, that is, by its sequence of amino acids.

### 1.3 Kinds of Amino Acids

The amino acids that occur in natural living organisms are of four kinds. Seven are nonpolar: alanine (ala), valine (val), phenylalanine (phe), proline (pro), methionine (met), isoleucine (ile), and leucine (leu).They avoid water and are said to be *hydrophobic*. Four are charged: aspartic acid (asp) and glutamic acid (glu) are negative, lysine (lys) and arginine (arg) are positive. Eight are polar: serine (ser), threonine (thr), tyrosine (tyr), histidine (his), cysteine (cys), asparagine (asn), glutamine (gln), and tryptophan (trp). The four charged amino acids and the eight polar amino acids seek water and are said to be *hydrophilic*. Glycine falls into a class of its own.

### 1.4 Protein Geometry

When two amino acids are joined to make a dipeptide, first the hydroxyl group OH attached to the carbonyl carbon C of the first amino acid combines with one of the two hydrogen atoms attached to the nitrogen N of the second amino acid to form a molecule of water H<sub>2</sub>O, and then a peptide bond forms between the carbonyl carbon C of the first amino acid and the nitrogen N of the second amino acid. A peptide bond is short, 1.33 Å, and resists rotations because it is partly a double bond.

To a good approximation, the six atoms C<sub>α1</sub>, C$`{}_{1}{}^{}{}_{}{}^{}`$, O<sub>1</sub>, N<sub>1</sub>, H<sub>1</sub>, and C<sub>α2</sub> lie in a plane, called the peptide plane. If a third amino acid is added to the carbonyl carbon C$`{}_{2}{}^{}{}_{}{}^{}`$ of the second amino acid, then the six atoms C$`{}_{\alpha 2}{}^{}\mathrm{}`$C<sub>α3</sub> also will lie in a (typically different) plane. Exceptionally, the peptide plane of proline is not quite flat because the side chain loops around, and its third carbon atom forms a bond with the nitrogen atom of the proline backbone.

### 1.5 The Protein Backbone

The protein backbone consists of the chain of triplets (N C<sub>α</sub> C)<sub>1</sub>, (N C<sub>α</sub> C)<sub>2</sub>, (N C<sub>α</sub> C)<sub>3</sub>, $`\mathrm{}`$, (N C<sub>α</sub> C)<sub>N</sub>. Apart from the first nitrogen N<sub>1</sub> and the last carbonyl carbon C$`{}_{n}{}^{}{}_{}{}^{}`$, this backbone (and its oxygen and amide hydrogen atoms) consists of a chain of peptide planes, C$`{}_{\alpha 1}{}^{}\mathrm{}`$C$`{}_{\alpha 2}{}^{}\mathrm{}`$ C$`{}_{\alpha n1}{}^{}\mathrm{}`$C<sub>αn</sub>. Since the angles among the four bonds of the C<sub>α</sub>’s are fixed, the shape of the backbone of peptide planes is determined by the angles of rotation about the single bonds that link each C<sub>α</sub> to the N that precedes it and the C that follows it. The angle about the N<sub>i</sub>-C<sub>αi</sub> bond is called $`\varphi _i`$, that about the C<sub>αi</sub>-C$`{}_{i}{}^{}{}_{}{}^{}`$ bond is $`\psi _i`$. The 2N angles $`(\varphi _1,\psi _1)\mathrm{}(\varphi _N,\psi _N)`$ determine the shape of the backbone of the protein. These angles are the main kinematic variables of a protein. The principal properties of proteins are discussed in the classic article by Jane Richardson .

## 2 PROTEIN FOLDING

The problem of protein folding is to predict the natural folded shape of a protein under physiological conditions from the DNA that defines its sequence of amino acids, which is its primary structure. This difficult problem has been approached by several techniques. Some scientists have applied all-atom molecular dynamics . We have used the Monte Carlo method in a manner inspired by the work of Ken Dill *et al.* .

Our Monte Carlo simulations are guided by a simple potential with three terms. The first term embodies the Pauli exclusion principle. Because the outer parts of atoms are electrons which are fermions, the Pauli exclusion principle requires that the side chains of a protein not overlap by more than a fraction of an angstrom. In our present simulations, we have represented each side chain as a sphere centered at the first carbon atom, the $`C_\beta `$, of the side chain or at the hydrogen atom that is the side chain in the case of glycine.

The second term represents the mutual attraction of nonpolar or hydrophobic amino acids. In effect the water electric dipoles, the free protons, the free hydroxyl radicals, and the other ions of the cellular fluid attract the charged and polar amino acids of a protein but leave unaffected the nonpolar amino acids. The resulting net inward force on the nonpolar amino acids drives them into a core which can be as densely packed as an ionic crystal.

The third term is a very phenomenological representation of the effects of steric repulsion and hydrogen bonding. For a given amino acid, this term is more negative when its pair of angles $`\varphi _i`$ and $`\psi _i`$ are in a zone that avoids steric clashes between the backbone and the side chain and that encourages the formation of hydrogen bonds between NH<sup>+</sup> and O<sup>-</sup> groups. One of these Ramachandran zones favors the formation of $`\alpha `$ helices, others favor $`\beta `$ structures.

We incorporate these zones in a Metropolis step with two scales, which we call zoning with memory. Each Monte-Carlo trial move begins with a random number that determines whether the angles $`\varphi _i`$ and $`\psi _i`$ of residue $`i`$ will change zone, *e.g.,* from its present zone to the $`\alpha `$ zone, the $`\beta `$ zone, or to the miscellaneous zone. If the zone is changed, then the angles $`\varphi _i`$ and $`\psi _i`$ revert to the values they possessed when residue $`i`$ was last in that zone. The trial move is then modified slightly and randomly.

### 2.1 Rotations

We have derived a simple formula for the 3$`\times `$3 real orthogonal matrix that represents a right-handed rotation by $`\theta =|\stackrel{}{\theta }|`$ radians about the axis $`\widehat{\theta }=\stackrel{}{\theta }/\theta `$:

$$e^{i\stackrel{}{\theta }\stackrel{}{J}}=\mathrm{cos}\theta Ii\widehat{\theta }\stackrel{}{J}\mathrm{sin}\theta +(1\mathrm{cos}\theta )\widehat{\theta }(\widehat{\theta })^𝖳$$

in which the generators $`(J_k)_{ij}=iϵ_{ikj}`$ satisfy $`[J_i,J_j]=iϵ_{ijk}J_k`$ and $`𝖳`$ means transpose. In terms of indices, this formula for $`R(\stackrel{}{\theta })=e^{i\stackrel{}{\theta }\stackrel{}{J}}`$ is

$$R(\stackrel{}{\theta })_{ij}=\delta _{ij}\mathrm{cos}\theta \mathrm{sin}\theta ϵ_{ijk}\widehat{\theta }_k+(1\mathrm{cos}\theta )\widehat{\theta }_i\widehat{\theta }_j.$$

In these formulae $`ϵ_{ijk}`$ is totally antisymmetric with $`ϵ_{123}=1`$, and sums over $`k`$ from 1 to 3 are understood.

### 2.2 Distance

A conventional measure of the quality of a theoretical fold is the mean root-mean-square distance $`d`$ between positions $`\stackrel{}{r}(i)`$ of the $`\alpha `$ carbons of the folded protein and those $`\stackrel{}{x}(i)`$ of the native structure of the protein,

$$d=\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{}}\left(\stackrel{}{r}(i)\stackrel{}{x}(i)\right)^2}.$$

The native states of many proteins are available from http://www.rcsb.org/pdb/. We have derived a formula for this distance in terms of the centers of mass $`\stackrel{}{r}=(1/n)_{j=1}^n\stackrel{}{r}(j)`$ and $`\stackrel{}{x}=(1/n)_{j=1}^n\stackrel{}{x}(j)`$, the relative coordinates $`\stackrel{}{q}(i)=\stackrel{}{r}(i)\stackrel{}{r}`$ and $`\stackrel{}{y}(i)=\stackrel{}{x}(i)\stackrel{}{x}`$, their inner products $`Q^2=_{i=1}^n\stackrel{}{q}(i)^2`$ and $`Y^2=_{i=1}^n\stackrel{}{y}(i)^2`$, and the matrix that is the sum of their outer products $`B=_{i=1}^n\stackrel{}{q}(i)\stackrel{}{y}(i)^𝖳`$. If $`(B^𝖳)_{kl}=B_{lk}=_{i=1}^nq(i)_ly(i)_k`$ denotes the transpose of this 3$`\times `$3 matrix $`B`$ and tr denotes the trace, then the rms distance $`d`$ is

$$d=\left\{\frac{1}{n}\left[Q^2+Y^22\mathrm{tr}\left(BB^𝖳\right)^{\frac{1}{2}}\right]\right\}^{\frac{1}{2}}.$$

### 2.3 Two Proteins

We have performed simulations on a protein fragment of 36 amino acids called the villin headpiece (1VII). We begin by rotating the $`2n`$ dihedral angles $`\varphi `$ and $`\psi `$ of the protein to $`\pi `$, except for the angle $`\varphi `$ of proline. In this denatured starting configuration, the average rms distance $`d`$ is 29 Å. Our best simulations so far fold the villin headpiece to a mean rms distance $`d`$ that is slightly less than 5 Å from its native state.

Our second protein is a 56-residue fragment of the 63-residue protein cole1 rop (1ROP). From a denatured configuration with $`d=55`$ Å, our code folds this protein to a mean rms distance $`d`$ of slightly less than 3.2 Å from its native state.

## ACKNOWLEDGEMENTS

We wish to thank Ken Dill for many key suggestions; Charles Beckel, John McIver, Susan Atlas, and Sorin Istrail for helpful conversations; Sean Cahill and Gary Herling for critical readings of the manuscript; and Sau Lan Wu and John Ellis for their hospitality at CERN. We have performed our computations on two (dual Pentium II) personal computers running Linux; we are grateful to Intel and Red Hat for reducing the cost of scientific computing.