INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Proteins are dynamic entities so that any attempt to understand their mechanism of action requires an analysis of their dynamic behavior while they are performing their allotted function, e.g., ligand binding. Having an accurate crystallographic or nuclear magnetic resonance (NMR) structure is a prerequisite for this, but it is certainly not enough. Although limited information about motions in the protein can be derived from certain types of crystallographic or NMR experiments, the picture that emerges from the use of these techniques is, to all intents and purposes, a static one. Molecular dynamics is a computational technique (McCammon et al. 1977; Brooks et al., 1983; Van Gunsteren and Berendsen 1987; Allen and Tildesley 1987) that enables the protein scientist to simulate the dynamic behavior of proteins and examine which motions are most critical for maintaining structure, unfolding, or refolding. Later applications of the technique have addressed issues such as the effect of mutations, docking to other proteins, or ligand binding (Peters et al., 1997; Peters and Bywater, 1999).

At the outset, one might consider what kind of dynamic behavior to expect from different types of protein having different functions. Some proteins like keratin in hair or silk fibroin do not bind a ligand at all but play a purely structural role. Other proteins have a more complicated function; enzymes, for example, not only bind ligands but perform chemical reactions on them. This is usually accompanied by changes in the conformation of the protein and the ligand (Peters et al., 1997; Peters and Bywater, 1999). Allosteric enzymes, which respond to changes in conditions such as varying concentration of ligands or other effectors by changing their shapes and even function, will exhibit even more complicated motions. Somewhere between the functionally simpler structural proteins and the enzymes is a class of proteins that perform only a binding or transporter role for ligands. It would seem likely that the complexity of motions within these proteins would be intermediate also. The complexity of motions would most likely correlate with the degree of complexity of the function of the protein.

To quantify the degree of complexity and to be able to compare motions between, for example, cases where different ligands are bound to the same protein, or a mutation has been made, we developed a tool, described below, for analyzing the trajectories obtained from MD simulations. In this paper, we consider the motions in members of different families of proteins, starting with the peptide-binding protein SH3, which is present as a distinct domain, or folding unit, in many intracellular proteins and which mediates protein-protein interactions important in downstream signaling and cytoskeletal function (Kato et al., 1986; Potts et al., 1988; Nemeth et al., 1989; Seidel-Dugan et al., 1992; Koch et al., 1991; Musacchio et al., 1992; Pawson and Gish 1992; Pawson and Schlessinger 1993; Cicchetti et al., 1992; Weng et al., 1993; Barfod et al., 1993; Gout et al., 1993; Liu et al., 1993; Taylor and Shalloway 1994; Fumagalli et al., 1994; Feng et al., 1994).

SH3 domains are small, water-soluble globular proteins typically constructed from about 60 amino acid residues consisting of 450 atoms. They have a compact, all- structure (SCOP fold type 21 [Murzin et al., 1995] and CATH class 2.30.30.10. [Orengo et al., 1997]) approximately spherical in shape, but with a discernible ligand-binding region on the surface.

There is a very voluminous literature associated with SH3, and correspondingly, a large number of biophysical studies have been carried out, culminating in a number of crystal and NMR structures, some in the unliganded form and others in the ligand-bound form. Unfortunately, it is not always clear from these published studies whether the peptide ligand being studied is in a neutral ("capped," e.g., by acetylation at the N-terminus or amidation at the C-terminus) form or whether it is a zwitterion, as would seem to be the case. But the latter is not a good model for the interaction between SH3 and its target protein SH2, because the stretch of polypeptide in SH2 that is involved in binding will be flanked by a continuous peptide chain (i.e., the peptide is capped). Nevertheless, we have compared the dynamics and binding of both capped and zwitterion forms.

SH3 domains have previously been the subject of extensive MD simulations (Van Aalten et al., 1996) that were aimed at the development of protein-solvent simulation methodology and for studying protein dynamics using the technique of essential dynamics (Amadei et al., 1993). The conclusion was that a full solvent simulation was the recommended procedure for all simulations of water-soluble globular proteins, and, in this work, a cubic box of TIP3 waters was used. This, of course, requires higher overheads in the calculation than when using simpler approximations, but computational techniques embodied in the MACSIMUS (Macromolecule Simulation Software) package (Kolafa 1999), along with the use of a multiprocessor computer, reduce the computation time so that inclusion of a full solvent model is not a problem.

In earlier work, the SH3 protein 1shg (Noble et al., 1993), an unliganded form of the protein was studied (Van Aalten et al., 1996). The goal of the present studies is to compare and contrast the mode of binding and dynamic behavior of different peptide ligands bound to a similar protein, to study the essential motions in the ligand and ligand-binding site on the SH3.

One very useful objective of MD and ED studies on protein-ligand complexes is to be able to decide, based on the analysis of the dynamics, where in the ligand molecule it may be advantageous to "engineer in" or "engineer out" features such as molecular bulk, or flexibility, to improve the binding constant or to alter the binding kinetics in some desired way. ED has potentially great value for this purpose because it allows the ligand engineer to focus on individual components of the dynamic behavior of the ligand.

Further, it would be very useful to be able to partition the total binding energy among fragments of the ligand, for example, residue by residue in a peptide ligand. Unfortunately, it is impossible to partition binding entropy, and hence the free energy, in this way, as has been pointed out by Mark and van Gunsteren (1994). But the internal energy can be partitioned, at least formally and for pairwise additive forces, between fragments, and we proceed to do this so that at least one component of the free energy can be studied in this way. In a future extension of this work, we shall study cases where ligand mutants have been made and their binding energies determined experimentally, while we propose to calculate the corresponding changes in the partitioned internal energy. Apart from computing trajectories for the various SH3-peptide ligand complexes, we analyze these trajectories using an essential dynamics algorithm (Ichiye and Karplus, 1991; Amadei et al. 1993) implemented within the MACSIMUS package.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Choice of starting coordinates

There exist a number of published structures for SH3 domains, both in the unliganded form and as complexes with various peptides. We have chosen to study one such complex determined by x-ray crystallography, the Ab1 tyrosine kinase SH3 domain, 1abo, which was determined to 2.0-Å resolution by x-ray crystallography (Musacchio et al., 1994). The 1abo file was edited to generate a starting structure file containing a single SH3 (chain A) and peptide ligand (chain C).

Two versions of the ligand were considered. The version in the Protein Data Bank (PDB) file is a zwitterion with the N-end positively charged and the C-end negatively. This version is given the letter Z in this work. The neutral version (referred to as N) the N-end is acetylated (CH₃-C:O-) and the C-end is substituted by methylamine (-NH-CH₃). These conversions are performed as one of the options in the PDB  MACSIMUS converting utility, which includes energy optimizations of the newly added atoms.

Force field

The forces between protein molecules are assumed to consist of the following components:

1.	Short-ranged nonbonded forces between individual atoms modeled by the Lennard-Jones potential with smooth cutoff, additive diameters, and the Slater-Kirkwood formula for the energy parameters (Brooks et al., 1983).
2.	Coulomb nonbonded forces between effective atomic partial charges qi and qj.
3.	Bond potential between a pair of bonded atoms approximated by a harmonic oscillator or a constrained bond.
4.	Bond angle bending forces between triplets of particles taking into account nonrigidity of bond angles.
5.	Torsional forces between quadruples of bonded atoms to take account of
	a.   the ability of two parts of a molecule connected by a bond to rotate around this bond (proper torsion or a dihedral potential), and b.   the need to preserve planarity of sp2 hybridized group or chirality of carbon atoms.

Apart from the Coulomb forces, the mathematical forms of these contributions are invariant across various academic and commercial implementations, but different implementations may use different numerical values of the atom-atom specific parameters appearing in those formulas that have usually been found semi-empirically elsewhere in the literature. For the Coulomb forces, the numerical values of the partial charges may likewise vary from one implementation to another. In addition, there are a number of approaches to the problem of how to deal with the long-ranged nature of the Coulomb forces. These include ignoring the problem by just cutting off the potential; introducing an effective dielectric constant into the Coulomb potential; and using truly periodic potential and replacing the bare Coulomb potential by the periodic Ewald (Allen and Tildesley 1987; de Leeuw et al., 1980) or fast multipole approximation (Greengard 1988)

To elucidate the role of the electrostatic approximations, we implemented both the cutoff version (referred to as C) and the Ewald summation (E). The particular form of the cutoff potential for MD includes smoothing and shifting,

(1)

where U, A, and B were set according to the conditions of continuous energy U_elst(r) and forces dU_elst(r)/dr.

The Lennard-Jones potential has a smooth cutoff as well (but not shifted) with standard corrections added. Because this potential decays rapidly, there may hardly be any doubt regarding the validity of this approximation in bulk liquid simulations. For details see Kolafa (1999).

Along with the Coulomb force approximations, there are several alternative strategies for dealing with the solvent (water): treating the solvent as vacuum or dielectric continuum; Brownian dynamics, where the dynamic effects of the solvent are modeled by random forces; and full molecular model of the solvent.

It is usual to speak of the set of numerical parameters appearing in the interatomic forces formulae, often with the method used to handle the Coulomb contributions, solvent, etc., used in package X as the "X force field." The values of partial charges and to some extent also the values of the parameters in the non-Coulomb forces, however, reflect the version of the Coulomb and solvent methodology used. Any comparison of results obtained using different methodology (e.g., cutoff electrostatics vs. Ewald summation) must be made with caution.

To make comparison with other works simple, we used the same values of bonded and nonbonded forces as in the CHARMM 19 force field², (hyperlink http://www.pharmacy.ab.umd.edu/~alex/, http://www.pharmacy.ab.umd.edu/~alex/) with minor modifications. Aliphatic hydrogens were considered implicitly via the extended (united) atom model for groups CH_n.

Simulation conditions

Simulations were performed on a single SH3-ligand complex in a cubic periodic box of 3000 TIP3 waters at density 1 g cm³, which guarantees a large enough box to prevent the protein from interacting with its periodic images. The start-up algorithm contained five steps: 1) energy minimization with all known atom positions fixed and optimized positions of hydrogens (not present in the original PDB files) and atoms in the acetyl and methylamine groups added to cap the termini, 2) random shooting of water molecules into the box with excluded water-protein overlaps, 3) Monte Carlo simulation of the water subsystem (i.e., with the configuration of the protein fixed) to remove water-water overlaps, 4) molecular dynamics equilibration of the water subsystem at temperature T = 300 K, and (5) assigning the atoms of the protein random velocities according to T = 300 K.

Both simulations (N and Z) were carried out twice, once with the use of the Ewald-sum method of handling electrostatics (E) and once using cutoff method (C). Optimal parameters for the Ewald summation were determined according to the error formulae by Kolafa and Perram (1992) with maximum error of force in real space of 0.024 cal/mol/Å and in the reciprocal-space 0.24 cal/mol/Å. This defines two conditions for three parameters, the remaining condition was performance optimization. The final values of parameters are real-space cutoff r_c = 18.0 Å, reciprocal-space cutoff K = 8.6, and the separation parameter  = 0.19/Å. The dielectric constant of continuum surrounding the periodic sample at infinite distance was 80. The value of the cutoff in Eqn. 1 was r_c = 12 Å.

All four simulation runs were carried out for at least 1 ns with 2 fs timesteps. Bond angles containing hydrogens and all bond lengths were constrained using the SHAKE algorithm with the Verlet integrator. The MD simulations were performed using the MACSIMUS package (Kolafa 1999). This shareware package contains code suitable for MD of large and small molecules in various boundary conditions as well as numerous tools for data analysis and visualization (X11 and DOS). It supports spatial decomposition of the simulation box using a variant of the linked-cell list method (Allen and Tildesley, 1987), parallel execution on multiple processors based on this decomposition, and efficient implementation of the Ewald summation. The most demanding calculations using the Ewald summation were performed at the Mærsk McKinney Møller Institute for Production Technology on SGI Onyx parallel supercomputer with 24 R4400 processors, whereas the version with cutoff electrostatics was performed on a standard workstation.

Principal radii of gyration

Let us consider the inertia tensor I and its diagonalization,
where r_CM is the center of mass and U is an orthonormal matrix describing the orientation of the principle axes. The eigenvalues I_diag,i may be used to define the principle radii of gyration, R_i = (I_diag,i/M)^1/2, where M is the total mass. The (total) radius of gyration is then given by
The physical meaning of R_i is that a "cloud" of mass M and Gaussian distribution of density with standard deviations R_i in the directions of the principle axes has the same inertia tensor. The three values of R_i thus give the first estimate of the shape of the molecule.

Essential dynamics

The essential dynamics (Amadei et al., 1993) method represents a powerful tool to separate the essential structural changes during the time development from thermal noise. It is based on the analysis of the positive definite, symmetric covariance matrix,

(2)

where q denotes the 3n_ess-dimensional vector of (x, y, z) coordinates of a certain subset of 3n_ess atoms chosen from the studied molecule and · is the ensemble average approximated by the average over the MD trajectory.

In principle the essential dynamics algorithm should include all atoms of the studied molecule. For simplicity, a reduced set of atoms is often used in the calculations. We have tried four different subsets of atoms: C,  carbons, for the capped (N) ligand also methyl carbons in acetyl and methylamine; B, backbone atoms; H, all heavy atoms; and A, all atoms (with the exception of aliphatic hydrogens). The coordinates, q, entering Eq. 2 are raw coordinates from the simulation only for the first configuration (frame). All subsequent frames were first translated and rotated as to superpose the protein onto the protein in the first frame. This was accomplished by minimizing the mean square distance of both configurations,

(3)

where R and T represent the operators of rotation and translation, respectively. The minimization thus takes place over six variables (three in the vector of translation T and three angles in R) and was implemented by a Monte Carlo method.

The eigenvalues _i and eigenvectors Q_i of the covariance matrix are given by

(4)

and are found numerically using the theshold Jacobi diagonalization method (Ralston 1965). The eigenvectors are orthogonal in the 3n_ess-dimensional space. Normalization factor n_ess makes comparison of results with different atom sets easier.

The eigenvectors Q_i corresponding to the highest (or several highest) eigenvalues give the direction of the essential motion relative to the center (average) position q. The quantity of physical interest is the time evolution of the projection of q  q to Q_i, i.e., the amplitude of the ith essential motion as the function of time,

(5)

The normalization factor n_ess^1/2 comes from averaging the (squared) amplitudes over all atoms.

It is useful to be able to compare two essential eigenvectors, either corresponding to different trajectories (simulations runs), or to the same trajectory with different atom sets. The simplest criterion is to compare 3D vectors of essential motions of one chosen atom. The cosine of the angle of these vectors is a measure of similarityvalues close to 1 or 1 correspond to parallel motions (note that the sign cannot be distinguished. However, it makes sense to compare signs of cosines calculated for different atoms). Before doing this, it is necessary to match the averaged configurations q for both trajectories, which is done in the same way as in Eq. 3.

Alternatively, it is possible to measure the cosine of the angle in the whole 3n_ess space,

(6)

If the sets of atoms are different for both vectors, it is necessary first to exclude superfluous atoms to obtain vectors of the same length. This criterion is much stronger than the first one because there is much lower probability that two randomly chosen vectors in a many-dimensional space are parallel. This formula can be interpreted also as the correlation coefficient of the set of all coordinates of the motion in the direction of Q_i⁽¹⁾ correlated with the corresponding coordinates of Q_i⁽²⁾. Correlation coefficients between different simulation runs may be defined in the same way. Note that we cannot distinguish Q_i from Q_i and thus the sign of c_i is irrelevant.

Energy partitioning

To cast light on the ligand-receptor binding mechanisms, we attempted to partition the binding internal energy of individual ligand residues to the receptor in three ways.

1.  The "bare" energy is a sum of all pair (Lennard-Jones plus charge-charge) energy contributions between the selected ligand residue and the receptor; residue-residue interactions and receptor-receptor interactions are not included, neither are interactions including water.

2.  The "solvated" energy is the bare energy plus the sum of all pair energies between the ligand residue and water.

3.  To determine the role of solvation, we also ran short (50 ps) trajectories of pure ligands in water (1000 molecules) and calculated the solvation energies of the ligand residues and water. Both the zwitterionic and capped forms were considered severally with the Ewald summation and cutoff electrostatics.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Analysis of trajectories

The ligand bound to the SH3 protein 1abo is the very proline-rich decapeptide APTMPPPLPP. To judge from the atom types present in the PDB file, this peptide was unprotected at the C- and N-termini, which means that the peptide is a zwitterion. The significance of choosing this particular construct may be questioned if the intention was to mimic the binding of an SH3 domain to a proline-rich stretch of polypeptide chain in, e.g., an SH2 domain in a downstream signaling complex. Such a polypeptide would be "protected" in the sense that it will be flanked by peptide residues of the rest of the main chain of the protein. As such, the peptide groups joining the decapeptide to the rest of the protein would be uncharged, whereas in the unprotected decapeptide (Musacchio et al., 1994), the N-terminal -amino group will be positively charged and the C-terminal carboxyl group will be correspondingly negatively charged. We have carried out the simulation of the peptide-protein complex in this form, even though these should presumably be capped. For comparison, we also ran the capped forms.

From the runs carried out here, a first impression is that the charged ends of the peptide ligand may have some effect on the ability of the peptide to fully bind to the protein. In both the E and C runs, the C-terminus seems to want to detach itself from the protein, and water molecules are seen to intervene between the partly detached end of the peptide and the protein. This trend was more clearly marked in the capped form. The N-terminus moves during the simulations from the vicinity of TRP110 in the protein to GLU86 (EZ run) or ASN78 (CZ run), whereas in the capped form, it releases. In the CZ run, the unprotected charged ends of the peptide get close together, which inevitably affects the whole configuration. The largest rearrangement of the peptide occurs for the EN run. This is also reflected by the mean-squared distance of the starting and ending configuration, see Table 1.

Time development of the protein structure is characterized by convergence profiles of selected quantities, which are collected, for the four trajectories considered, in Fig. 1, frames 1-4. Shown variables are total energy E; principle radii of gyration R_i; amplitudes of the first three essential motions, A₁, A₂, and A₃; and partial energies of the first and last residue E₁ and E₁₀.

View larger version (61K): [in this window] [in a new window]   FIGURE 1   Convergence profiles of selected quantities for systems EZ, EN, CZ, and CN. E, total energy in kcal/mol; Rg, radius of gyration in Å; Ai, amplitude of ith essential motion in Å (solid, backbone atoms; dotted, heavy atoms); Ei, partial residue energy (solid, bare; dotted, solvated). Sampling rates for E and Rg are 0.05 ps and averages over 1 ps are shown. Sampling rates of Ai are 1 ps and of Ei are 10 ps.

Essential Dynamics

The first question, though very technical, is how the choice of the atom set in the essential analysis affects the results. The simplest information is seen in Fig. 2, which shows the spectra of essential eigenvalues for four sets of atoms. For the lower indices, up to 5-10, the curves for _i versus i for all four sets are almost parallel. This suggests that the essential motions affect the whole backbone and are sufficiently described by the sparsest atom set C. For higher indices, curves C and B gradually begin to deviate from curves A and H, exemplifying the role of motions of sidechains in the latter case. Interesting and not understood is the decrease of _i(C) at i about 130-140. Note that the C-spectrum ends at i_end = 3(n_ess  2) = 198 for Z systems and i_end = 204 for N; it holds _i = 0 for i_end < i  3n_ess because these last six eigenvalues correspond to rotations and translations that have been filtered out. To elucidate this question, we plotted in Fig. 3 cosines of angles between eigenvectors for different atom sets, Eq. 6; all vectors entering this equation are reduced to positions of the C_a atoms.

View larger version (23K): [in this window] [in a new window]   FIGURE 2   Logarithms of essential eigenvalues sorted from the highest to the lowest. Labels A, H, B, and C denote atom sets. Solid line, system EZ; dotted line, EN; long-dashed line, CZ; short-dashed line, CN.

View larger version (37K): [in this window] [in a new window]   FIGURE 3   Absolute values of the cosines of angles, Eq. 6, between essential motions calculated using different atom sets. , pair C-B; , C-H; ×, C-A; +, B-H; , B-A, , H-A.

It is seen that the first several (2-4) essential motions are always highly correlated, in other words, the C set is sufficient. The correlations C versus B (only backbone) and A versus H (side chains) are higher than for any combination of an atom set containing only backbone (C and B) versus an atom set with side chains (A and H). For system CZ, the third most significant motion takes place in side chains.

Along with Fig. 2, numerical values of three most significant eigenvalues are also shown in Table 1. The highest value is reached for EN system. This can be interpreted such that the motions in this system occur mainly along one direction, or that one (linear) mode of motion is by far most significant. Note that, for this system also, the averaged mean distance d reaches maximum, i.e., this system is subject to the biggest structural changes.

The time development of eigenvalues and of the partial radii of gyration reveal one important feature of the simulations, namely, that the complex undergoes a gradual change of conformation rather than fluctuations from a well-defined equilibrium value. In addition, the results of all four runs seem to be uncorrelated, see Table 2 and Fig. 1. This may reflect 1) inaccurate force field, 2) inappropriate simulation methodology, 3) inappropriate ligand model, or 4) the system has not reached equilibrium. As regards 1, there is not too much space left for future improvements in the current state-of-the-art of complex molecule simulation technology. One might include explicit aliphatic hydrogens, try a different force field, or include a better water model. As regards 2, we have tested the Ewald summation. This is supposed to be more accurate than cutoff electrostatics, provided that the values of partial charges are adjusted accordingly, which is questionable. As a consequence of omitting some parts of the Coulomb interaction, the cutoff method destabilizes the structure as can be deduced from the larger fluctuations of the energy curves in Fig. 1. As for 3, the short decapeptide is an approximation of the SH2 ligand loop, which, in complexes involving the intact protein, may be much more stable. The capped version of the model peptide (presumably more appropriate as a model of the real complex) tends to detach more easily than the zwitterionic form. This may reflect the fact that SH2-SH3 complexes need to dissociate at some stage, in addition to binding. The charges at the end of the zwitterionic form introduce attractive interactions not present in the natural SH2-SH3 complex. As regards 4, which is suggested by the convergence profiles in Fig. 1, correlation times for protein motions are of the order of milliseconds. Therefore, the only way to assess the accuracy of short trajectories is to compare several independent runs.

In terms of accuracy, i.e., ability to reproduce experimental structural data, the x-ray crystal structure may differ from the structure (or rather set of dynamic structures) found in solution. First, crystal structures are a time average and second, by their very nature, the molecular structures are subject to perturbations arising from crystal packing forces. Another factor that must be born in mind when comparing crystal structures and those obtained by molecular simulations is that the crystal structures are refined using energy minimization software, which may have a different force field than the one used in the simulations.

Partitioning of the energetic contributions from each residue

The results for the four cases EZ, CZ, EN, and CN are shown in Fig. 4. The two zwitterionic cases behave very similarly, with much lower energy contributions from the terminal residues. This is most apparent when water-ligand interactions are included. When the water contribution is absent, the C-terminal end has an energy contribution more in line with the other residues. This may reflect the observation above that, in the trajectories, the C-terminus appears to want to detach itself and admit water. For the residues along the chain, only threonine, in the third position, shows any behavior that deviates markedly from the other residues, and this is only apparent when water is included. EN and CN behave similarly to each other, but, compared with the zwitterionic pair, show much more marked differences as one proceeds along the chain. When water is included, the termini show a low energy contribution and the threonine behaves as before with the leucine in position 8 also deviating somewhat.

View larger version (24K): [in this window] [in a new window]   FIGURE 4   Energy partitioning of ligand-receptor interactions. The x axis shows the residue. Solid lines, bare residue-receptor energies; dashed lines, solvated energies; and dotted lines, solvated energies minus solvated energies of the respective free ligand in water.

The dotted lines in Fig. 5 are obtained by subtracting solvation energies of the respective residues of a free ligand in water from the solvation energies of a bound ligand. Even though the solvation of a free ligand may be affected by random conformational changes (there results were obtained by independent MD runs), this difference represents "clean" binding energies of individual ligands in water environment. If we consider the size of the residues and a special role of the ends, the results support the contention that the prolines contribute most to the binding.

View larger version (101K): [in this window] [in a new window]   FIGURE 5   Stereo picture of the backbone essential motion corresponding to the highest eigenvalue for systems from the top EZ, EN, CZ, and CN. Color changes from yellow to orange (the ligand) and from cyan to green correspond to the same direction of motion. The N-termini are blue and green. The amplitude of the motion is given by the condition of the same mean square displacement as the real motion in the direction of the essential eigenvector. The left and right figures are for the left eye and the central figure for the right eye. Watch the left pair with parallel eyes (the cross-point is behind the paper) or the right pair with eyes crossed at a point between the paper and eyes.

We have not yet investigated what happens if residues in the peptide are mutated. This kind of simulated experiment will be carried out in other peptide-protein systems where there is experimental data (binding energies) for the binding of different mutants. Some data of this kind has recently been published by Nguyen et al. (1998), who showed that not only proline, but other N-substituted residues, can bind very efficiently to SH3. It would be of interest to include these mutants in future simulation studies.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Analysis of trajectories

The ends of the peptide tend to bind more loosely in the capped form than in the zwitterionic form, due to the lack of the electrostatic interactions operating in the latter case. Inasmuch as we surmise that the capped form is more relevant to the real situation where SH3 recognizes SH2, this suggests that the end residues in the peptide fragment are not important for binding, whereas the middle residues, including proline, are responsible for most of the binding energy.

Essential dynamics

It is possible to observe significant motions in the bound ligand and to distinguish between motions of the peptide backbone and of the side chains. This could be useful in designing ligands, which fit optimally to the binding protein. The most favorable model for observing these motions is probably EN, i.e., the capped peptide studied using Ewald summation. From the practical point of view, it is important that the first essential motions are sufficiently well described by the atom set reduced to the C atoms only.

Partitioning of the energetic contributions from each residue

We show that it is possible to determine the different contributions of each residue in a peptide to the enthalpy of binding. Proline is probably the major contributor to binding enthalpy, in keeping with the observed propensity for this family of proteins to bind proline-rich peptides.