Molecular dynamics simulations of native and substrate-bound lysozyme

J. Mol. Biol. (1986) 190, 455-479

Molecular

Dynamics Simulations of Native and Substrate-bound Lysozyme

A Study of the Average Structures and Atomic Fluctuations Carol Beth Postl, Bernard R. Brooks’, Martin Karplus’ Christopher M. Dobson2, Peter J. Artyrniuk3 Janet C. Cheetham3 and David C. Phillips3 Harvard

1 Chemistry Department University, Cambridge, MA 02138,

2 Laboratory of Molecular Biophysics Oxford University, Oxford OX1 3PS,

Department

of Zoology,

3 Inorganic Chemistry Laboratory Oxford University, Oxford OX1 3&R, (Received

U.S.A. U.K.

U.K.

12 August 1985, and in revised form 19 February

1986)

Molecular dynamics simulations of hen egg-white lysozyme in the free and substrate-bound states are reported and the nature of the average structures and atomic fluctuations are analyzed. Crystallographic water molecules of structural importance, as determined by hydrogen-bonding, were included in the simulations. Comparisons are made between the dynamics and the X-ray results for the atomic positions, the main-chain and side-chain dihedral angles, and the hydrogen-bonding geometry. Improvements over earlier simulations in the potential energy function and methodology resulted in stable trajectories with the C” co-ordinates within 1.5 A of the starting X-ray structure. Structural features analyzed in the simulations agreed well with the X-ray results except for some surface residues. The Asx x2 dihedral distribution and the geometry of hydrogen bonding at reverse turns show differences; possible causes are discussed. The relation between the magnitudes and time-scales of the residue fluctuations and secondary structural features, such as helices /?-sheets and coiled loops, is examined. Significant differences in the residue mobilities between the simulations of the free and substrate-bound states were found in a region of the enzyme that is in direct contact with the substrate and in a region that is distant from the active-site cleft. The dynamic behavior of the structural water molecules is analyzed by examining the correlation between the fluctuations of the water oxygens and the lysozyme heavy-atoms to which they are hydrogen-bonded.

1. Introduction

enzyme and enzyme-substrate complexes. It is known that the binding of cofactors, substrates and inhibitors can alter the enzyme conformation. An understanding of both the structural rearrangements and changes in motional properties induced by the binding is therefore important for the analysis of enzymatic reactions. Of particular interest is the relation of the fluctuations to the system entropy, which can have an important role in substrate binding, in the relative free energies of species along the reaction path, and in product dissociation (Sturtevant, 1977). One of the approaches to the study of motional properties of proteins is provided by molecular dynamics simulations. The extensive experimental results available on the structure and the

One of the current interests in the study of proteins is the characterization of the internal motions (Karplus & McCammon, 1983; Petsko & Ringe, 1984). Since the native structure is stabilized by a large number of weak non-bonded interactions, substantial fluctuations in the atomic positions are expected to occur. Analysis of the fluctuations can provide insight into the nature of these interactions. Variations in the protein of properties such as hydrophobicity or hydrogen-bonding polarity, density modulate the stability of different regions and can lead to differences among residues in the fluctuation amplitudes and time-scales. Moreover, any attempts to interpret enzymatic catalysis require knowledge of the dynamic properties of the 455 0022%2836/86/150455-25

$03.00/O

0 1986 Academic

Press Inc. (London)

Ltd.

C. B. Post et al. enzymology of hen egg-white lysozyme make it an attractive enzyme for study by such simulations. Crystallographic determinations have been made for different crystal forms (Blake et aZ., 1965; Hodsdon et al., 1974; Kurachi et al., 1974; Moult et al., 1976; Hogle et aZ.. 1981) and at different temperatures (Berthou et al., 1983), as well as for lysozyme with various inhibitors (Ford et al., 1974; Kelly et al., 1979; Blake et aZ., 1967). A refined, high-resolution structure exists for free lysozyme that includes refined positions of interior and surface waters (Handoll et al., unpublished results). In the solution state. structural and dynamic information has been obtained by nuclear magnetic resonance studies (Poulsen et aZ., 1980; Olejniczak et al., 1981, 1984a; Delepierre et al., 1984), and the binding properties and enzymatic hydrolysis rates of numerous glycoside substrates and inhibitors have been measured (Banerjee et al.. 1975; Schindler et al.. 1977; for references, see also Imot’o et al.. 1972). In this paper, we describe the results of molecular dynamics simulations of free lysozyme and the substrate-bound enzyme. The results for free lysozyme provide information on the differences in mobilit,y of various regions in the molecule. An analysis is made of how the structural stability, as evaluated by the extent of deviation between the X-ray and the dynamics structures, and the positional fluctuations vary within the molecule, and of how their variations relate to secondary structural features such as helices, /?-sheet and loop regions. Lysozyme is of particular interest in this regard, since all three types of protein structural elements are present (Fig. 1). The dynamic properties of lysozyme with and without the bound substrate hexa-N-acetylglucosamine, (GlcNAc),?, are compared. Examination of modifications in the enzyme and substrate conformations that occur durmg the simulation suggest a possible reaction pathway for the enzymatic catalysis (Post & Karplus, 1986). A more detailed study of t’he lysozyme-(GlcNAc), structure and dynamics, and the implications for the reaction mechanism, is in progress. Three molecular dynamics simulations of lysozyme are described. The present work complements an earlier simulation of lysozyme (Ichiye et al., 1986), which has been used to analyze experiments concerned with fluorescence anisotropy (Tchiye & Karplus, 1983), nuclear magnetic resonance relaxation (Olejniczak et al., 19846) and nuclear magnetic resonance coupling constant-s (Hoch et aZ., 1985). The simulations reported here were performed with a version of the CHARMM program (Brooks et al., 1983) modified for use on the CRAY-IS supercomputer. Several aspects of the simulations are improvements over the earlier

t Abbreviations used: GlcNAc, N-acetylglucosamine; r.m.s., root-mean-square; CM, center of mass; MurKAc, N-acetylmuramic acid.

simulation. The functional form used for hydrogen bonding (Reiher 62 Karplus. unpublished results) has been modified, as have the pararneters for \-an der Waals’ interactions (States: 1984) and the functions for smoothly t,runcating the non-bonded interactions to obtain a reduced number of non-bonded pairs (Brooks et al., 1983). Polar hydrogen atoms were explicitly included to ac:commodat,tx the new hydrogen-bond pot,ential: the extended atom model for incorporating the hydrogens int,o thr heavy-atoms t!o which they are bonded was used only for aliphat’ic hydrogen at’oms. The &art)ing structures for the simulations are the recentI)refined and higher resolution X-ray co-ordinattss for free lysozyme (Handoll rf a,l.. unpublished results) and the (GlcNAc),-lysozyme complex (Art,ymiuk Pt al., unpublished results). Furthermore. water molecules that appear to play a struc>t,ura,i role werta included in the simulation (Blake et a./.. 1983). Sillc>ca the substrate replaces a number of’ these water molecules, their presence in the free lysozyrne simulation allows a more nieaningf1~1 comparison ot the free and substrate-bound trajecatories. as well as providing hydrogen-bonding interactions that, ar(b significant for stabilizat’ion of t)hr ?i-ray nat ivtb conformation. Finally, the lenggh 01‘ two of’ t)he simulations is greater (55 ps and .I 01 px. cvmpared with 33 ps for the earlier simulation). so that thchrti is a better basis for the exarninat,ion of structural and dynamic properties. Most of the results presented here are from ttw analysis of two of the simulations: a 101 -picosecond trajectory of free lysozyme with 53 water molecules and a 55.picosecond trajectory of bound lysozyrnr t,he substrat,e with 35 water molecules plus (GlcNAc),. The third simulat,ion is a 27.picosecond trajectory of free lysozyme. using the same initial co-ordinates as those used for the 101 -picosec~ond trajectory. Since an isolat,ed (va,cuum) system without boundaries or a crystalline cnvironmcni was used in the simulat,ions. a constraining potenlJial was applied t,o the water oxygen aboms to prevc’nt the water molecules from separating frorn tht> as found in an early bovine panureatic~ protein, trypsin simulat,ion (Levitt. 1980). The 27.picosecond simulation of free lysozyme wa,s c~alculated using an alternative potential for constraining tlitb molecules. The purpose of this short water trajectory was to compare the two methods used t’o constrain the water molecules. We begin by describing the computational procedure followed in the molecular dynamics calculations; a brief discussion of t,he potential function and starting co-ordinat,es is given. ‘l’hr results of the simulations are then examined by first comparing the average dynamics struct.ures with t,tlcA starting X-ray st’ruct’ures. i’iec‘ond. the atomic. fluctuations in the simulations are analyzed and the mobilities of the various structural regions of lysozyme are compared. Differences found in thr motions of lysozyme between the unbound and bound stat,es are outlined and their possible significance is discussed.

Molecular

Dynamics

Simulations

2. Methods (a) Simulations The simulations were carried out on the CRAY-1S computer, at the SERC Daresbury Laboratory, England, with a modified version of the CHARMM program (Brooks et al.. 1983). For 1531 atoms, with approximately 1.1 x lo5 non-bonded interactions and 1200 hydrogenbond pairs, 1 ps (1100 integration steps for a time-step of 0.98x lo-“’ s) required 15 min of CRAY-1S central processing unit time. The empirical energy function used to obtain the forces on the atoms includes harmonic potentials that govern changes in bond lengths and angles, an improper dihedral potential, a 2 and 3-fold torsional potential, a van der Waals’ and a coulombic potential, and a hydrogenbonding potential. The methodology and the general form of the potential functions were as described (Brooks et al., 1983). Special aspects of the potential functions used in the reported trajectories are outlined below. The form of the hydrogen-bond potential (Reiher & Karplus, unpublished results) was as follows:

x sw(r~~, where the switching

rod eaa-a-H. cos* Wa-“-‘&

function,

(1)

SW, has the form:

sww = (x,,,-x)*(x,,,+~~-~X,,)/(X,,,-X,,)~; A and B are the parameters for each hydrogen-bond type; a, aa and d are heavy-atom acceptor, acceptor antecedent and donor, respectively; and rad is the a,d distance. The well depth for N,N pairs, N,O pairs and 0,O pairs are 3.85, 4.0 and 4.25 kcal/mol, respectively (1 cal = 4.184
sf Lysozyme

457

0.98 x lo- l5 s. In calculating the trajectories, the lengths of the bonds involving hydrogen atoms were constrained according to the SHAKE algorithm (Ryckaert et al., 1977; van Gunsteren & Karplus, 1980). All other degrees of freedom of the atoms treated explicitly were included. Lists containing the non-bonded and the hydrogenbonded pairs were determined every 0.02 ps and 0.01 ps, respectively, using a cutoff distance of 9-O A and 7.0 A for the non-bonded and hydrogen-bonded lists, respectively, and a minimum angle of 80” for the hydrogen-bond list. These update frequencies and cutoff limits conserved energy within 0.01 cal/mol between updates of t’he lists. No adjustments of the atomic velocities were required to maintain a constant energy. The water molecules were modeled according to the ST2 formulation of Stillinger & Rahman (1974), except that the oxygen-oxygen van der Waals’ interactions were switched off at long range as well as at short range. To maintain consistency with the ST2 model for water, explicit hydrogen-bond interactions were not included between the water molecules and the protein (Rossky et al., 1979; Rossky & Karplus, 1979). In addition, the dielectric constant used in the protein-ST2 electrostatic terms was set equal to 1.0 to have nearly equal potential well depths for the protein-protein hydrogen-bond interactions and the ST2-protein interactions. To eliminate the possibility of a water molecule becoming separated from the protein during the trajectory, a constraining potential, referenced to the initial position of the oxygen, was placed on each water molecule. The form of this potential Tires,, was chosen to be a quartic function in the oxygen position, R,(t), relative to its initial position R,(O):

f’r,,, = k,[Ro@)- Ro(0)14.

(2)

Use of a quartic dependence provides a nearly constant force at distances R,(t)-R,(O) I 1.0 ip and a steeply rising constraint at distances R,(t) - R,(O) > 1.0 A. In the 101 ps trajectory of free lysozyme and the 55 ps trajectory for substrate-bound lysozyme, the value of the force constant was chosen to give an energy of 1.0 kcal/ mol at a distance of 3.0 A, this being 3 times the r.m.s. displacement calculated from the average crystallographic Debye-Waller factor of the 53 water oxygen atoms (Handoll et al., unpublished results). With this small value for k,, the constraining potential has little influence on the water dynamics until a water molecule has moved nearly 3 A from its crystallographic position. The third trajectory, a 27 ps trajectory of free lysozyme, utilized a different potential to constrain the water molecules in the neighborhood of the protein. The reference position for the quartic potential was the center of mass of the system rather than the initial oxygen position. This trajectory of free lysozyme will be referred to as the CM native trajectory, distinguishing it from the 101 ps native trajectory. The CM constraining potential was applied to all atoms in the systems and used a force constant resulting in a 0.3 kcal/mol energy for an atom 30 A from the center of mass. The most, &stant protein atom in the initial structure was approximately 28 A from the center of mass. (b) Starting

co-ordinates

The initial heavy-atom co-ordinates for the native trajectory are the crystallographic co-ordinates of lysozyme refined to 1.6 A resolution (Handoll et al., unpublished results). Polar hydrogen atoms were included explicitly in the simulations. The co-ordinates for these

C. B. Post et al. atoms were determined using standard geometric parameters, combined with a search for a minimum-energy site in the case of non-unique configurations (Brunger & Karplus, unpublished results). From the 150 water molecules determined in the X-ray structure, 53 were included in the simulation because of their apparent structural importance; a water molecule was included if it had at least 2 hydrogen bonds to the protein either directly or through intermediate water moledules. Fig. 2(a) shows the backbone of lysozyme and the location of the 53 water molecules. The water molecules were included both to help maintain structural integrity and to give a more reliable comparison between the motional properties of free lysozyme and those of the substrate-bound enzyme, i.e. certain water molecules present in the active site of the free enzyme are displaced by the substrate. For the bound trajectory, the starting co-ordinates are the crystallographic co-ordinates of the complex between lysozyme and the trisaccharide inhibitor (GlcNAc),, plus model-built co-ordinates for the remaining 3 sugar molecules. The initial stage of the crystal structure determination was to calculate an electron density map for the complex with phases obtained from the 2-O A resolution, refined co-ordinates of free lysozyme (Grace, 1980) and the diffraction intensities from the enzymeinhibitor complex. The co-ordinates of lysozyme plus (GlcNAc), were built into this map and then refined by the Hendrickson-Konnert method (Hendrickson, 1985) using an estimated occupancy of 0.7 for (GlclC’Ac), (Atymiuk et al., unpublished results). The cis conformer of the acetamido side-chain was mistakenly built into the electron density at site A in the active-site cleft, and remained cis throughout the 55 ps trajectory. Since there are few contacts with the enzyme for either the methyl or the carbonyl oxygen and no hydrogen bond, the incorrect geometry is not thought to be of consequence to the analysis presented here. From the full set of water molecules present in the crystal structure, the 35 water molecules for the bound trajectory were chosen by the same criterion as that used for the native trajectory. This resulted in an almost 1 : 1 correspondence between the water positions of t,he bound and unbound trajectories. except for the water molecules present in the native trajectory that were excluded by the hexamer substrate (Fig. 2). To determine co-ordinates for the sugar molecules in sites D, E and F, a GlcNAc monomer was built into each site using a computer graphics system. Starting in site D. a Glch’Ac monomer in a regular chain conformation was added with a B-linkage to the terminal oxygen atom of t,he sugar in site C. The bonds of the glycosidic linkage were rotated until the sugar fitted the site without unreasonably close contacts with the protein. Further fitting was done by rotating the hydroxyl and acetamide side-groups of the sugar to optimize hydrogen-bond formation. The sugars in sites E and F were built sequentially from site D in a similar fashion. Removal of bad contacts in some cases required rotation of amino acid side-chains; no rotations of backbone dihedral angles were performed. The side-chains of Glu35, Asn44, Asn46, Asn59 and Vall09 were rotated from the free structure by the values given in Table 1. Building the hexasaccharidelvsozyme complex by allowing rotation of amino acid sjde-chains did not require distortion of the chair conformation of the ring in site D. Favorable interactions at all sites were found without introducing steric strain in the substrate (see Results). Pincus & Scheraga (1981). using a different approach involving a combination of

Table 1 Side-chain dihedral angle changesmade to build (GlcNAc) into sites D, E and F

Residue Glu35

Asn44

Native

Bound

Difference (‘1 --.

Xl x2 x3

-i7 -60 -- 15

-67 -82

10 -22

x1 x2

Asn46 Asn59 Va1109

Angle (“)

Dihedral angle

x2 x2 XI

- 24

-. 9

-72 106

-143 -54

- il - 160

IO6

--ti4

--- 170

-94

55

-127

- 33

-176

LPI)

conformational search plus minimization, also predicted a low-energy binding conformation in which there is no ring distortion. With respect to the orient,ation of ClcNAc in sit’es E and F for the different enzyme--substrat,e models. the structure used in this work is intermediate bet,weeu the original wire-model structure (Blake et nl.. 1967) and the structure report’ed by Pincus & Scheraga (1981). In t,he original wire-model structure. the 2 terminal sugar molecules lie closer to the right, side of the clrft. making contact with residues 37 and 114. while the complex determined by Pincus & Scheraga interacts with residues 45 and 46 on the far left side of the cleft (Fig. l(a) and (b) of Pincus & Scheraga? 1981). As shown in Fig. 3, the complex used in this work falls bet,ween these 2 possibilities, interacting wit,h residues 34. 35 and 44. A crystallographic structure of another t,risaccharide lysozyme complex. in which the sugars occupy sites B. C’ and D. also indicates that binding in site D occurs without distortion (Kelly et al., 1979). That (GlcNAc), can bind to lysozyme without steric strain makes the possibility of a planar distortion of the saccharide ring contributing to the rate of enhancement of the chemical reaction unlikely. lt is possible, however, t,hat t,he transition state could involve some distortion of the ring in site D away from the chair conformation (Imoto PI ~1.. 1972). The original enzyme-substrat’e model (Phillips. 1966) was concerned with the rxploration of possible functional contacts between the sugar in site D and the potentially catalytic residues Glu35 and Asp52. and involved deeper penetration into site D. The dist,ortion of the sugar molecule in this site is consistent) wit,h t,hr crystallographically observed complex involving tetrasaccharide-lactone inhibitor (Ford it nl.. 1974). With the 5 center ST2 model (Stillinger & Rahman. 1974) for the water molecules and the protein polar hydrogen atoms capable of hydrogen bonding. the number of centers included in the native trajectory was 1531 (1001 prot,ein heavy-atoms. 265 protein hydrogen atoms and 265 ST2 centers) and in the bound trajeator!, was 1546 (1001 protein heavy atoms. 265 protein hydrogen atoms, 85 substrate heavy atoms. 20 subst,rate hydrogen atoms and 175 ST2 centers). (c) Computational

procedure

The X-ray co-ordinates of free lysozyme were adjusted to remove large initial forces calrulat,ed from the empirical energy function used for the molecular dynamics and to improve the positions of the hydrogen at,oms built into the X-ray structures. Fewer than 15

Molecular

Dynamics

Simulations

heavy-atom contacts, excluding hydrogen-bonded pairs and 1,4 pairs, were closer than -3.0 A in the X-ray structure. To relax the water molecule positions, 30 steps of steepest, descent minimization of the energy were performed with the protein atoms fixed. Twenty cycles of conjugate-gradient minimization were done subsequently, with weak harmonic constraints applied to all atoms, followed by an additional 60 cycles of conjugate-gradient minimization without constraints. The r.m.s. gradient in the potential energy and the potential energy for the initial free structure were 18 kcal/A per mol and 700 kcal/mol, respectively, and were 0.9 kcal/A per mol and - 1100 kcal/mol after minimization. The r.m.s. difference in the atomic co-ordinates of the 1001 lysozyme heavy-atoms before and after minimization was 0.37 A (Table 2). Using a similar minimization procedure, corresponding results were obtained for the starting coordinates for the bound lysozyme trajectory. For the 1001 protein heavy-atoms of bound lysozyme, the r.m.s. difference between the initial X-ray co-ordinates and the optimized co-ordinates was 0.34 A. To reach a temperature near 300 K, the temperature was gradually increased from zero by 5 deg. K increments at intervals of 0.06 ps (61 steps) for the native trajectory and 0.09~~ (92 steps) for the bound trajectory. A particular temperature was obtained by randomly

459

of Lysozyme

assigning individual velocities from a Gaussian distribution with a variance corresponding to the desired temperature. After heating to the final temperature of 300 K, the system was allowed to equilibrate, whereby the temperature of the system averaged over a period of 05 ps was checked, and if the temperature was more then 5 deg.K away from 300 K, the velocities were reassigned randomly. The equilibration continued in this manner until no further drifts in the temperature occurred. The native trajectory required a 17 ps heating and equilibration period with a final temperature averaging 295 K, while the corresponding values for the bound trajectory were 28 ps and 304 K. During the equilibration, the 2-acetamido group of the sugar molecule in site B underwent a configurational change, moving from the equatorial position to the axial position. Such inversion at a chiral center (ring C-2 atoms) is due to exclusion of the bonded hydrogen atom, although this drawback of the extended-atom model could have been circumvented by inclusion of an improper dihedral potential term (Brooks et al., 1983). The analysis period of the simulations began after the last temperature adjustment. The trajectories were run without further interference (there was no scaling of energies or velocities). The ratio of the fluctuations in the

Table 2 Root-mean-square co-ordinate diferences of lysozyme X-ray structures, optimized structures and average dynamic8 structure8 Structures

Atoms

(sample

number)

r.m.s.

difference

Native lysozyme 101.ps

average

nera~a X-ray

Heavy-atoms (1001) C” (129) N, C”, c (387) Side-chain atoms1 (614)

2.2 15 1.5 2.5

Heavy-atoms C” N, C”, c Side-chain

atoms

2.2 1.5 1.4 2.6

atoms

2.0 1.4 1.4 2.3

atoms

0.9 0.3 0.3 1.1

atoms

1.8 1.2 1.2 2.2

atoms

0.37 0.25 0.26 0.43

atoms

0.34 0.22 0.22 0.39

CM native

lysozyme 27.ps average WLWS X-ray

Bound lysozyme 55-ps average

tzr~ua

X-ray

X-ray lysozyme Native UWRM bound

Dynamics Native

lysozyme V~MUS bound

Native lysozyme Optimized versus

X-ray

Bound lysozyme Optimized versus X-ray

Heavy-atoms C” N, c”, C Side-chain Heavy-atoms C N, c”, C Side-chain Heavy-atoms c” N, C”, c Side-chain Heavy-atoms C” N, c”, C Side-chain Heavy-atoms C N, c”, C Side-chain

t r.m.s. difference of atoms after least-squares $ Includes heavy-atoms 0, C8 and beyond.

fitting

of all heavy-atoms.

(A)?

460

C. B. Post et al

total energy to those in the kinetic energy remained less than 0.05, an acceptable value for stability (Brooks et al., 1983: van Gunsteren & Berendsen, 1977). The temperature. averaged over 0.5 ps. stayed within a 5 deg. K window of the simulation average temperature. There was. however, a very slow decrease in the total potential energy of the system amounting to about 1 kcal/mol in 3 ps. The results reported here were obtained from coordinate sets taken everv 0.1 ps. unless otherwise

indicated. In order t,o a&id motional effects due to overall net rotation or translation of the protein. an approximation to removing the rigid-body motion was made by rotating each dynamics co-ordinate set into the optimized X-ray co-ordinates by a least-squares fitting

procedure with mass weight,ing. The rotation angle determinedin this manner wasnever more than 6”.

3. Results ln the analysis of the molecular dynamics simulations of free and bound lysozyme, the deviations from the X-ray co-ordinates and the fluctuations in the atomic positions are examined. A comparison is made between the atomic positions and the C$and $ dihedral angle values of the timeaveraged dynamics structures and the X-ray structure. The distributions of side-chain dihedral angles, many of which undergo transitions during the simulations, are obtained from individual coordinate sets in the trajectory; these distributions are compared with distributions in dihedral angle values in the X-ray structure. The main-chain hydrogen-bonding patterns and the variations in the donor-acceptor distances and angles during the trajectory are also analyzed and related to the X-ray results. With respect to positional fluctuaCons, variations among residues and differences between the free and the bound simulations are described. The fluctuations computed for the native simulation are compared to the values estimated Debye-Waller factors. from crystallographic Finally. the dynamics of t,he st’ructural water molecules included in the simulations are examined. To facilitate the discussion of the relation between the dynamics results and the structural features of the lysozyme structure? we present a brief overview of the X-ray struct,ure (Figs I and 2). The enzyme is divided into two domains separated by the cleft that binds the substrate. One domain, consisting of residues 1 to 39 and 85 t,o 129, is composed of four a-helices plus a five-residue, oneturn 3,0 helix. The short stretches of backbone loops and turns connecting the helices are between five and nine residues in length. One a-helix, helix K (residues24 to 36), is totally buried, while the three helices A (residues 4 t,o 15). C (residues X8 to 99) and D (residues 108 to 115) lie on the protein surface. partially exposed to solvent; the 3,, helix (residues 119 to 124) is partially exposed to solvent. The second domain of lysozyme, consisting of residues 40 to 84, includes a three-stranded, antiparallel P-sheet (residues 42 to 60), a second small B-sheet (residues 1 to 2 and 39 t,o 40). and a single-turn 3,, helix (residues 79 t,o 84). Bet)ween

Figure 1. Schematic drawing of’ hen epp-white lysozpme(from Make rt al.. 1965). the large ,&sheet and the 3,, htalix IS 8 long c~oil~~tl loop region strrt,ching from residws 61 to TX. Residues lining the cleft include the first anal SNYIIK~ strand of the p-sheet (residues 13. 14. 16. 5:! and 56 to 59). the carboxyl terminus of helix I< (rrsidrw 35). the loop connecting helix (’ and helis I) (residues 98 a,nd 101 to 103). the amino terminus ot helix D (residues 107 to 110 and 112j, and wsidur~ 62, 63 and 73. (a) Pornparison

with th,r ?i-ray

structwe,q

(i) Radius of gyratiorr An overall measure by which to cwnpare similar structures is the radius of gyration:

where rCM is the posit’ion of t.he ~ent~rr of mass. ri is the position of atom i and rei its mass.The radius of gyration was determined for dynamics structures The averaged over fire-picxosrcond int.rrvals. average

of

these

five-picosecond

valurs

for

t 11~1

native, CM native and bound simulations are given in Table 3. There is no significant difference between the various calculated values and those for the two X-ray st)ruct,ures. which hare essentialI>

Table 3 Radius of gyration of the X-ray stmcture and the average !value for that of the Jive-picosPcorLd-avr~ag4d dynamics structures i
of gyration

(21)

Molecular Dynamics Simulatiow of Lysoxyme

461

(bi Figure 2. Stereo drawing of hen egg-white lysozyme backbone atoms N, C” C and 0 made with the HYDRA program (R. Hubbard, unpublished results). Three-center water molecules with the hydrogen atoms constructed by the CHARMM program (Brooks et al., 1983) and the substrate, (GlcNAc),, are shown. (a) Lysozyme and 53 water molecules. (b) Lysozyme, 35 water molecules and (GlcNAc),.

Figure 3. Stereo drawing of the lysozyme residues in the active site and the substrate (GlcNAc),. The sugar molecules in sites A. B and C are from the crystallographic co-ordinates, while sugar molecules in sites D, E and F are from model-building (see the text).

462

C. B. Post et al.

identical radii. This result for the radius of gyration indicates that the present choice of van der Waals’ radii is such as to prevent the protein contraction found in earlier vacuum simulations. In part, the larger van der Waals’ radii compensate for the neglect of solvent (van Gunsteren & Karplus, 1982). (ii) Differences in atomic co-ordinates A comparison of the full simulation timeaveraged structures with the corresponding X-ray structure was made after fitting the 1001 protein heavy-atom positions by a least-squares rotation (Table 2). The r.m.s. differences in the co-ordinates of the 1001 lysozyme heavy-atoms between the dynamics average structures and the X-ray structures are 2.2 A, 2.2 A and 2.0 A for the 101. picosecond native, the CM native and the bound simulations, respectively. The r.m.s. co-ordinat’e deviation for all heavy-atoms from the earlier lysozyme simulation is 2.8 A (Ichiye et al., 1986), significantly larger than that for the simulations reported here. For the C’ atoms, the r.m.s. deviations are 1.5 A, 1.5 w and 1.4 A for the three sets of structures obtained in the present simulations. Similar values are found for the r.m.s. deviation for all main-chain atoms N, C”, and C. Some idea of the significance of these deviations can be gained by noting that the values are close to those found when comparing the C” at’oms of homologous proteins such as hen egg-white and human lysozyme (0.7 8: Artymiuk & Blake, 1981), two bacterial species of dihydrofolate reductase (1.l 8; Bolin et aZ., 1983) or the reduced and oxidized forms of tuna cytochrome c (0.9 to 1.0 -4; Mandel et al., 1977). Simulations of other proteins report C” deviations of 2.2 A4 (van Gunsteren & Karplus, 1982) and 1.5 A for bovine pancreatic t’rypsin inhibitor (Levitt, 1983) and 1.5 A for ferrocytochrome c (Northrup et al.. 1980). As expected, the side-chain atoms have larger r.m.s. deviations than do the backbone atoms; the values of 2.5 .&, 2.6 A and 2.3 a were found for the three simulations (Table 2). For the most part, the side-chain residues with differences greater than 2.5 A were residues near the surface of the protein and, of these surface residues, a substantial portion were charged or polar. Side-chains for all arginine residues showed large differences, as did the carboxyl-terminal leucine. Co-ordinates for the (GlcNAc), atoms from the 55-picosecond-average bound structure were also compared with those from the minimized modelbuilt structure. When the lysozyme heavy-atoms are least-squares fitted, the deviation is 1.8 A for the (GlcNAc), heavy-atoms, intermediate between the main-chain and side-chain values. Fitting the (GlcNAc), heavy-atoms, the deviation is 1.3 A. similar to the size of the deviation found in the lysozyme backbone atoms. The co-ordinate deviations of the sugars in sites A, B and C (i.e. those obtained from the crystallographic results) are 2.1 A with fitting lysozyme heavy atoms and 1.1 with fitting only (GlcNAc), atoms.

Table 4 Root-mean-square main-chain co-ordinate differences of secondary structural elements between the X-ray and average dynamics structures r.m.s. Structure

Residues

P-Sheet

Nativef

72

14

1.2,39-60

cc-Helix A I3 (‘ I) 3,,

No. atoms

4-15

X--84

Loops

120-124

IX 15

16. 23 61 78 85- 87 100-107

24 54 9 24

1~ 129

Whole

(A)t

Ikmndg: 0.5

36 39 36 24

24.-36 88-99 108-115

helix

difference

387

The main-chain atoms S, C” and 0 are included. t r.m.s. difference of atoms after least-squares heavy-atoms. $ 101.ps-average structure afxwus X-ray structure. 5 55-ps-average structure versus X-ray structure.

fitting

of all

The discrepancies between t)he time-averaged dynamics and X-ray C” co-ordinates are not the same for different regions of the molecule. The a-helices and P-sheets have the smallest deviations (see Table 4). As shown in Figure 4 for both t,he native and bound simulations, the a-carbon atoms of the buried helix B, of residues 39 and 40, and of the two interior P-sheet strands (residues 49 to 54 and 57 to 60) differed in atomic positions by less than 0.8 A from the X-ray positions. The bound-

40

2 c a

3.0 2.0 I.0 0

1 (b)

I G c

a

#

?

II

3.0

2.0 I.0 0

20

40

60 Residue

SO

100

120

number

Figure 4. Magnitude of (I” co-ordinate deviations between the time-averaged dynamics structure and the X-ray structure. (a) lOl-ps average of the native simulation.

(b) 55-ps

average

of the bound

simulation.

Molecular

Dynamics

Simulations

trajectory average structure, in addition to these groups of atoms, showed similarly small deviations in the a-carbon atoms of helix C (residues 88 to 99) and the first strand of the b-sheet (residues 42 to 47). It is noteworthy that P-sheet residues 43, 44 and 46 make contact with (GlcNAc),, and it is suggested that the interactions between this exposed p-strand and the substrate stabilize the crystallographic conformation. The largest deviations when comparing the C” positions of either dynamics-average structure with the corresponding X-ray structure were in the large coiled-loop region from residue 66 to 79, as well as at the carboxyl terminus; the co-ordinate shifts and motions in this loop region are considered below (see Discussion). The short loop (85 to 87) connecting a 310 helix and the C helix has a considerably lower deviation in the free than in the bound trajectory; the origin of this difference is not evident, though it may be related to the inverse behavior of helix C (see Table 4). To examine the time development of the deviation of the dynamics co-ordinates from the X-ray results, the co-ordinate sets obtained from the trajectory were averaged over intervals of five picoseconds. The r.m.s. difference between the fivepicosecond average atomic positions and the crystallographic positions is plotted as a function of time in Figure 5. The r.m.s. differences remain nearly constant for both native trajectories over the entire period. It should be noted that the first fivepicosecond-average difference in Figure 5 is taken from the beginning of the analysis period, so that 17 picoseconds of dynamics corresponding to heating and equilibration of the structure have already taken place. This suggests that the molecule moves away from the X-ray structure fairly rapidly to a new equilibrium structure, and that the dynamics samples a region of co-ordinate space in the neighborhood of the equilibrium structure. Comparison of the five-picosecond co-ordinates amongst themselves yields r.m.s. deviations as large as 1.1 A for the set of all protein heavy-atoms. In the bound trajectory, the time-course of the coordinate differences shows a slight increase at the beginning of the trajectory, and appears to be leveling off near the end. This suggests that the enzyme-substrate structure had not been equilibrated completely at the beginning of the analysis, even though the behavior of total energy and of temperature was similar to those for the unbound trajectory. The bottom curve in Figure 5 shows the deviations in the (GlcNAc), co-ordinates for the same five-picosecond intervals; these are essentially constant with time. To interpret the structural differences between the free and bound trajectories, it is useful to compare

the crystal

structures.

The lysozyme

main-

chain conformation appears not to be altered extensively by the binding of (GlcNAc),. The r.m.s. difference in the C” positions from the two X-ray structures is 0.26 A; the very small difference may be in part due to the fact that the (GlcNAc),

463

of Ly.sozyme

2.50,

I

I

I

1

I.25 I 0

I 20

1 40 Srnulatlon

I 60 time (ps)

I SO

I 100

Figure 5. Time-dependence of the r.m.s. co-ordinate deviations between the 5ps dynamics time-averaged coordinates and the X-ray co-ordinates. The atoms included in the average are the 1001 lysozyme heavy-atoms from the (0) lOl-ps native simulation, (0) 27-ps CM native simulation, and ( x ) 55-ps bound simulation and the 85 (GlcNAc), heavy-atoms from the (0) 55-ps bound simulations.

refinement started with the free lysozyme coordinates. The C” atoms with the largest deviations are located at the two chain ends, the exposed, coiled-loop region (68 to 75), and the two tryptophan residues 62 and 63. In these regions, (GlcNAc), makes contact with residues 62, 63 and 73. From the simulations, structural differences between the free and the (GlcNAc),-bound dynamics-average structure are found in the 66 to 77 loop region, the exposed P-strand 44 to 47, the loop residues 101 to 102, and in the carboxylterminal residues 127 to 129. All but the carboxylterminal region make contact with (GlcNAc),. Main-chain dihedral angles The main-chain dihedral angles were analyzed in two ways. The individual 4 and $ angles for pairs of structures were differenced, and the distribution of (4,II/) pairs on a Ramachandran plot was examined. Differences in the 4 and Ic/ angles between the average dynamics structures and the X-ray structures were calculated. The r.m.s. deviation of the I#Jand Ic/angles from the X-ray values is 30” and 32”, respectively. This is of the order of the accuracy with which the dihedral angles are determined by protein crystallographic methods; i.e. for an uncertainty of 0.25 A in the atomic coordinates for the four dihedral atoms, the dihedral angles are known only to within approximately 30”. The larger deviations are less than 50” for most residues, as shown in Figure 6. This corresponds to the fact that the dihedral angles’tend to remain in the same minimum and preserve the character of (iii)

the

secondary

structure.

.There

are anticorrelated

changes in many of the dihedral pairs; i.e. the angle $ of residue i and 4 of residue i+ 1 rotate in opposite directions, so that the intervening peptide group rotates but there is only a small perturbation on the rest of the main-chain (McCammon et al., 1977). An example is provided by the peptide group

464

C. B. Post et al.

-1201

i I

I I

I

I

/

I

1

I

I

/

20

40

60

/ 80

Residue

/ 100

1 120

number

Figure 6. Differences for the C#J and I) angles of each residue. The differences were taken between the 101.ps codynamic average co-ordinates and the X-ray ordinates.

are located in t>he same regioli of angles conformational space for the two structures, in accordance with the analysis given above. The largest changes are found for glycine residues, particularly those with values of 4 greater than zero. That larger deviations are found for glycine residues is not surprising because of the greater rotational freedom present in t,he main-chain dihedral angles due to the lack of a side-chain. The latter also makes it more difficult to determine accurate values for the dihedral angles in the X-ray structure. Several of t,he residues with values of 4 greater than zero in Figure 7 are not glycine residues. These non-glycyl residues! which occur in a high-energy region of the Ramachandran plot, are listed in Table 5. Some of the larger values seen in Figure 6 for A$ or A+ come from such non-glycine residues in the X-ray structure. They include AsnlS, ArgPl

Bound

dynamic

Bound

X-ray

500

involving $21 and 422, one of the larger deviations in Figure 6. For Gly71, Asn74 and Leu75, the large changes in c$,$ values are of interest and may be related to the existence in this part of the main of more than one conformation, as chain determined by crystallographic studies (Berthou et al., 1983). The main-chain transitions for the peptide group between residues 74 and 75 result in the loss of a reverse-turn hydrogen-bond interaction between 74 and 77. -4 Ramachandran plot of the 4,II/ angles for each residue from the lOl-picosecond dynamics average structure and the X-ray structure is shown in Figure 7. It can be seen that, in general, the $,$

t

I300

c

i

6

J-l Free

I_-

nr

IBC

l-l

dynamic

9c

630 7 : D

I n c

5

Free X-ray 8

-9c

q 6

X

L

-I

I EC

q X

I

I

-90

0

+

I

90

Figure 7. Ramachandran plot of the $,$ pairs. Cp,$ values obtained from the X-ray structure (+, x ) and from the 101 ps dynamic average structure (0.0). Glvrine reiidues

residues are by + and 0.

I

1,

(deg.)

indicated

by

x

and

0:

ot)hrr

-180

_I

J

il I

I

0

90

XI

(a)

Fig. 8.

l-l

Molecular

Dynamics

Table 5 Main-chain dihedral angle values for the X-ray and native trajectory average structures for the nonglycine residues that map in the glycine region of the Ramachandran plot in the X-ray structure X-ray

Dynamic

Residue

4

ti

AsnlS Arg21

71 59 63 66 43 51 53

18 28 34 18 60 48 52

Asn37

Phe38 (k.57 Asn74 Asn77

Bound

f#J 50 50 67 70 71 -82 68

average Ic, 70 -45 32 71 88 118 38

Simulations

of Lysozyme

and Asn74, of which only Asn74 moves out of the “glycine” region during the simulation. Arg21 adopts a conformation close to a C-‘lax configuration with (50”,-44”), as compared with the X-ray values of (59”,28”). The newly formed hydrogen bond between the main-chain Try20 oxygen CO and Gly22 NH is only a weak bond, and does not break the reverse-turn hydrogen bond between Tyr20 CO and Tyr23 NH. A possible cause for the structural change in the main chain at Arg21 is the hydrogenbonding of its side-chain. The guanidinium group of Arg21 is not hydrogen-bonded in the X-ray structure, but during the equilibration period becomes bonded to the carbonyl group of Ser99, displacing a water molecule from Ser99. This water

I

dynamic

465

Bound

t

Bound

X-ray

Free

dynamic

n

Bound

Free

240

Free

dynomlc

250

X-ray

dynamic

rlh

-

Free

75

X-roy

X-ray

t

0

90

180

X2 (b)

Figure 8. Side-chain dihedral angle distributions. (a) x1 an&es with unbranched tetrahedral b-carbons (Asx, Glx, Arg, Lys. Leu): (b) x2 angles with tetrahedral y-carbons (Arg, Glx, Lys, Met); (c) x2 angles with trigonal y-carbons (Asx). The X-ray distributions are constructed by including all the residues of the given type. The dynamics distributions are cxonstructed from individual co-ordinate sets at I-ps intervals for all of the residues of a given type. In each case, the ordinate gives the number of residues of a given type in the X-ray structure or summed over the dynamics co-ordinate sets.

466

C. B. Post et al.

molecule, having multiple hydrogen bonds, moved to bind more strongly with Gly104. As such, the absence of solvent in the simulation may not be compensated fully by the few included water molecules, resulting in the co-ordinate shifts near this exposed loop region involving residues 20 to 23. It is possible, as well, that the dynamics results represent an alternative solvent-mediated structure. Averaging of two or more conformations can be responsible for unusual angle values and some of the differences seen in Figures 6 and 7. A possible example is provided by Arg5 at (3”, -59”) in the dynamics structure, relative to (-49”, -61”) in the X-ray structure. Substantial dynamics averaging has occurred for C$ of Arg5. In the trajectory, 4 undergoes transitions between -50” and + 50”, giving an average close to o”, which represents a very improbable value for C$ in any individual coordinate set. One of the t,wo conformations, that with 4 = -50”, is close to the X-ray results. In some cases where multiple sites are occupied, the X-ray refinement procedure gives a stereochemically reasonable model by picking out one of the possible sites or some average stru&ure (Kuriyan et al., 1985). The dynamics suggests that Arg5 has two conformations corresponding to C$ of -50” and +50”. Examination of the X-ray map at low contour levels does not give any indication of an alternative conformation for Arg5 at +50” and, if such a conformation exists, it must have a low occupancy. Another example of averaging in the dynamics is Gly102. In the simulation, the C$angle of Gly102 jumps between values of +90” and - 120”, and the calculated average is - 140”. This suggests that the reason for the discrepancy between the dynamics average value of - 140” and the X-ray value of 150” could be incomplete sampling in the lOO-picosecond simulation; i.e. since the dihedral angle convention yields angles between 0” and 180” and between 0” and - 180”, - 140” and 150” correspond to somewhat different averages of 90” and -120”. (iv) Side-chain dihedral angles The individual dynamic co-ordinate sets at onepicosecond intervals are used to examine the sidechain dihedral angles. This is necessary because, unlike the main-chain dihedral angles, the sidechains undergo many transitions, making average angle values meaningless. The dihedral angle distributions for residue types obtained from the full trajectory are compared with the initial lysozyme X-ray structure distribution. The distributions from the lysozyme X-ray structure, after crystallographic refinement, are similar to other empirical distributions tabulated from many crystallographic structures (Janin et al., 1978). Histograms are made for particular sets of dihedral angles based on the distance alone: the chain (x1yxZ,x3) and on the type of atoms”forming the bond about which dihedral rotation occurs (branched or unbranched carbons and trigonal,

tetrahedral or aromatic carbon atoms). Examples of the distributions obtained are shown in Figure 8 and include x1 (Fig. 8(a)) for unbranched b-carbon atoms (Asx, Glx, Arg, Lys and Leu) and x2 (Fig. S(b)) for tetrahedral y-carbon atoms (,4rg, Glx. Lys and Met). The prominent, peaks in the X-ray distributions are preserved (Janin et al., 1978; Gelin & Karplus, 1979). but with the expected motional broadening (Hoch et al., 1985). Examination of other dihedral angle classes show that t)he X-ray side-chain conformational distributions are generally maintained during the simulation, even though many of the individual residues undergo dihedral angle transitions. An exception to t’he agreement cited occurs for the x2 angles of bhe trigonal II-carbon in asparagine and aspartate. We consider these results in some detail because they illustrate questions t)hat arise in any comparison between the ?i-ray average structure and dynamics simulation results. The distributions for the 21 Asx residues are shown in Figure 8(c). Empirical tabulations and the distribution from t’he lysozyme X-ray structure for this dihedral angle type exhibit little preference for an) particular orientation about x2. However, the distribution obtained from the molecular, dynamics co-ordinat,e sets shows a definit,e peak at 90”, the staggered conformation of t,he planar carboxylate or amide group. The result)s for Glx-xJ are similar, t,hough the small number of five Glx residues makes the X-ray distribution less significant, and the peak of t,he dynamics dist,ributions at 90” is less pronounced. It is important) to note that the alt#ered dihedral dist,ributions in thr dynamics generally did not eliminate the hydrogen~ bonding of the Asx side-chains found in the X-ray structure. Most, of the interactions were maintained, although some hydrogen-bond pairs with 061 and 062 (Asp) or N62 (Asn) were lost and others wclrf’ formed. Several fa.ct,ors could lead tcr thr discrepancy in the x2 distributions. Of t.he nine Asx residues (18. 37. 48, 74, 77, 87. 93. 113 and 119) with x2 values that, are eclipsed in t)he X-IX> structure but, that have x2 distributions closer t>o 90” in the simulat.ion. five residues (IX. 37, 77. X7 and 119) have Debye-Wailer factors greater than 50 for at> least, one of the three carboxyl (COO) or amide (CON) heavy atoms. A large Debye-Wailer factor is consistent wit,h a highly mobile functional group, commonly found for surfacae residues, a,nd can be a refleckon of non-uniqueness in t,hta xZ angle. Janin et al. (1978) noted that Asx side-chain positions are of%en poorly determined due to t htb weak electron density for these residues and that a large uncertainty in the dihedral angle result,s. Asn93, which has a low temperature factor, forms hydrogen bonds in the X-ray no side-chain struct,ure. Thus, rotation of xZ is expected t’o bfb relat,ively free, as found in the dynamics. The possibility of multiple sites for Asn93 is not. contradicted by its low temperature fact,ors (Kuriyan et aE., 1986). Intermolecular contacts in the crystal can be the reason for the discrepancy

Molecular

Dynamics

between the X-ray result and the dynamics distribution for two of the nine residues, Asn37 and Asnl13 (Gelin & Karplus, 1979). These two residues are in crystal-contact regions. Finally, xZ values of residues 48 and 74 are changed in the simulation because these residues are in a part of the protein where the main chain is quite flexible. Large mainchain fluctuations in the simulation and large Debye-Waller factors are found at the p-bend, 48 to 50, and large co-ordinate differences occur in the main chain of the loop 73 to 79. That the molecular dynamics simulation, with only a few structural water molecules, shows shifts in surface residue conformation relative to a crystal environment is to be expected (Gelin & Karplus, 1979). Potential energy maps for xZ rotation in the X-ray structure were constructed, holding all atoms except the two polar atoms of either carboxylate (Asp) or amide (Asn) group fixed. The resulting rigid rotation maps, which are dominated by nonbonded contacts, have minima at the xZ values found in the X-ray structure for most of the residues. When analogous rotation maps are made for structures from the dynamics trajectory, it is found that all but five (27, 46, 66, 87 and 93) of the 21 Asx residues have altered xZ potential surfaces. These changes are either a diminution in the rotational barrier to a nearly flat surface or 30” to 50” shifts toward 90” in the position of the minimum. Surfaces of only two (87 and 93) of the nine Asx residues with altered xZ values in the simulation were unchanged in the dynamics structure. The X-ray values for these two angles do not fall at potential minima and move to minimumenergy values in the simulation. The potential surface for the dihedral angle x2 is determined by a balance of hydrogen-bonding interactions and non-bonded contacts with neighboring atoms, as well as by the usual dihedral torsion potential (Gelin & Karplus, 1975). Small structural changes occurring during the simulation remove the interactions responsible for the less favorable eclipsed x2 values found in the X-ray structure, and allow xZ to adopt the usual staggered conformation. None the less, certain eclipsed x2 angles, involving close 1,4 interactions between C” and the N6 or 06 atom of Asx, are preserved. Two such cases are Asn39 and Asp66, whose x2 value is maintained near zero in the simulation. Asn39 strongly in the eclipsed hydrogen-bonds conformation and Asp66, a buried residue, is held in the eclipsed configuration by the surrounding protein matrix. Two Gln residues (57 and 121) have eclipsed x3 values in the simulations. Although side-chains with large positional shifts in the dynamics are typically solvent-exposed residues, the buried Ile98 side-chain deviates by an amount substantially larger than the r.m.s. residue average for the native trajectory. The displacement of Ile98 is due to the motions of xZ, which undergoes several transitions. The behavior of Ile98 is interesting, in that the rotation of xZ is anticorrelated with x1 of Lys96, both in the dihedral

Simulations

467

of Lysozyme

-180 0

20

40 Slmulatton

60 time (ps)

80

100

Figure 9. Time series for x1 of Lys96 and xZ of Ile98 value and in the fluctuation magnitude. As shown in Figure 9, the dihedrals exhibit two states; from 0 to 63 picoseconds, x196 has small oscillations near -75” and x298 has very large oscillations between the gauche and trans wells. In the second state, xi96 has large oscillations between the gauche-minus and tram wells, while x298 has small oscillations near 160”. The dihedral angles observed in the X-ray structure correspond to the values for which x1 and x2 have the small oscillations (i.e. -75” for x196 and 160” for x298). Apparently, only the more rigid conformation is observed in the X-ray refinement, while the alternative conformation of higher flexibility, and therefore less well-defined density, is not. The dynamics of residues where side-chain transitions occur and the effect of the motions on X-ray structure analyses are being examined (Kuriyan et al., 1986). (v) Main-chain hydrogen bonding Hydrogen-bonding interactions are one of the dominant features stabilizing different secondary structural elements in proteins. There are many energetic factors, including steric constraints, packing of side-chains, electrostatics and solvation, which contribute to the stability-of the CO-HN interaction, in addition to the hydrogen-bond potential (eqn (1)). To analyze the interactions involved, a comparison of the CO-HN geometry from the X-ray co-ordinates with that from the dynamics co-ordinates was made. Specifically, the heavy-atom O-N distances, the donor (O-H-N) angles and the acceptor (C-O-H) angles are examined. Both the full set of main-chain-mainchain hydrogen bonds and the subsets occurring in the various secondary structural elements are discussed. Of particular interest is the distribution of the 8C-O-n angle, since crystallographic results for small molecules (Taylor & Kennard, 1984) and for proteins (Artymiuk & Blake, 1981; Glusker & Murray-Rust, 1984; Baker & Hubbard, 1984) show &-O-H to be frequently in the neighborhood of 120”; e.g. the lysozyme X-ray structure puts the acceptor angle near 120” for many reverse turns (hydrogen-bonding between the peptide carbonyl of residue i and the peptide amide of residue i+ 3).

468

C. B. Post

et al.

------/

180

360

I 60

120

16

I6

30

12

12

8

8

4

4

90

120 N-H-O

150 angle

140\(b)

180

I-

‘7 20 I 1

90

(deg.)

120 C-O-H

Dynamic

‘$

150 angle

180

2

3

(deg 1

4

5

R,,,ix)

60

80

90

120 N-H-O

150 angle

180

90

(deg:)

I20 C-O-H

150 angle

180 (deg.1

r-----120

Dynamic

(c)

--

60 45

80

30 40 15 :I X-ray

90

120 N-H-O

1

150 180 angle (deg.)

90

120 C-O-H

150 angle

180 (dep.~

240 c

60

60

X-ray

81

/

IO 8

6

6

4

4

2

2

l----A 90

120 N-H-

150 180 0 angle (deg.)

90

120 150 C-O-H angle (deg.)

2

3

4

5

Ro,N( fi )

Figure 10. Hydrogen-bond geometry distributions for CCHS. The distributions from the compared with the distributions from 1 -ps-averaged dynamics struct’ures for the native simulation. type are indicated on the ordinate. Only hydrogen bonds with energies < -0.2 kral/mot (eyn (1)) distribution. Geometric parameters are the heavy-atom distance 0-h’. the hydrogen angle O-HPK. angle (lFOPH for main-chain-main-chain hydrogen bonds from (a) all residues: (b) r-helix residues

X-ray

strucaturr

WV

the numbers of eacah were included in the and t,hr ac~ceptot 1 to 15. 21 to 36. 88

Molecular

Dynamics

Simulations

This raises the question of whether the angular dependence on 8C-O-H in equation (1) may be too strong. &uantum mechanical calculations (Reiher & Karplus, unpublished results) indicate that the acceptor angle well is broader than that given by equation (1) but does not have a minimum of 120” in the lone-pair directions. For hydrogen-bonded CO-HN pairs, distributions in the O-N separation, the O-H-N angle and the C-O-H angle were determined from the native X-ray structure and from the one-picosecondaverage co-ordinate sets for the native simulation. All possible pairs &thin a 6.0 A O-N distance and with angles greater than 100” were examined and a pair was included in the distribution if the hydrogen-bond energy was more negative than - 0.2 kcal/mol. Distributions for all main-chainmain-chain hydrogen bonds calculated from the X-ray structure and from the simulation structures are shown in Figure 10(a); the results have not been divided by sin 0. In the dynamics distributions, there is a small tendency for the angles to become more linear than in the X-ray structure, particularly for the acceptor angle. The radial separation of 0 and N does not change significantly. Distributions were calculated also for subsets of atoms grouped by secondary structure. Main-chainmain-chain hydrogen bonds were evaluated from residue groups corresponding to a-helical, P-sheet’ or 3,0 reverse-turn structures determined from the crystallographic co-ordinates; the specific groups are listed in the legend to Figure 10. For the helical group, all main-chain hydrogen bonds of energy < -0.2 kcal/mol were included, without restriction t’o i and i+4 pairs. All a-helix type and P-sheet hydrogen-bond pairs in the X-ray structures were present during the dynamics simulation, except for one a-type hydrogen bond lost between residues 111 and 115 at the carboxyl end of helix D. The peptide hydrogen-bond distributions for the four a-helices and /?-sheet in lysozyme determined from the dynamics structures (Fig. 10(b) and (c)) are very similar to those determined from the X-ray coordinates. The distributions for hydrogen bonds of reverse turns (Fig. 10(d)) show larger discrepancies between the X-ray and the simulation values. There is i peak at 120” in the acceptor angle C-O-H distribution from the X-ray structure, while in the dynamics the distribution is close to being level between 120” and 150”. The donor angle O-H-N is also more linear in the dynamics distribution. Furthermore, there is more dynamic broadening in the N-O distance distribution for the reverse-turn hydrogen bonds than may be seen for that in a-helices or in P-sheets. During the simulation, several reverse-turn hydrogen bonds were broken, including pairs 13-16,

of Lysozyme

469

19-22, 36-39, 54-57, 69-72, 74-77, 81-84, 98-101, 104-107, 115-118 and 119-122. The 54 to 57 turn has an X-ray C-O-H angle of 111” and is a buried turn in the P-sheet. In the simulation, the local structure remains close to the X-ray structure, even though there is no hydrogen-bonding interaction. Two a-helical hydrogen bonds (not included in the distributions for Fig. 10(b)) were formed in the simulation; these were between pairs 12-16 and 80-84, which in the X-ray structure were reverseturn hydrogen bonds 13-16 and 81-84. The contrary situation, a reverse-turn hydrogen bond becoming a stronger interaction at the expense of an a-type bond, also occurred in the simulation; at the amino-terminal of the A helix, the hydrogen bond energy between residues 4 and 7 was more favorable, while that between residues 4 and 8 was weaker. Factors involved in the structural changes of the reverse turns are described below (see Discussion). (b) Fluctuations in atomic positions In this section we consider the magnitude and time development of the atomic fluctuations that occur during the simulations and compare their dependence on residue position with magnitudes calculated from the temperature factors determined in the X-ray refinement. (i) Time development The r.m.s. fluctuations calculated for 20picosecond intervals and averaged over all heavyatoms of lysozyme were 0.55 A, 0.63 a and 0.61 a for the native, CM native and bound trajectories, respectively. A previous simulation for lysozyme (Tchiye et al.. 1986) yielded a value of 0.76 a for a 30-picosecond interval. The present values for lysozyme, shown in Table 6, are somewhat smaller than the Wpicosecond fluctuations reported for a simulation of pancreatic trypsin inhibitor (0.7 A; van Gunsteren $ Karplus, 1982: Swaminathan et al., 1982) and the 20-picosecond fluctuations of cytochrome c (0.85 8; Morgan et al., 1983). The molecular dynamics simulations of all three proteins were calculated with similar potential functions. with similar potential functions. Positional fluctuations were calculated for timeintervals of varying lengths and then averaged for these intervals over the trajectory (Swaminathan et al., 1982). The time-development of the fluctuations illustrates the rate at which an atom explores the available configuration space. Intervals of only a few picoseconds do not allow a complete sampling of the space, and the r.m.s. fluctuations for short intervals are small. As the time-averaging period increases, larger fluctuations are observed. Values for the native and bound simulation fluctuations

to 100 and 108 to 115: (c) /?-sheet residues 1 and 40 to 60; (d) reverse-turns13 to 16, 17 to 20, 19 to 22. 20 to 23, 24 to “7, 36 to 39, 39 to 42, 46 to 49, 69 to 72, 74 to 77, 79 to 82, 80 to 83, 81 to 84, 98 to 101, 103 to 106, 104 to 107. 105 to 108. 108 to 11 I. I15 to 118. 119 to 122, 120 to 123, 121 to 124 and 124 to 127. (Other reverse-turn interactions in the X-ray structure with R,,. < 4.0 A but with either angle < 100” and therefore not included in the distribution are 12 to 15. .54 to 57 a.nd 94 to 97.)

C.

B. Post et al.

Table 6 Average root-mean-square Jluctuations for varying time intervals r.m.s. fluctuation8 (A) Atom type

20 ps

Native trajectory 1001 lysozyme heavy-atoms Lysozyme N, c”, C Lysozyme 0, CB and beyond Water oxygen atoms

0.55 0.45 0.62 0.69

CM native trajectory 1001 lysozyme heavy-atoms Lysozyme N, c”, C Lysozyme 0, CB and beyond Water oxygen atoms

0.63 0.52 0.71 1.20

Bound trajectory 1001 lysozyme heavy-atoms Lysozyme N, c”, C Lysozyme 0, Cr and beyond Water oxygen atoms (GlcNAc), heavy-atoms

0.61 0.50 0.69 0.92 0.65

(0.63)1 (0.51) (0.70) (0.71)

50 Ps

100 ps

0.61 0.50 0.67 0.74

0.65 0.53 0.72 0.78

0.68 0.56 0.76 1.09 0.65

t Values in parentheses are calculated from a single 20 picosecond period of the native trajectory for comparison with t,hc CM trajectory.

calculated from intervals of three different lengths of time are shown in Figure 11. The r.m.s. fluctuations were avera ed over main-chain (N,C”,C) and side chain (0,C P, and beyond) atoms of each residue. (The peptide oxygen is included in the side-chain averaging because its motional properties are

similar to those of CD.) The motions of main-chain atoms are less than those of side-chain atoms. Residue-specific motions are apparent in the fluctuations but only after the time interval over which the averaging is carried out becomes sufficiently long. The lowest curves in each panel of Figure 11 correspond to one-picosecond fluctuations and show only a slight variation with residue number (Swaminathan et al., 1982). As the averaging time increases and larger amplitude motions contribute (upper curves), a pattern emerges, reflecting the residue-specific mobilities. Examination of the main-chain fluctuations from the entire lOl-picosecond native and the 55 picosecond bound trajectories (top curves in Fig. 11(a) and (c)) shows that there is a general correspondence between fluctuation magnitudes and secondary structure; helical and P-sheet regions are more rigid than coiled-loop regions. The ol-helices and b-sheet residues, marked along the bottom of Figure 1I (a) and (c), have the smallest fluctuations. particularly helix R (residues 24 to 36), the internal helix, and the /?-sheet. In contrast, the largest main-chain atom fluctuations are from the b-bend residues 47 to 49 and the coiled-loop regions 16 t,o 23, 67 to 75, 83 to 87 (Fig. 11(c) only), 102 to 104 (Fig. 1l(a) only), and 116 to 120. The time-development of the atomic displacement,s also varies with the type of residue (Post, et al., unpublished results); i.e. the spacing between the curves in Figure 11 is not uniform. In the case I.._

I-50

__~1-_.-7..--

-----..~

~_ , ,

!C) I

G

0

20

40

60 Residue

00 number

100

120

.25 c

0

,>55 ps

I

20

40

60 Residue

80

$00

.i

120

number

Figure 11. Residue average of the positional fluctuations of (a) main-chain atoms (N, C”, C) and (b) side-chain atoms (0, Cfl and beyond) for the 101.ps native simulation, and (c) main-chain atoms and (d) side-chain atoms for the 5%ps bound simulation. Averaging times over which the fluctuations were calculated are indicated for each curve. Residue ranges marked in (a) and (c) indicate (m) a-helix, (ES=) /?-sheet or (-) 3,” helix.

Molecular

Dynamics

Simulations

of the lOl-picosecond native simulation (Fig. 11(a) and (b)) the main-chain fluctuations of a substantial number of residues, mostly in the b-sheet and a-helices, have nearly converged in that the values for 50-picosecond averaging almost equal those for 101.picosecond averaging. For other residues, particularly with respect to the side-chain, there is a substantial increase in the lOI-picosecond fluctuations over the 50-picosecond fluctuations, reflecting continued sampling of configuration space. The peaks of Figure 11 generally have the largest spacing between the two top curves, indicating that large fluctuations often have long timedevelopments. This trend, however, does not hold for all residues with high mobility; fluctuations of the /?-bend residues 47 to 49 in the bound trajectory are large yet of equal value for 20-picosecond and 50-picosecond averaging. The large amplitude motion at this P-bend is therefore of sufficiently high frequency that it converges during the simulation. The extent of convergence in the fluctuation time-development, as well as that of statistical errors arising from the use of single trajectories, can be examined by comparing the 50-picosecond fluctuations calculated from the two halves of the lOl-picosecond trajectory (Fig. 12). The fluctuations of main-chain atoms show only small changes

471

of Lysozyme

in magnitudes and, in some cases, these are oneresidue shifts in the relative values (e.g. residues 100 to 105). The side-chains show greater variation between the two halves of the trajectory. (ii) Comparison of the native and bound trajectories In comparing the native and bound simulations, one must recognize that some disparities in fluctuations exist because of the limited sampling of configuration space in each simulation, such that not all differences are due to the bound substrate. The expected statistical error may be estimated from the comparison in Figure 12 (see above). It indicates that small changes in magnitudes or shifts in peaks of a few residues cannot be interpreted as resulting from the presence of (GlcNAc),. Fluctuations from the lOl-picosecond native and the 55-picosecond bound trajectories are shown in Figure 13. The values are the main-chain and sidechain residue average r.m.s. fluctuations over 50 picoseconds. For the bound trajectory, this corresponds to the first 50 of 55 picoseconds and, for the native simulation, the first and second 50picosecond fluctuation values are averaged. Two sets of differences in the residue fluctuations between the two trajectories are significantly

(b)

1

60 Residue

SO

100

120

number

Figure 12. Residue-averages of the r.m.s. fluctuation from the 2 halves of the lOl-ps native simulation: (a) mainchain atoms (N, c”, C) and (b) side-chain atoms (0, C@and beyond). The average r.m.s. fluctuations for 0 to 50 ps (---) are very similar to those for 50 to 100 ps (. . . .).

C. B. Post et al.

. . . .

(b)

. .

. t

V”

Residue

Y”

iL”

I””

number

Figure 13. Comparison of residue-averages of t’he 50-ps r.m.s. fluctuations from the native simulation ( ) and from the bound simulation (.. . . .). (a) Main-chain atom (N. C”. C) and (b) side-chain atoms (0. CB and beyond). r.m.s. fluctuations from the two 50-ps intervals from the native simulation (see Fig. Id) were averaged. Two regions with substantial differences between (a) and (b) are residues 67 to 88 and residues 101 to 105.

greater than those found in Figure 12. The chain from Asp101 to Met105, the loop connecting helix C and helix D, shows quite large amplitude motion in the native trajectory, but in the bound trajectory this part of the protein is rather rigid. The difference is almost certainly due to substrate binding, since AsplOl, Gly102 and AsnlO3 form hydrogen bonds with t’he N-acetylglucosamine residues occupying sites A and B in the binding cleft. Thus, these interactions between the enzyme and substrate reduce the large motions that this otherwise flexible loop undergoes in the free enzyme. Not all large deviations in the fluctuations of residue main-chain atoms with and without (GlcNAc), present are correlated with proximity to the substrate. This is demonstrated by the second region (residues 67 to 88), where there are substantial differences between the two trajectories. Unlike the region between 101 and 105, this stretch of residues is mostly on the side of lysozyme away from the binding site. The region includes the coiled-loop from Gly67 to Ile78, the 3,,-helix turn (Pro79 to Ile84), and the connecting loop to helix C (Ser85 to Ile88). Motions in this part of the

molecule are larger in the bound trajectory than in the free trajectory. The reason for the different, motional behavior is not readily a,pparent,, since in this part of the chain, only the side-chain of Arg73 makes substantial contact with (GlcXAc),. Tt is of interest that a study of the hinge bending modes of free lysozyme (Brooks & Karplus, 1985) has shown that some residues far from the cleft in the back of the molecule are involved in the motion. (iii) Comparisons X-vay

of molecular temperature factors

dynamics

nvvd

of lyst,zyme. In R crystallographic study individual temperature factors for all heavy-atoms were determined (Handoll et rcl.. unpublished results). This makes possible a comparison of the fluctuations calculated from the temperature factors with those from molecular dynamics. If contributions from lattice disorder or experimental uncertainty are assumed to be negligible (Sternberg Qt al.. 1979; Petsko & Ringe, 1984), isotropic, harmonic atomic displacements may be estimated from crystallographic temperat,ure factors according to the expression (Willis & Pryor, 1975:

Molecular Karplus

& McCammon,

Dynamics

Simulations

1983):

473

of Lysozyme

than the actual positional fluctuation. The direction of the error depends on the nature of the atomicposition distribution (Kuriyan et al., 1986; Ichiye & Karplus, unpublished results). Crystallographic refinement, including individual Debye-Waller factors, is being carried out on the inhibitor (GlcNAc),-lysozyme complex (Artymiuk et al., unpublished results). A direct comparison between the simulation results and the crystallographic results cannot be made, of course, since the simulation involves the substrate (GlcNAc),lysozyme complex. It is nevertheless interesting to note that the largest differences in the main-chain Debye-Wailer factors relative to native lysozyme are increases for residues 70 and 71 (Artymiuk et al., unpublished results), two residues in the region with fluctuations in the substrate-bound simulation larger than those in the free simulation.

B = &r=(AR;>/3.

(3) Residue r.m.s. displacements from the native simulation are compared in Figure 14 with those calculated from temperature factors according to equation (3). The agreement in the main-chain r.m.s. fluctuations between the X-ray and dynamics results is quite good. Somewhat larger discrepancies are seen when comparing the two curves for the side-chain atoms. The extent of the agreement in Figure 14 is similar to that found between temperature factors of homologous proteins (e.g. tetragonal and triclinic lysozyme; Artymiuk & Blake, 1981). Some of the side-chains with large differences are surface residues whose motions are not well modeled by the simulation, due to the lack of solvent (e.g. residues 16, 20, 71 and 72). Multiple occupancy of side-chains can also lead to discrepancies in the fluctuations calculated from the simulations and the temperature factors. If an atomic-position distribution has more than one high-density peak, crystallographic refinement assuming single-site occupancy can result in temperature factors that are either smaller or larger

(c) Dynamics

of the bound

water molecules

In this section we compare first the results for the two native enzyme simulatBions with different constraints on the water molecules. We then consider in some detail the interactions between the

60 Residue

number

Figure 14. The residuevariation in the r.m.s. positional fluctuations calculated from crystallographic temperature factors (upper curve, to scale)and from the lOl-ps native dynamic simulation (lower curve, displaced-0.25 A). Residue averagesof the r.m.s. fluctuations were taken over (a) main-chain atoms (N, C”, C) and (b) side-chainatoms (0, Cfi and beyond). secondary structures indicated asin Fig. 11.

C. B. Post et al. water molecules and lysozyme. cross-correlation calculations dynamic coupling.

and make use of to analyze the

(i) Positional changes and jhetuations The two types of constraints on water molecules, used to compensate for the absence of a full solvent environment, are described in Methods. On average, water molecules constrained to their initial position showed lower mobility and smaller co-ordinate deviations from the X-ray structure than did those constrained to the center of mass of the protein. The 20-picosecond r.m.s. fluctuation in the ST2 oxygen position, averaged over all 53 water molecules, was 0.69 A (average value of the 5 20picosecond intervals) for the lOl-picosecond trajectory, substantially less than the fluctuations of 1.20 L% for the CM trajectory. Correspondingly, the r.m.s. positional shift for all water oxygen atoms from the 20-picosecond-averaged structure was 1.8 A in the lOl-picosecond trajectory and 3.5 A for the CM trajectory. This is as expected, since the initial positional constraint does not permit the water molecules to move away from the X-ray result, while they are free to do so with the center-of-mass constraint. Such diffusional or hopping motion of the water molecules is indicated also by the increase in the average co-ordinate deviation over the course of the CM trajectory. As shown in Table 7, the co-ordinate deviation for the five-picosecond-averaged structures steadily increases, in contrast to constant five-picosecond values for the lOl-picosecond trajectory, of which the first four are given in Table 7. With regard to the overall effect on the protein from the two water-constraining potentials, no differences were found in the co-ordinate deviations of lysozyme (in Table 2 compare native lysozyme versus X-ray with CM native versus X-ray) or in the r.m.s. fluctuations (in Table 6 compare 20-ps

Table 7 Comparison of fluctuations and co-ordinate deviations for two water moleculeconstrain&g potentials Native? Water

oxygen

atoms

20 ps fluctuations Co-ordinate 20.~9

deviations11 average wwus

5-ps average o-5 ps 6-11 ps 11-15 ps 16-20 ps

w-su.9

X-ray

CM

nativeS

(4

(4

0.75

1.2

1.8

3.5

1.8 1.8 1.8 I.8

3.1 3.4 3.7 4.0

X-ray

t Water molecules constrained to initial optimized X-ray ordinates by a quartic potential. $ Water molecules constrained to the protein center-of-mass co-ordinate by quartic potential. 5 Average of r.m.s. fluctuations from 5 20-pa intervals. 11 Time-averaged dynamics structure r.m.s. co-ordinate deviation from the X-ray structure.

co-

fluctuation values for the native trajectory with fluctuations for CM native trajectory). The small differences in the 20-picosecond fluctuation values in Table 6 appear to involve statistical sampling; e.g. the fluctuations for one of the 5 20-picosecond periods in the native trajectory, shown in parentheses, are nearly identical with the CM native fluctuations. (ii) STSlysozyme interactions The water molecules that interact strongly with a protein heavy-atom were examined by comparing two dynamic properties. First, the magnitude of the 20-picosecond fluctuations of the water oxygen atom was compared with that of the protein heavyatom with which it interacts and, second, the CPOSScorrelations in the fluctuations of the water oxygen atom and the protein heavy-atom were examined for the strongest water-protein interactions that are common to both trajectories. To determine the strength of the interaction between a water molecule and a lysozyme heavyatom? the energy of a protein-water hydrogen bond was obtained from equation (1). Although an explicit hydrogen-bond potential between the protein and the water molecules was not included in the simulation (see Methods), equation (1) can be used to examine the relative strengths of proteinwater interactions. The time-dependence of the hydrogen-bond energies for all protein-water pairs was calculated from one-picosecond-averaged coordinates. As expected, a variety of patterns in the time series of hydrogen-bond pairs was found, including strong interactions that persisted throughout a simulation, short-lived but strong interactions and weak interactions. Examples of time series of hydrogen-bond strengths of lysozyme-water pairs are shown in Figure 15. The time average of the hydrogen-bond energy was determined for each lysozyme-water hydrogen bond, and the lysozymr heavy-atom that is most strongly hydrogen-bonded was found for each water molecule. Of the 53 water molecules included in the simulation, 40 interacted most strongly in both trajectories with the same lysozyme heavy-atom and these water-lysozyme pairs were used for the calculation of the water-protein fluctuation crosscorrelations. Nine water molecules included in the simulation did not hydrogen-bond to the protein. but instead hydrogen-bonded to other water molecules. How the fluctuations of a water oxygen atom correlate with those of the most strongly hydrogenbonded lysozyme heavy-atom was analyzed hy plotting the 20-picosecond fluctuations of a water oxygen atom against those of t)he protein heavy atom (Fig. 16); results for both the positionally constrained and the CM constrained simulation are given. A similar plot has been reported for r.m.x. displacements estimated from temperature-factors for water molecules and protein in the crystallographic studies of tortoise and human lysozyme (Blake et al.: 1983). The plot in Figure 16(b) for the

Molecular

Dynamics

Simulations

475

of Lysozyme

I

0,

I

I

I

I

I 2.0

I 2.5

(b) 1,5-

LO-

-+otb-0

I

,

20

40

b

I

I

60 Time

80

1

100

molecule hydrogen bonds. The top panel shows persistent hydrogen bonds of water molecule 8 with (------) NH-41, (----) NHT-1 and (......) 07-39. The bottom panel shows transient hydrogen bonds of water molecule 7 with (---) NH,-114, (------) H-l-water molecule 10 and (......) H-2-water molecule 10.

0

I 0.5

I I.0

I I.5 ST2 r.m.s.

3.0

(8)

Figure 16. Correlation of the 20-ps r.m.s. fluctuations of the water molecules with that of the lysozyme heavyinteracts for (a) the center-of-mass constrained trajectory and (b) the initial position constrained trajectory. The strength of interaction is determined by averaging the hydrogen-bond energy obtained from eqn (1) using 1 -psaveraged co-ordinates of either the lOl-ps native or the CM native simulation.

is given

lOl-picosecond simulation, with initial

positional constraints, is similar to that obtained by Blake et al. (1983); that is, the values cluster around the equivalence diagonal, with a somewhat larger fraction falling to the right of the line in the simulation than in the X-ray results. In the CM trajectory, the water fluctuations tend to be larger than the protein fluctuations, so that points fall further below the equivalence diagonal than is found for the native simulation with initial positional constraints or the crystal results. Of the points in Figure 16(a) for water fluctuations equal to 1.5 A or greater, eight interact with either lysine or arginine side-chains, which themselves are highly mobile. Some of these side-chains have larger fluctuations in the CM than in the lOl-picosecond trajectory. The cross-correlations in atomic fluctuations of lysozyme heavy-atoms with the water oxygen atoms were calculated for the strongest waterlysozyme pairs common to both trajectories, as stated above. We wished to determine whether the motions of a water molecule and the protein atom to which it is hydrogen-bonded are correlated, and potential

gc.:

atom with which the water molecule most strongly

(ps)

Figure 15. Time series of the energy of certain water

whether the constraining the correlation. The

0.5-

has an effect on

cross-correlation between displacements fr_om the mean positions of lysozyme heavy-atoms (ARL) and water oxygen atoms (AR,)

by:

(AR,. AR,> (AR;) ‘/=(AR;> I/=’

(4)

The value of equation (4) varies between - 1.0 and 1.0, where an absolute value of 1.0 results from perfectly correlated motion of atom L with 0, and a value of 0.0 is the result if the motions are uncorrelated. Values for the fluctuation cross-correlations between the water oxygen atom and the protein heavy-atom with which the water molecule interacts most strongly are given in Table 8. For both constraining potentials, the cross-correlations are large. (Only about 10% of the cross-correlations between water oxygen atoms and any lysozyme heavy-atom are greater than Oe25.)In many cases, the cross-correlation behavior is similar in the two trajectories, although there are also significant differences; e.g. N of residue 5 has a positive water correlation in the native trajectory and a negative correlation in the CM trajectory. On the average, the cross-correlation values are greater for the CM native trajectory (average cross-correlation equals 0.47) than for the native trajectory (average crosscorrelation equals 0.36). The larger cross-correlation values indicate that the less-confining CM constraint allows a water molecule to move somewhat more freely with the protein.

C.

B. Post et al.

Table 8 Atomic jluctuation cross-correlations between a water oxygen atom and the lysozyme heavy-atom to which it is most strongly hydrogen-bonded (‘rowcorrelation Lywzpmc residue

Atom

0 N Cl NH1 XHl SH% SD2 ODl OEl ?1’ N s s OF, 1 NE2 s 0 SE1 N NH2 N (I(: NH2 OG N N N NZ NZ N NH1 SH:! SHl SE NH1 XH2 NZ N 0 N

LVater residue

Native

I 0.26 040 04% 040 0~51

0.2

2 4

3x 26 28

26 13

x 12 16 PO 45 27 30 1x 35 43 43 48 53 49 .50 1 I7 31 33 3” 3X 42 46 34 41 7 9 9 II 22 21

I0 ;5

-

CM native

0.39 0.27 041 0.09 0.09 0.56 0.43 0.30 0.30 0.35 0.6 I 0.34 0.49 0.1 1 0.44 0.S 0.31 0.33 0.14 0.49 041 0.12 0.33 0.65 040 0~20 0.28 0.16 0.18 0.59 0.59 0.36 0.32 0.57

0.37 - 0.26 046 0.84 0.60 0.80 0.66 0.33 0.45 0.16 0.32 0.65 041 040 0.25 WI 0.55 0.38 0.56 0.04 0.57 0.30 0.52 0.10 0.34 0.61 0.69 0.10 0.55 0.77 0.45 0.32 W-17 0.64 0.63 0.37 0.75 0.81 0.32 0.57

0.36

047

Values for the initial position constrained trajectory (Native) and the center-of-mass constrained traject,ory (CM native) are shown. Only pairs common to both trajectories are listed.

4. Discussion In the current’ molecular dynamics study of the motions of both native lysozyme and substratebound lysozyme, several improvements have been made over previous simulations (Northrup et al., 1980; van Gunsteren & Karplus, 1982: Tchiye et aZ.. 1986). Larger van der Waals’ radii, a new hydrogenbond potential involving explicit polar hydrogen atoms, and careful computational procedures for heating, equilibration and updating pair-interact,ion lists resulted in long, stable trajectories of over 1500 atoms that required no external adjustment,s to conserve energy. Crystallographic water molecules found to be twice hydrogen-bonded to the prot,ein (53 for the native simulation and 35 for the bound simulation) were included because of their

structural importance inferred from the strong hydrogen-bonding. These water molecwlns sa,tisf> in part the hydrogen-bonding rrquiremrnt,s for I~soz)m~ side-chains and. most importani.1~. occupy sites in the ac%ive-site cleft in the frrtl enzyme a,nd are displaced 1)~ the substrate in the bound cnzymr: their presence in tht> free enzyme permits a more meaningful comparison of the dynamicas of free and subst,rate-bound lysozymr. .1 longer simulation t i tne (100 ps) provides dynamical information about lower-freyuencv motions and a more rigorous test of the simulat;on method t)han did earlier st udirs. The average dynamics structures are close to the X-ray structures in most) regions (r.m.s. differeucc of 2.2 A for all protein heavy-at,oms and 1.5 A for main-chain atoms). Moreover. the st,ability of thtb simulations \va\2 achieved without artiticially adjusting the paratnet,ers. such as forc(l con&ants or well dept,hn. for the particular prot,ein undrr caonsideration (Levit’t. 1983: set also Aqvist of ~1,. 1985). The energies of interactions are standarcl. few-body potential terms and t hr paramet,ers originate from small-molecule studies (Brooks c~f~1.. 1983). (‘onstraining potentials. however. wele applied to the simulated structura,l wa,ter molecules t)o prevent tjhem from becoming separated from the prot,ein (Levitt. 1980). Comparison of the prot,ein behavior using a constraining potential with reference either to initial oxygen positions or to t’he protein c*ent)er of mass showed no differences in t’he overall protein fluctuations or (to-ordinate deviations from the S-ray strnc%ure. The wa.ter oxygen ~iom tluctuat ions a4 co-ordinat,r deviations. on the ot,hrr hand. n-err larger overall in the case of the center of mass c&onstraint and the c+ordinate dtlviation of the wat’er molecbules increased throughout the simulat ion. \Vit,h regard to water-lysoe~rnr interactions. the

initial

positional

constraint

gives

a

cmr’elatiorl

the wat,er oxygen atjon fluctuation and thr fluc+uation of’the Iysozymr heavy-atom with whitoh it most strongly hydrogen-bonds t)hat is more consistent with the crystallographic caorrrlat)ion of I< factors (Hlakr et a./.. 1983). On this basis. it is cGonc4uded that the initial positional c:onst,raint is pt~f~ferrt4, Thtb sequence of residues from 66 t)o 79. a long. exposed coiled-loop region. exhibits larger t,han average fluctuations and devia,tions in atomic, positions in the simulation relative to the crystal structures. (‘ertain peptide groups along t)his part, of the chain are quite flexible-in bot,h the k~ountl anti free trajectories: the peptide links Gly67 -Ser68. Arg73-Asn74 and Asn77FArg78 have large anticaorrelated oscillations in the $ of residue i and 4 of residue i+ 1 Common changes in the main-chain dihedral angles art’ found in the two simulations WK will: shifts between the tirne-averaged structurt, and thr S-ray st~ructures occurred for $66, $68. $70; 471, $74 and t)htb pair $73--474. The st rnctural c*hangt‘s associated wif,h the altered dihedral angles arv localized in thr cxpos~d IOO~J. between

Molecular

Dynamics

Simulations

Residues 80 to 83 and residues with numbers below 66, which lead back to the interior of the protein, retain positions close to those in the X-ray structure. Crystallographic data are consistent with the polypeptide chain in this loop being rather flexible and capable of adopting several conformations; more than one conformation for residues 70 to 75 has been observed crystallographically for lysozyme in different crystals (Berthou et al., 1983 and references therein; Ford et al., 1974; Kelly et al., 1979). Furthermore, some of the largest differences in the main chain between free lysozyme and inhihit’or complexes are in residues 68 to 75, as reported for the trisaccharides (GlcNAc), (see Result’s) and MurNAc-GlcNAc-MurNAc (Kelly et al., ‘19791, and for the monosaccharide GlcNAc (Perkins et al.? 1978)) again demonstrating the rather flexible nature of this loop region. The X-ray geometry of the main-chain-mainchain CO-HN hydrogen-bonding interactions is generally reproduced during the simulation, as shown by similar distributions in the N-O distances and the N-H-O and C-O-H angles. For reverset’urn hydrogen bonds, however, there was a tendency for the CO-HN geometry to become more linear in the simulation and, in some cases, the interaction was completely lost, reverse-turn suggesting that the form of the hydrogen-bond potential (eqn (1)) may have too strong an angular dependence on 0C-O-H. As already mentioned, detailed quantum mechanical calculations of the hydrogen-bonding interaction have shown that the C-O-H acceptor angle corresponds to a broad minimum (Reiher & Karplus, unpublished results). None the less, it is unlikely that the source of the discrepancy lies solely in the angular dependence of the hydrogen-bond potential. We note that one reverse-turn hydrogen bond lost in the simulation, 115-l 18, in fact has comparatively linear donor and acceptor angles (37” and 36”, respectively) in the X-ray structure. In view of the complexity of a protein molecule, other influences may be involved in maintaining the integrity of the reverse-turn st’ruct’ure. For instance, it is thought that the role of solvation is of major importance in stabilizing this type of interaction (Rose et al., 1983) for which there are unfavorable energy factors both in hydrogen bonding (main-chain hydrogen bonding of the intervening residues cannot be satisfied by other protein donors or acceptors) and in the dihedral angles of the intervening residues. The value for in one type of reverse turn is near O”, as is true $i+2 for all but one of the reverse turns lost in the simulation. The change seen in the dynamics average structure results in either an increase or a decrease in Ii/i + 2 away from the potential barrier at 0”. This shift in + values can be seen in Figure 7 as the decreasein the number of points near $ = 0” for the dynamics structure. Thus, the loss of most of the reverse-turn interactions appears to be associated with a move to a region of lower energy on the dihedral-angle potential surface. Differences between residues in both the

of Lysozyme

477

magnitude and the time-dependence of the positional fluctuations are found to be related to structural features; e.g. helices and /I-sheet structures generally have more restricted mobilities than do coiled-loop regions. The main-chain atoms of the buried helix B have the smallest fluctuations, and the fluctuation magnitudes reach a plateau value for intervals of 50 picoseconds and greater. Both the stabilizing forces of the main-chain hydrogen bonds and the confinement) of the surrounding protein atoms contribute to the rigidity of helix B. Other helices have larger fluctuations and take longer to reach their limiting values. The coiled loops near the protein surface show the largest fluctuations, and the fluctuation magnitudes increase up to the full lOl-picosecond period of the simulation because of either long timescale events, such as concerted motion of a large section of the chain, or rare events such as dihedral rotations. For these regions, therefore, longer simulations are required to obtain full convergence. A detailed examination of the relation between secondary structure, and fluctuation magnitude and its time-development shows some exceptions to the general correspondence just described (Post et al., unpublished results). A large group of residues (20 to 60 in Fig. 11(a)), which includes both a-helix, P-sheet and coiled loop, has main-chain fluctuations that plateau after 50 picoseconds, indicating that only higher frequency motions are sampled by this part of the protein even though all types of protein structure are present. It is noteworthy that residues 20 to 60 form a large fraction of the binding cleft for sites C, D and E, including the position at which hydrolysis occurs. The substrate alters the motions of lysozyme in the simulation such that both larger and smaller atomic fluctuations are observed in different parts of the enzyme. That the structural fluctuations in lysozyme are perturbed by the binding of inhibitors has been demonstrated by proton exchange (Cassels et al., 1978). The exchange rates of the indole protons of Trp62, Trp63, Trp108 and Trplll are slowed in the presence of GlcNAc. Trp62, Trp63 and TrplO8 are close to GlcNAc (Perkins et al., 1978), but Trplll is in the hydrophobic box and more distant from the binding site. Two regions of the enzyme differed in their fluctuation magnitudes between the native and bound simulations by an amount significantly greater than statistical errors, as estimated by the differences between the two halves of the native simulation. One of these regions, Asp101 to Met105 makes several contacts with the substrate in sites A and B and shows a substantial decrease in mobility in the presence of (GlcNAc),. Unlike the more rigid lower cleft region of residues 20 to 40 (see above), the upper cleft is quite flexible in the free enzyme state. The second region, Gly67 to Ile88, which is for the most part distant from the active site and makes only a single contact with the substrate at Arg73, has larger fluctuations in the presence of (GlcNAc),. The altered motional behavior of residues 67 to 88 could

478

C. B. Post et al.

be in responseto structural adjustments around 73, whose side-chain moves down to hydrogen-bond with the substrate. Another possible cause for changes in the motional behavior of residues remote from t,he active site may be an indirect response to the presence of (GlcNAc), in the cleft. The starting structure for the bound trajectory was based on the X-ray co-ordinates for the trisaccharide (GlcNAc), Consequently, motions leading to complex. structural adjustments that result in more favorable interactions in the hexamer (GlcNAc), complex, might occur in the simulation. The larger fluctuations distant from the active site would be in response to substrate contacts, of which the effects are transmitted through more rigid structures such as the P-sheet or helix C to the less constrained loop region 67 to 88. A recent analysis of the hingebending mode in lysozyme has shown that, some residues distant from the active-site cleft are involved in the motion (Brooks & Karplus, 1985). The fluctuation differences between the native and bound trajectories for residues Asp101 to Met105 are interesting in view of the observed binding energies of the inhibitor (GlcNAc), (Imoto et aE., 1972), which binds mainly in sites A, B and C. The side-chein of Asp101 hydrogen-bonds to the inhibitor and protonation of the carboxylate should lower the enthalpy of binding. However, a decrease in pH from 5.3 to 2, which neutralizes the carboxyl group of AsplOl, results in only - 1 kcal/mol reduction in the observed binding free energy due to compensating enthalpic and entropic factors (Banerjee & Rupley, 1973). The enthalpy of binding at pH 2 is decreased by -6 kcal/mol, consistent with the loss of hydrogen-bonding by AsplOl, while the entropy of binding becomes more favorable by approximately 15 cal/mol per K. These data can be interpreted qualitatively in terms of the molecular dynamics results. We propose that the entropic effect, as well as the enthalpic effect, involves the interaction of the bound molecule with lysozyme residues 101 to 105. The protonation of Asp101 weakens the hydrogen-bonding strength, reducing not only a favorable enthalpic factor but also an unfavorable entropic factor. Comparison of the free and bound simulations shows that the fluctuations of residues 101 to 105 are very large in the free enzyme but are substantially reduced by the presence of (GlcNAc),. The hydrogen bonds of AsplOl, Gly102 and AsnlO3 with (GlcNAc), anchor this part of the enzyme, thereby decreasing the entropy in the bound state; this decrease is expected to be reduced when Asp101 is protonated. The balance of the two effects leads to the small net, change in the binding free energy observed in the experiment.

5. Conclusion Molecular dynamics simulations of heri egg-white Iysozyme in the free state and with a bound substrate have been analyzed. The variation in

fluctuations of residues in secondary structural elements has demonstrated a relationship between the structure and dynamics of this protein. Such structure-dynamics characterizations provide insight into the nature of protein motions that may be useful for interpreting a variety of phenomena, such as substrate binding and hydrogen exchange. Comparison of the simulation results with the structure and fluctuations obtained from X-rag crystallography helps to establish the degree of success of the simulation in representing the behavior of the protein in a, crystal. Differences betwc:n the average dynamics and crystal structure point to regions in the protein where multiple conformations may occur. Evaluation of the nat,ure of the interactions important in maintaining the integrity of the protein may lead to improved energy functions for simulations and can be useful in X-ray structure analysis. Molecular dynamics provides the opportunity to study systems, such as enzyme-substrate complexes, that are not readily accessible to most experimental techniques for determining structure or dynamic properties. Comparison of t,he free and bound simulations has suggested that specific residues contribute to t,he entropy of binding. Moreover, some of the disparities in the dynamics of the two structures suggest novel effects of substrate binding. Although t,hese features could result from artifacts of the simulations, it is hoped that the) will be explored by future experiments. Changes in structure and mobilities in different regions of the protein found in the two simulations can be examined by techniques, such as X-ray diffraction. nuclear magnetic resonance and site-directed mutagenesis.

We thank W. Reiher for advice concerning the hydrogen-bond potential, and C. I,. Brooks III, T. Ichiye and W. Pulford for helpful discussions. C.M.D. and D.C.P. are members of the Oxford Enzyme Group, which is supported by the Science and Engineering Research Council. The work at Harvard was supported by the National Institutes of Health. Support is also acknowledged from NIH (to C.B.P.), EMBO (to P.,J.A.). and the MRC (to P.J.A. and J.C.C.). The dynamics were performed on the CRAY IS at Daresbury. supported 1)) the Science and Engineering Research Council.

References Aqvist, J.. van Gunsteren, W. F., Leijonmarck, M. & Tapia, 0. (1985). J. Mol. Viol. 183, 461-478. Artymiuk, P. J. & Blake, C. C. F. (1981). .J. Mol. B&Z. 152, 737-762. Baker, E. N. & Hubbard, R. E. (1984). t’rog. Biophys. Mol. Biol. 44, 95-179. Banerjee, S. K. & Rupley, J. A. (1973). J. Bid. (lhenz. 248, 2177-2124. Banerjee, S. K., Holer, E., Hess, G. I’. & Rupley. J. A. (1975). J. Biol. Chem. 250, 4355-4367. Berthou, J., Lifchitz. A.. Artymiuk, P. t Jolles, P. (1983). Proc. Roy. Sot. ser. B, 217. 471-489.

Molecular

Dynamics Simulations

Blake, C. C. F., Koenig, D. F., Mair, G. A., North, A. C. T., Phillips, D. C. & Sarma, V. R. (1965). Nature (London), 206, 757-761. Blake, C. C. F., Johnson, L. N., Mair, G. A., North, A. C. ‘I’., Phillips, D. C. & Sarma, V. R. (1967). Proc. Roy. Sot. ser. B, 167, 378-388. Blake, C. C. F., Pulford, W. C. A. & Artymiuk, P. J. (1983). J. Mol. Biol. 167, 693-723. Bolin, J. T., Filman, D. J., Matthews, D. A., Hamlin, R. C. & Kraut, J. (1983). J. Bio2. Chem. 257, 1365013662. Brooks, B. R. & Karplus, M. (1985). Proc. Nat. Acad. Sci., 1J.S.A. 82, 8458-8462. Brooks, B. R., Bruccoleri, R. E., Olafson, B. D., States, D. J., Swaminathan, S. & Karplus, M. (1983). J. Comput. Chem. 4, 187-217. Cassels, R., Dobson, C. M., Poulsen, F. M. & Williams, R. J. P. (1978). Eur. J. B&hem. 92, 81-97. Delepierre. M., Dobson, C. M., Howarth, M. A. $ Poulsen, F. M. (1984). Eur. J. B&hem. 145, 389395. Ford, L. 0. Johnson, L. N., Machin, P. A., Phillips, D. C. & Tjian, R. (1974). J. Mol. BioE. 88, 349-371. Gelin. B. R. & Karplus, M. (1975). Proc. Nat. Acad. Sci., U.S.A. 72, 2002-2006. Gelin. B. R. t Karplus, M. (1979). Biochemistry, 18, 1256-1268. Glusker, J. & Murray-Rust, P. (1984). J. Amer. Chem. Sot. 106, 1018-1025. Grace, D. E. P. (1980). Ph.D. thesis, Oxford University. Hendrickson, W. A. (1985). Methods EnzymoE. 115, 252270. Hoch, J. C., Dobson, C. M. & Karplus, M. (1985). Biochemistry, 24, 3831-3841. Hodsdon, J., Seiker, L. C. & Jensen, L. H. (1974). Amer. Crystallogr. Assoc. Abstr. 2, ser. 2, 78. Hogle, J., Rao, S. T., Mallikarjunan, M., Beddell, C., McMullans, R. K. & Sundaralingham, M. (1981). Acta. Crystallogr. sect. B, 37, 591-597. Ichiye, T. & Karplus, M. (1983). Biochemistry, 22, 28842893. Ichiye, T.? Swaminathan, S., Olafson, B. 8: Karplus, M. (1986). Biopolymers, in the press. Imoto, T., Johnson, L. N., North, A. T. C., Phillips, D. C. & Rupley, J. A. (1972). The Enzymes, 7, 665-868. Janin, J., Wodak, S., Levitt, M. & Maigret, B. (1978). J. Mol. Biol. 125, 357-386. Karplus, M. & McCammon, J. A. (1983). Annu. Rev. Biochem. 53, 263-300. Kelly, J. A., Sielecki, A. R., Sykes, B. D. 6 James, M. N. G. (1979). Nature (London), 282, 875-878. Kurachi, K., Sicker, L. C. & Jensen, L. H. (1974). J. Mol. Biol. 101. 11-24. Kuriyan. J., Petsko, G., Levy, R. M. & Karplus, M. (1986). J. Mol. Biol. in the press. Levitt, M. (1980). In Protein FoEding (Jaenicke, R., ed.), pp. 17-39, Elsevier/North-Holland, Amsterdam. Levitt, M. (1983). J. Mol. Biol. 168, 595-620. Mandel, N.. Mandel, G., Trus, B. L., Rosenberg, J.,

Edited

of Lysozyme

479

Carlson, G. & Dickerson, R. E. (1977). J. Biol. Chem. 252, 4619-4636. McCammon, A., Gelin, B. R. & Karplus, M. (1977). Nature (London), 267, 585-590. Morgan, J., McCammon, A. & Northrup, S. (1983). Biopolymers, 22, 1579-1593. Moult, J., Yonath, A., Traub, W., Smilansky, A., Prodjorny, A., Rabonovich, D. & Saya, A. (1976). J. Mol. Biol. 100, 179-195. Northrup, S. H., Pear, M. R., McCammon, J. A., Karplus, M. & Takano, T. (1980). Nature (London), 287, 659-660. Olejniczak, E. T., Poulsen, F. M. & Dobson, C. M. (1981). J. Amer. C%m. Sot. 103, 6574-6580. Olejniczak, E. T., Poulsen, F. M. & Dobson, C. M. (1984a). J. Magn. Reson. 59, 518-523. Olejniczak, E. T., Dobson, C. M., Karplus, M. & Levy, R. M. (1984b). J. Amer. Chem. Sot. 106, 1923-1930. Perkins, S. J., Johnson, L. N., Machin, P. ,4. & Phillips, D. C. (1978). B&hem. J. 173, 607-616. Petsko, G. & Ringe, D. (1984). Annu. Rev. Biophys. Bioeng. 13, 331-371. Phillips, D. C. (1966). Sci. Amer. 215, 78-90. Pincus, M. R. & Scheraga, H. A. (1981). Biochemistry, 20, 3960-3965. Post, C. B. & Karplus, M. (1986). ,I. Amer. Chem. Sot. 108, 1317-1319. Poulsen, F. M., Hoch, J. C. & Dobson. C. M. (1980). Biochemistry, 19, 2597-2601. Rose, G. D., Young, W. B. 6 Gierasch. L. M. (1983). Nature (London), 304, 654-657. Rossky, R. J. & Karplus, M. (1979). J. Amer. Chem. Sot. 101, 1913-1937. Rossky, P. J., Karplus, M. & Rahman, A. (1979). Biopolymers, 18, 825-854. Ryckaert, J. P., Ciccotti, G. & Berendsen, H. J. C. (1977). J. Comput. Phys. 23, 327-341. Schindler, M., Assaf, Y., Sharon, N. & Chipman, D. M. (1977). Biochemistry, 16, 423-431. States, D. J. (1984). Ph.D. thesis, Harvard University. Sternberg, M. J. E., Grace, D. E. P. & Phillips, D. C. (1979). J. Mol. Biol. 130, 231-251. Stillinger, F. H. & Rahman, A. (1974). J. Chem. Phys. 60, 1545-1557. Sturtevant, J. M. (1977). Proc. Nat. Acad. Sci., U.S.A. 74, 2236-2240. Swaminathan, S., Ichiye, T., van Gunsteren, W. & Karplus, M. (1982). Biochemistry, 21, 5230-5241. Taylor, R. & Kennard, 0. (1984). Act. Chem. Res. 17, 320-326. van Gunsteren, W. F. & Berendsen, H. J. C. (1977). Mol. Phys. 34, 1311-1327. van Gunsteren, W. F. & Karplus, M. (1980). J. Comput. Chem. 1, 266-274. van Gunsteren, W. F. & Karplus, M. (1982). Biochemistry, 21, 2259-2274. Willis, B. T. M. & Pryor, A. W. (1975). Therm& Vibrations in Crystallography, Cambridge University Press, Cambridge.

by A. KEug