A simple way to compute protein dynamics without a mechanical model

SHORT COMMUNICATION

A Simple Way to Compute Protein Dynamics Without a

Mechanical Model

Chien-Hua Shih, Shao-Wei Huang, Shih-Chung Yen, Yan-Long Lai, Sung-Huan Yu, and Jenn-Kang Hwang* Institute of Bioinformatics, National Chiao Tung University, HsinChu 30050, Taiwan, Republic of China

ABSTRACT We found that in proteins the av-erage atomic fluctuation is linearly related to the square of the atomic distance from the center of mass of the protein. Using this simple relation, we can accurately compute the temperature factors of proteins of a wide range of sizes and folds, and the correlation of the fluctuations in proteins. This sim-ple relation provides a direct link between protein dynamics and the static protein’s geometrical shape and offers a simple way to compute protein dynam-ics without either long time trajectory integration or any matrix operations. Proteins 2007;68:34–38.

V

VC 2007 Wiley-Liss, Inc.

Key words: protein dynamics; thermal fluctua-tions; molecular dynamics; normal mode analysis

INTRODUCTION

The computation of protein dynamics is usually based on a mechanical model. For example, molecular dynam-ics computes the protein’s trajectory based on a sophisti-cated molecular mechanical model of bond stretching, bond angle bending, bond twisting, van der Waals and electrostatic interactions.1–3The recently developed elas-tic network model,4–6 which has been successfully applied to analyzing large-scale protein motion, is based on a much simpler mechanical model that each atom in proteins are connected to its surrounding atoms that are within a certain cut-off distance by a one-parameter har-monic spring. Here we report a simple method to com-pute protein dynamics directly from the static protein geometrical shape without any mechanical models. Our method is based on the observation that the deeper an atom is buried inside a protein structure, the less it will fluctuate around its equilibrium position. We found that this observation goes beyond a mere qualitative descrip-tion. We found that in proteins the atomic fluctuation is in fact linearly related to the square of the atomic dis-tance from the center of mass of the protein. Using this simple relation, we can accurately compute the tempera-ture factors and the correlation of the fluctuations in proteins.

METHODS Center of Mass Distance

Let X0 be the center of mass of the protein, that is,

X0¼

P

kmkXk=

P

kmk, where mk and Xk are the mass

and the crystallographic position of atom k, respectively. The distance of atom i from the center of mass of the protein is computed by

r2

i ¼ ðXi
X0ÞðXi
X0Þ; ð1Þ

Each protein of size N will have its distinct distribution given by r2

1; r22; . . . ; r2N



, referred to as the r2 profile. In this work, we computed the r2

-profiles of the Ca atoms of several proteins of sizes ranging from 54 to 736 and different folds including all-a, all-b, aþb, and a/b folds. The descriptions of these proteins are listed in Table I. We will show in later sections that the r2 profile is in fact closely related to the temperature factors.

Correlation Between Atoms in Proteins

We can generalize Eq. (1) to compute the correlation between the center of mass distances of atom i and j, that is,

cij¼ ðXi
X0ÞðXj
X0Þ; ð2Þ

Note that when i¼j, cij¼ r2i. In other words, the

autocor-relation reduces to the square of the distance from the center of mass of the protein. We will show later that the cijis closely related to the correlation between

fluctu-ations of atom i and j, which is given by

Cij¼ dX idXj   ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dXi dXi h i dXjdXj   q

Grant sponsor: National Science Council, Taiwan.

*Correspondence to: Jenn-Kang Hwang, Institute of Bioinfor-matics, National Chiao Tung University, HsinChu 30050, Taiwan, ROC. E-mail: jkhwang@faculty.nctu.edu.tw.

Received 4 December 2006; Accepted 16 January 2007

Published online 13 April 2007 in Wiley InterScience (www. interscience.wiley.com). DOI: 10.1002/prot.21430

V

where dXiand dXjare the fluctuations of the atom i and

j, respectively, around their equilibrium positions. The correlation between fluctuations can be computed using the normal mode analysis,7–9which requires the evalua-tion of the second derivative matrix of the total potential function and the diagonalization of the matrix. In this work, the Enzymix10 force field is used for the normal mode analysis.

RESULTS Temperature Factors

Figure 1 compares the computed r2_{profile and the}

cor-responding crystallographic temperature factors (or ther-mal B factors) of these proteins. The agreement between them is excellent. The good agreement is surprising since our method [i.e., Eq. (1)] is very simple and does not requires either a long time trajectory integration11_or

any matrix operations.4–6

Since the B factor is given as Bi¼ (8p2/3)hdXidXii, our

results suggest the following interesting relation,

dXidXi

h i  ðXi
X0ÞðXi
X0Þ ð3Þ

Equation (3) provides a direct link between protein’s dy-namics properties and its static geometrical shape. On the practical side, Eq. (3) also offers a very simple way to compute the temperature factors of proteins.

Correlation Maps

Figure 2 compares the correlation maps computed by Eq. (2) and those computed using the normal mode

analysis. Again, the agreement between cij and Cij is

excellent.

DISCUSSION

The dynamic properties of a protein result from a com-plex network of various interactions; however, our results indicate that they can be computed directly from the pro-tein’s geometrical shape without assuming any mechanical models. Now, the question is why equations as simple as Eqs. (1) and (2) will work for a molecule as complex as a protein. It will be instructive to compare our method with the recently developed elastic network model.4–6 The elas-tic network model assumes that the fluctuations of an atom are related to the positions of its surrounding atoms that are within a certain cut-off distance. This is in certain way similar to our center of mass distance—fluctuation relation [i.e., Eq. (3)], since intuitively the deeper an atom is buried inside the protein, the more number of atoms it will be surrounded, and the less it will fluctuate. The elas-tic network model assumes a mechanical model in which the atom and its surrounding atoms are connected to each other through a single-parameter harmonic potential. To compute the temperature factor, the elastic network model inverts a N3N connectivity matrix (or the Kirchhoff ma-trix), where N is the size of the protein. This matrix is con-structed based on the positions of the surrounding atoms of each Ca atom in the protein.

Originally, it was rather surprising that a simple method like the elastic network model could describe protein dy-namics so well4–6 _{when compared with other more}

sophis-TABLE I. Information of Proteins Studied in this Work

PDB IDa _{Length Protein name Fold}b

1PD3:A 54 Influenza A NEP

M1-binding domain

ROP-like (a)

1U0S:A 87 Chemotaxis kinase CheA P2

domain

Ferredoxin-like (aþb)

1VJH:A 118 At1g24000 TBP-like (aþb)

1MIJ:A 148 Homeo-prospero Domain DNA/RNA-binding 3-helical

bundle (a)

1F35:A 157 Olfactory marker protein Olfactory marker protein (b)

1WUB:A 178 Polyisoprenoid-binding protein Hypothetical protein (b)

1B12:A 239 Singnal peptidase I LexA/signal peptidase (b)

1IOW 306 D-Ala-D-Ala ligase PreATP-grasp domain and ATP

grasp (a/b)

1NVM:A 340 Bifunctional

aldolase-dehydrogenase

TIM barrel and RuvA C-terminal domain-like-(a/b)

1DKQ:A 410 Phytase Phosphoglycerate mutase-like (a/b)

1FUO:A 456 Fumarase C Orthogonal bundle and up-down

bundle (a)

1FCE 629 Endocellulase CelF a/a toroid (a)

1LF9:B 674 Glucoamylase (1AYX) Six-hairpin glycosidases and

galactose mutarotase-like (a)

1J7N:B 736 Anthrax toxin lethal factor Zincin-like and ADP-ribosylation

(aþb)

a_{The first four letter is the PDB ID. The fifth letter (if any) after the colon mark is the chain designator.} b_{The definition of the fold follows that of SCOP.}

35

Fig. 1. The computed (bold line) and the crystallographic (dotted) temperature factors (orB factors) of Ca atoms of 1PD3:A, an influenza A NEP M1-binding domain with a ROP-like all-a fold; 1U0S:A, a chemo-taxis kinase CheA P2 domain with a ferrodoxin-like aþb fold; 1VJH:A, a hypothetical protein (At1g24000) with an aþb fold; 1MIJ:A, a homeo-prospero domain with a 3-helical bundle; 1F35:A, an olfactory marker protein with an all-b fold; 1WUB:A, a polyisoprenoid-binding protein with an all-b fold; 1B12:A, singnal pepti-dase I with an all-b fold; 1IOW, ligase with an a/b fold; 1NVM:A, bifunctional aldolase-dehydrogenase with a TIM barrel (a/b) fold; 1DKQ:A, phytase with an a/b fold; 1FUO:A, fumarase C with an all-a fold; 1FCE, endocellulase CelF with an all-a fold; 1LF9:B, glucoamylase with an all-a fold; 1J7N:B, anthrax toxin lethal factor with an aþb fold. The computed B factors are normalized to the scale of the crystallographic B factors.

ticated approaches such as molecular dynamics11 and the normal mode analysis.7–9Hence, it is even more surprising that an even simpler model [i.e., Eqs. (2) and (3)] can be developed, which does not assume any mechanical model, does not have any adjustable parameters, and does not require any matrix operations.

Both the elastic network model and our method do not require the knowledge of amino acid sequence. This suggests that protein dynamics of a folded protein are dictated mainly by its folded structure rather than by its chemical properties. Our results further suggest that the folded protein behaves like a rotating sphere Fig. 2. The cross-correlation maps for some of the proteins examined in Figure 1. For each protein, the map on the left is computed by Eq. (2) and the map on the right by the normal mode analysis coupled with energy minimization. The colors of the map ramp from red (positive correlation) to blue (negative correlation).

37

around its center of mass such that the angular de-pendent components of its dynamic properties (such as the temperature factor or the correlation of fluctua-tions) are averaged out. As a result, they can be well approximated by a function depending only on the atomic distances from the center of the protein. To fur-ther illustrate this point, Figure 3 shows a radar plot overlapping the r2 profile and B values of chemotaxis kinase CheA P2 domain with a aþb domain. The excellent overlap between these two indicates the atoms with the similar radius from the center of the protein will generally have similar thermal fluctua-tions (or B factors). In the real three-dimensional space, it is as if that there are many concentric spheres, all of which are centered on the center of mass of the protein but are of different radii, and that the atoms lying on the same surface will tend to have the same thermal fluctuations, regardless of their chemical properties.

With the growing number of protein structures depos-ited in PDB, our method may provide a very efficient way to systematically characterize protein dynamics in the study of the structure-dynamics relationship of proteins.

REFERENCES

1. Levitt M, Warshel A. Computer simulation of protein folding. Nature 1975;253:694–698.

2. Warshel A. Bicycle-pedal model for the first step in the vision process. Nature 1976;260:679–683.

3. McCammon JA, Gelin BR, Karplus M. Dynamics of folded pro-teins. Nature 1977;267:585–590.

4. Tirion MM. Large amplitude elastic motions in proteins from a sin-gle-parameter, atomic analysis. Phys Rev Lett 1996;77:1905–1908. 5. Bahar I, Atilgan AR, Erman B. Direct evaluation of thermal

fluctuations in proteins using a single-parameter harmonic potential. Fold Des 1997;2:173–181.

6. Ming D, Kong Y, Lambert MA, Huang Z, Ma J. How to describe protein motion without amino acid sequence and atomic coordi-nates. Proc Natl Acad Sci USA 2002;99:8620–8625.

7. Brooks B, Karplus M. Harmonic dynamics of proteins: normal modes and fluctuations in bovine pancreatic trypsin inhibitor. Proc Natl Acad Sci USA 1983;80:6571–6575.

8. Levitt M, Sander C, Stern PS. Protein normal-mode dynamics: trypsin inhibitor, crambin, ribonuclease and lysozyme. J Mol Biol 1985;181: 423–447.

9. Kidera A, Go N. Normal mode refinement: crystallographic refinement of protein dynamic structure. I. Theory and test by simulated diffraction data. J Mol Biol 1992;225:457–475. 10. Lee FS, Chu ZT, Warshel A. Microscopic and semimicroscopic

calculations of electrostatic energies in proteins by the Polaris and Enzymix programs. J Comput Chem 1993;14:161–185. 11. Levitt M. Molecular dynamics of native protein. I. Computer

simulation of trajectories. J Mol Biol 1983;168:595–617. Fig. 3. The radar plot of ther2-profile (black line) andB values (gray line) of the Ca atoms of chemotaxis

kinase CheA P2 domain (1U0S:A). The residues are numbered on the outer rim of the radar in a clockwise direction and theB values are numbered in the radial axis. The r2_{values are normalized to the scale of the}