Protein Engineering vol.5 no.7 pp.703-714, 1992
A molecular dynamics simulation of bacteriophage T4 lysozyme
Gregory E.Arnold and Rick L.Ornstein' Molecular Science Research Center, Pacific Northwest Laboratory, PO Box 999, K2-18, Richland, WA 99352, USA 'To whom correspondence should be addressed
Introduction Experimental and theoretical studies have contributed to an increase in the understanding of the relationships between structure, function and dynamics in proteins and macromolecules in general. Molecular dynamics simulations provide a tool to examine both the time scale and magnitudes of atomic fluctuations. It has been suggested that fast thermal atomic fluctuations propagate slower motions and that these generate yet slower motions (Ansari et al., 1985; Dommair and Jahnig, 1989). Taken together, these architectonic levels of motion are thought
703
Downloaded from http://peds.oxfordjournals.org/ at University of California, San Francisco on December 5, 2014
An analysis of a 400 ps molecular dynamics simulation of the 164 amino acid enzyme T4 lysozyme is presented. The simulation was carried out with all hydrogen atoms modeled explicitly, the inclusion of all 152 crystallographic waters and at a temperature of 300 K. Temporal analysis of the trajectory versus energy, hydrogen bond stability, r.m.s. deviation from the starting crystal structure and radius of gyration, demonstrates that the simulation was both stable and representative of the average experimental structure. Average structural properties were calculated from the enzyme trajectory and compared with the crystal structure. The mean value of the Ca displacements of the average simulated structure from the X-ray structure was 1.1 ± 0.1 A; differences of the backbone $ and * angles between the average simulated structure and the crystal structure were also examined. Thermal-B factors were calculated from the simulation for heavy and backbone atoms and both were in good agreement with experimental values. Relationships between protein secondary structure elements and internal motions were studied by examining the positional fluctuations of individual helix, sheet and turn structures. The structural integrity in the secondary structure units was preserved throughout the simulation; however, the A helix did show some unusually high atomic fluctuations. The largest backbone atom r.m.s. fluctuations were found in nonsecondary structure regions; similar results were observed for r.m.s. fluctuations of non-secondary structure $ and ¥ angles. In general, the calculated values of r.m.s. fluctuations were quite small for the secondary structure elements. In contrast, surface loops and turns exhibited much larger values, being able to sample larger regions of conformational space. The Ca difference distance matrix and superpositioning analyses comparing the X-ray structure with the average dynamics structure suggest that a 'hinge-bending' motion occurs between the N- and C-terminal domains. Key words: computer simulation/protein dynamics/protein motions/T4 lysozyme
to give rise to lower frequency oscillations within the protein (Ansari et al., 1985). Hence, fast co-operative atomic fluctuations lead to side chain fluctuations, secondary structure fluctuations, backbone motions, domain movements and ultimately conformational changes. Thus, the collective dynamics of these regions are relevant to biological activity for proteins involved in catalysis and transport (Mao et al., 1982; Levitt, 1983; Ansari et al., 1985). Protein tertiary and secondary structures are stabilized by weak non-bonded interactions and, as a result, considerable atomic fluctuations should occur. The study and characterization of these fluctuations can lend insight into the nature of internal protein motions. X-ray crystallography and/or NMR studies primarily provide time average data but are unable to provide detailed information regarding the motional properties of a macromolecule. Assuming that the temporal behavior of proteins can be accurately modeled, molecular dynamics simulation methods can be used to investigate the collective oscillations within the molecule. Molecular dynamics simulation can provide information regarding both the equilibrium properties and the temporal characteristics of a particular molecule. We have used the method of molecular dynamics simulation to study the motional properties of bacteriophage T4 lysozyme (T4L). Other investigators have reported simulations for hen egg white lysozyme (Mao et al., 1982) but no trajectory has been reported and analyzed for T4L to date. Both enzymes are endoacetylmuramidases and both cleave the /3(1—4) glycosidic bond between A'-acetylmuramic acid and A'-acetylglucosamine, although T4L is ~ 250 times more active towards the enzymatic cleavage of Escherichia coli cell walls (Matthews et al., 1981). The two enzymes show little homology yet are thought to be evolutionarily divergent (Matthews et al., 1981). The vast amount of experimental data on the structure and function of T4L and therelativelysmall size of this enzyme make it an excellent system to study by simulation. T4L is a 164 residue, monomeric ~18.7kDa protein with no disulfide bridges. This enzyme has a globular bi-lobed structure with its active site cleft juxtaposed between the two lobes. The enzyme aids in the digestion of the bacterial cell wall facilitating release of the virus from the host and is produced late in the infection of E.coli by bacteriophage T4 (Travers, 1970). The wild-type structure of T4L has been refined to 1.7 A resolution, with R equal to 0.193, and the inclusion of isotropic thermal parameters (Weaver and Matthews, 1987). This structure serves as the starting point for our analysis. In the native enzyme structure, the catalytic cleft is occluded and thus obstructs substrate access to the active site. Consequently, conformational flexibility and dynamic fluctuations are necessary for substrate binding and enzymatic catalysis to occur (Faber and Matthews, 1990). A particularly important objective of our molecular dynamics study is to investigate the relationships of positional fluctuation relative to die movements of essentially 'rigid' structural units within the protein. In particular, we want to examine whether me methods of molecular dynamics simulation can model relative domain movements in
G.E.Arnold and R.L.Ornstein
T4L. Intra-domain movements have been reported for native T4L (Weaver and Matthews, 1987) and relative domain movements have been observed in five different crystalline allomorphs of the T4L M6I mutant (Faber and Matthews, 1990). In this paper the results of a 400 ps molecular dynamics simulation of T4L are presented. The analyses in this initial report focus on: (i) the global properties of the dynamics simulation; (ii) the relationship between atomic fluctuations and protein structural units; and (iii) the comparison of the X-ray structure to the time average dynamics structure. Of particular interest was the relative movement observed in the two lobes of the enzyme when comparing the time averaged dynamics structure with the X-ray structure.
The X-ray structure with added hydrogens was energyminimized using the method of steepest descents for 500 steps, with the positions of all heavy atoms fixed to remove any artifacts induced by the addition of explicit hydrogens. The structure was further minimized for 500 steps with only the positions of the heavy atoms of the protein fixed to allow the waters of crystallization to relax. Then, the structure was minimized for 5000 steps using the steepest descents method followed by 5000 704
Results and discussion To facilitate discussion of T4L, a brief description of the enzyme's structure is given below. T4L is a monomeric 164 amino acid, 18.7 kDa protein and is composed of two globular domains separated by an active site cleft (Figure 1) (Faber and Matthews, 1990). The N-terminal domain consists of residues 12-59 and comprises one relatively surface-exposed a-helix, helix B (39-50) and a four-stranded antiparallel /3-sheet. The first three strands of the sheet are contiguous along the sequence and are connected by /3-turns; the fourth strand is separated from the
Fig. 1. Ribbon trace of the X-ray crystal structure of T4L. Amino acid number* important to the discussion are labeled
Downloaded from http://peds.oxfordjournals.org/ at University of California, San Francisco on December 5, 2014
Materials and methods The coordinates of the reported X-ray structure of T4L at 1.7 A resolution (Weaver and Matthews, 1987), obtained from the Brookhaven Protein Data Bank (Bernstein et al., 1977), served as the initial molecular model for the molecular dynamics calculations. The system simulated comprised 164 amino acid residues and 152 crystallographic waters. All of the hydrogen atoms were explicitly modeled, bringing the total number of atoms in the system to 3099. The coordinates of the added hydrogens were generated according to idealized bond lengths and valence angles using the MOLEDT software package from Biosym Technologies (San Diego, CA). The molecular dynamics simulation was performed using the Discover simulation software package (version 2.41) from Biosym Technologies on a Cray X-MP. The default Discover model was used for the explicit (crystallographic) waters. No constraints were applied to the crystallographic waters and thus they were free to move from their crystallographically determined sites. No cross-terms were used in the energy expression and a simple harmonic potential was chosen for the bond stretching terms. All calculations were conducted with a group-based non-bonded cut-off of 10.5 A imposed over a switching distance of 1.5 A and the non-bonded pair list was updated every 20 time steps. A linear distance-dependent dielectric (equal to the interatomic separation) was used. The parameters used were those of the consistent valence force field (VFF) (Hagler, 1985; Dauber-Osguthorpe et al., 1988) except that the charges of acidic, basic, N-termini and C-termini 'functional groups' were made net neutral; although atoms were assigned net charges similar to their normally uncharged polar counterparts. This was done to compensate for the missing dielectric effect of bulk water, because only the crystallographic waters were included. This type of charge screening effect has been successfully used by others (Aqvist et al., 1985; Kruger et al., 1985; Makinen et al., 1989; Mehler and Solmajer, 1991; Paulsen et al., 1991; Solmajer and Mehler, 1991; Swaminathan et al., 1991; Braatz et al., 1992; Gu and Brady, 1992). In previous simulations of crambin, using net neutral residues with various dielectric models, a distance-dependent dielectric or a dielectric of 1 yielded nearly identical results and were superior to other dielectric models (Ornstein, 1990).
steps using the conjugate gradients method before beginning molecular dynamics. Dynamics was performed using the leapfrog algorithm with a 1 fs time step. A constant temperature for the simulation was maintained by weak coupling to a thermal bath (Berendsen et al., 1984). For the dynamics, 0.5 ps was first performed at 50 K. This was followed by elevating the temperature to 300 K, using an exponential approach (time constant of 2.0 ps), such that the target temperature was achieved after 10 ps. At 10 ps the time constant was adjusted to 0.1 ps and the dynamics was re-initialized at 300 K. The trajectory was continued for a total of 400 ps. The calculation was performed on a Cray X-MP computer and 1 c.p.u. h was required per 4.5 ps. The simulation was restarted at 100 and 200 ps with the previous structure used as the current starting structure and with new random velocities. The data set considered in this work comprised 801 structures saved at 0.5 ps intervals.
Molecular dynamics of T4 lysozyme
Time course of the simulation Several parameters have been examined that assess the stability and the convergence characteristics of the molecular dynamics trajectory. Some of the properties monitored over the time course of the simulation include the total kinetic and potential energies of the system, the potential energy of the protein and water molecules, the hydrogen bonding patterns and the radius of gyration of the protein. In Figure 2 the variations in the kinetic and potential energies of the protein plus the 152 waters of crystallization and the potential energy of just the protein have been compared. Initially the potential energy was quite low (Figure 2), indicating that the starting structure was sufficiently minimized. As the temperature of the system increased there was a concurrent and rapid rise in the potential energy, attaining a maximum value of 1598 kcal/mol at - 1 2 ps. Subsequently, there was a graded exponential decrease in the potential energy, reaching a convergent value at ~ 100 ps. The time course changes observed for the total system energy directly parallel those described for the potential energy of the system (Figure 2). In marked contrast, the system kinetic energy reached a steady state value of — 2700 kcal/mol following the initial system warming. Examination of the protein contribution to the potential energy showed little change in its potential energy after the initial system warm-up (Figure 2); however, there was a considerable decrease in the
190
200 290 Tine (put)
Fig. 2. Calculated potential energy ( • ) and kinetic energy ( • ) as a function of time for all waters and the protein in the molecular dynamics simulation. Calculated potential energy of the protein only (A). Energy units are in kcal/mol.
water contribution to the potential energy through the first 100 ps of the simulation (Figure 2). This decrease in water potential energy directly correlates with the loss in total potential energy of the system. Changes observed in the potential energy, through the first 100 ps, can be directly attributed to the changes seen in the water-water/water-protein hydrogen bonding network. For two atoms to be counted as a hydrogen bond pair they must meet the following criteria: (i) an acceptor atom hydrogen distance of 00
190
J00
190
400
Fig. 3. Number of total (A), conserved total ( • ) , water-water ( • ) , water-protein (A), backbone (O) and side chain (D) hydrogen bonds found in T4L as a function of time.
705
Downloaded from http://peds.oxfordjournals.org/ at University of California, San Francisco on December 5, 2014
others by the lone a-helix within the domain and short stretches of unordered structure on either side of the helix. The C-terminal domain consists of two disjunct regions, residues 3 — 11 and residues 81-162 (Faber and Matthews, 1990). This domain is predominantly a-helical, except for a /3-turn at position 159-162 and a 3 10 helix at position 108 — 113, helix F. The a-helices are A ( 3 - l l ) , D(82-90), E(93-106), G(115-123), H(126-134), 1(137-141) and J(143—155) (Nicholson et al., 1988). The helices are interconnected by short 2 — 3 residue sections of unclassified structure. Helices E, G, H and J pack together in a way reminiscent of an up-and-down four helix cluster (Richardson and Richardson, 1987). The two domains are joined by the 21 amino acid C-helix consisting of residues 60-80. The overall secondary structure composition of T4L is 68% a-helix, 11% /3-sheet, 10% /3-turns, 7% unclassified structure and 4% 3i0 helix.
G.E.Arnold and R.L.Ornstrin
706
TabJe I. Hydrogen bond comparisons between the time average and X-ray structure X-ray
Time
average structure Total number of hydrogen bonds Conserved hydrogen bonds Intermolecular hydrogen bonds Conserved intermolecular hydrogen bonds Water-water hydrogen bonds Water—protein hydrogen bonds Backbone-backbone hydrogen bonds Conserved backbone—backbone hydrogen bonds Side chain hydrogen bonds Conserved side chain hydrogen bonds
229 91 27 64 99 39
279 93 173 6 132 41 87 78 19 9
i y?
Fig. 4. Deviation (r.m.s.) in position away from the crystallographic coordinates for main chain atoms (—), heavy atoms ( ) and main chain atoms in secondary structure units ( ) as a function of time for T4L
Fig. 5. Time course of the calculated radius of gyration for the molecular dynamics simulation of T4L
warm-up portion of the simulation, the radius of gyration remained stable throughout the remainder of the simulation. The radius of gyration averaged over the last 350 ps simulation was 16.62 ± 0.06 A which is in good agreement with the X-ray crystal value of 16.4 A. Thus there is no decrease in packing
Downloaded from http://peds.oxfordjournals.org/ at University of California, San Francisco on December 5, 2014
of the simulation there were typically an average of 80 ± 4 backbone-backbone hydrogen bonds and an average of 73 ± 4 conserved backbone hydrogen bonds. Thus, on average, - 7 2 % of the original 102 backbone—backbone hydrogen bonds were retained through the last 350 ps of the simulation. However, if the weak hydrogen bonds that were lost during the initial 10 ps of the simulation are neglected, - 9 1 % of the average 80 backbone-backbone hydrogen bonds are conserved during the remainder of the trajectory. There is an - 6 % SD in both the mean backbone-backbone and the mean conserved backbone hydrogen bonds indicating that some backbone hydrogen bond exchange occurs during the simulation. Contrasting the backbone-backbone hydrogen bonds in Figure 3, there is only a slight discernible decrease in the number of side chain hydrogen bonds. There were 39 side chain hydrogen bonds in the initial energy-minimized structure which decreased to an average value of 32 ± 4 over the last 350 ps of the trajectory. The number of conserved side chain hydrogen bonds is markedly lower than the average number of side chain hydrogen bonds, with an average value of 11 ± 2 calculated over the last 350 ps of the simulation. This number is one-third lower than the average number of side chain hydrogen bonds and there is an 18% SD in its average. These values are indicative of a much higher hydrogen bond exchange rate associated with the protein side chains. An average dynamics structure for T4L was calculated by averaging the atomic coordinates of all atoms over the final 350 ps of the trajectory. In Table I the hydrogen bonding patterns between the time average dynamics structure and the crystal structure are compared. There is an increase in the total number of hydrogen bonds from 229 to 279 with only - 4 1 % of these hydrogen bonds being conserved between the two structures. The average dynamics structure shows a marked increase in both the number of intermolecular hydrogen bonds and the number of water-water hydrogen bonds. These results can be attributed to the clustering of water molecules into discrete groups around the protein. There is also a substantial decrease in side chain hydrogen bonds and conserved side chain hydrogen bonds which decrease by - 5 0 % and - 7 5 % respectively. The time average structure shows a slight decrease in the number of backbonebackbone hydrogen bonds, although - 7 9 % of these hydrogen bonds were conserved between the two structures. The majority of the non-conserved backbone-backbone hydrogen bonds, in the time average structure, were located, not surprisingly, in turns and regions of irregular secondary structure. In Figure 4 the time course changes in the r.m.s. deviation in a position away from the crystallographic coordinates of all non-hydrogen atoms, backbone atoms and backbone atoms in secondary structure units are examined. All three curves exhibit a rapid initial increase in r.m.s. deviation away from the starting crystal structure within the first 8.5 ps of the simulation. The heavy atom r.m.s. deviation increases in value rising from 0.96 to 1.75 A, the backbone atoms from 0.78 to 1.56 A and the backbone atoms in secondary structure units from 0.69 to 1.45 A. This is followed by several large changes in local maxima and minima values until a steady state value is reached at - 5 0 ps. The r.m.s. deviation away from the crystal structure averaged over the last 350 ps is 1.7 ± 0.1 A for the heavy atoms, 1.4 ± 0.1 for the backbone atoms and 1.2 ± 0.1 for the backbone atoms in secondary structure units. The radius of gyration was used as a general measure to compare the dynamic structures with that of the X-ray structure (Figure 5). After an initial 2.5% increase during the early
Molecular dynamics of T4 lysozyroe
In Figure 6, the reported X-ray temperature factors (Weaver and Matthews, 1987) are compared with the temperature factors determined from the simulation. Figure 6A displays the theoretical temperature factors calculated from the isotropic r.m.s. fluctuations of the heavy atoms and Figure 6B shows the isotropic
ihTiviiiiiii irwi iiiiiiiri>
Fig. 6. Thermal factor B of the heavy atoms (A) and backbone atoms (B) averaged per residue and plotted as a function of residue number. Dark trace is calculated from the dynamics simulation and the stippled trace is determined from the X-ray structure. Bars along the abscissa demarcate units of secondary structure.
r.m.s. fluctuations calculated from only the backbone atoms. Both sets of data were averaged on a per residue basis and plotted as a function of residue number. Consistent with several other molecular dynamics simulations, the experimental values are, by and large, greater than the theoretical values. The average value for the theoretical temperature factor for all heavy atoms is 16.06 A2 and for the backbone atoms is 10.25 A2, whilst the respective experimental values are 23.89 and 20.29 A2. There is, in general, good qualitative agreement between the experimental and theoretical values, demonstrated by the coincident peaks and valleys in both plots. The theoretical and experimental graphs exhibit a common pattern; however, there are some noteworthy exceptions. In the simulated molecule, Lys48 exhibits a peak in both the heavy atom and backbone atom temperature factor plots which is absent from the experimental data. The simulated backbone temperature factors are quite large from residue 47 to 57 and the majority of these residues are associated with a surface turn and unordered structure, except for residues 4 7 - 5 0 ; these residues are located at the C-terminal end of the B helix. Visual inspection of the trajectory revealed especially prominent fluctuations in Lys48 and a significant amount of fraying at the end of the B helix. The conspicuously large peaks observed for the theoretical values are all located on the protein surface. Also, there are many crystal contact points in this region; in particular residues 48, 49, 51 and 55 all make inter-protein contacts (Weaver and Matthews, 1987). Given these attributes it is not surprising that this particular area exhibits high mobility in the dynamics simulation. Interestingly, the calculated temperature factor of d u l l , in the modeled system, is quite large which is significantly different from the experimental data. This residue is the terminal amino acid in the A helix and is also situated at the junction between the N- and C-terminal domains. It has been suggested that a shift in the experimental position of Glul 1 is necessary for the enzyme to become catalytically competent (Anderson et al., 1981). Whether the high temperature factor calculated reflects some sort of an adjustment that the molecule is trying to make to become catalytically active, or whether this high value is merely an artifact of the residue being located at the edge of a secondary structure element, or whether crystal packing forces have induced a strained conformation is unclear. In general, the backbone atoms of the simulated molecule exhibit less fluctuation in well ordered secondary structure units than do the experimental backbone atoms; in the simulated structure, the largest fluctuations are found at the boundaries of highly ordered secondary structure segments. A case in point is the large C helix which connects the two domains; note the low and constant r.m.s. fluctuations through the helix proper and the spikes at either end of the helix. Contrast this to the experimental thermal values where a gradual decrease is seen from both ends of the helix, subsequendy reaching a minimum value in the middle of the helix proper. Also, note that many of the prominent peaks present in the heavy atom temperature factor plot are absent in the backbone atom temperature factor plot and that these discrepancies can be attributed to large fluctuations contributed by side chain atoms (e.g. Lysl35). 4> and * dihedral angle fluctuations were calculated in an attempt to examine the backbone motions of the protein. The dihedral angle fluctuations are plotted as a function of residue number in Figure 7. The average r.m.s. fluctuation in and * for all residues is 12.2° and the average r.m.s. fluctuation in w for all residues is 7.3°. Residues found in secondary structure elements exhibit the lowest r.m.s.fluctuations(e.g. those residues 707
Downloaded from http://peds.oxfordjournals.org/ at University of California, San Francisco on December 5, 2014
density of the protein using our molecular dynamics parameters, but rather, a slight increase upon addition of thermal energy to the system. Experimentally measured values obtained from diffuse X-ray scattering studies indicate that the radius of gyration of the molecule is ~ 18.6 A in solution. This value is 13% greater than our calculated value and the crystal coordinate value, and suggests there is a difference between the protein in solution and in the crystal (Timchenko et al., 1978). The small increase in the radius of gyration of T4L contrasts the slight decrease in the radius of gyration seen for hen egg-white lysozyme (Post et al., 1986) and carboxypeptidase A (Makinen et al., 1989); however, the slight increase in radius of gyration for T4L is consistent with other protein simulations (Swaminathan et al., 1991; Bass et al., 1992). The T4L simulation used an all-atom model, in contrast to the earlier hen egg white lysozyme and carboxypeptidase A simulations which used a united-atom approximation for non-polar hydrogens. Atomic fluctuations of the simulation
G.E.Arnold and R.L.Ornstein 60
r
-
Table II. Average r m s deviation of the secondary structure elements in T4L
mtdiM mmibtr
found in helix and /3-sheet structures have average r.m.s. fluctuations of 10.7° and 12.0°respectively).In contrast, residues which are neither in helix or /3-sheet structures have an average r.m.s. fluctuation of 18.9°. Examination of Figure 7 shows that the largest rotational r.m.s. fluctuations occur at the boundaries of secondary structure elements and all of these fluctuations occur in 3> and ¥ pairs between adjacent residues. Visual analysis of the trajectory shows that these large r.m.s. fluctuations arise from concerted anticorrelated rotations between adjacent ,¥ pairs, resulting in virtually no perturbation in the direction of the protein backbone. Similar observations have been noted in other simulations (Post et al., 1986; Paulsen et al., 1991). Secondary structure analysis In Table II the r.m.s. deviations away from the starting crystal structure, averaged over the final 350 ps of the simulation, are listed for the individual secondary structure elements. The mean r.m.s. deviation for the helical elements is 0.40 ± 0.08 A. The largest r.m.s. deviation for the helical units is found in helices A and C. There exist three distinct loci where hinge bending is thought to take place and, interestingly, both of these helices are localized at or near the putative hinge bending region of the molecule (Faber and Matthews, 1990). The first of these is in an extended section contiguous with the C-terminal end of helix A which links the C-terminal domain to the N-terminal domain. The other two sites are near the boundary regions at the two ends of the long 20 residue C helix which connects the two domains. In addition, there is a relatively large mean r.m.s. deviation of 1.12 exhibited by turn 4. This turn is located near the N-terminus of the C helix and is part of one of the purported hinge bending loci. In Figure 8 and Table III the time development of the overall atomic fluctuations in the simulated molecule is examined. The atomic fluctuations of the molecule as a whole will be used to benchmark the atomic fluctuations of individual structural elements. Average values of the r.m.s. isotropic displacements were calculated for all heavy atoms and for the backbone atoms N, C, Ca and O. Not surprisingly, the values calculated for the heavy atoms were significantly greater than those calculated for the backbone atoms because of the inclusion of side chain fluctuations. Other than the magnitudes of thefluctuations,the general characteristics of the two curves appear to be nearly identical, even though large fluctuations from side chain atoms are included in the heavy atom calculation. Initially, the increase in fluctuations are rapid and nearly linear; thefluctuationsreach - 82% 708
Structural unit
R m s deviation (A)
4-162 3-11 39-50 60-80 82-90 93-106 108-113 115-123 126-134 137-141 143-155 56-58 14-20 24-27 31-34 20-23 28-30 54-57 159-162
molecule helix A (a) helix B (a) helix C (a) helix D (a) helix E (a) helix F (3,0) heilx G (a) helix H (a) helix 1 (a) helix J (a) strand 1 strand 2 strand 3 strand 4 turn 1 turn 2 turn 3 turn 4
1 35 0 53 0 39 0.52 0 33 0 48 0 42 0 37 0.34 0.29 0.38 0.55 0.44 0.28 0.26 0 39 0 22 1 12 0 47
Fig. 8. Time development of isotropic atomic fluctuations in T4L. The data points indicate the time interval over which the fluctuations were calculated The calculations were averaged for all heavy atoms ( • ) and backbone atoms (O).
of their maximal value within the 10 ps average. Time average fluctuations after 50 ps increase slowly but constantly and the functions do not appear to be approaching an asymptotic limit. As the time period is increased, the time development of the fluctuations should increase smoothly and rapidly to a convergent value. Because an asymptotic limit is not achieved, the equilibriumfluctuationof the overall protein has not yet been reached at 400 ps. This can be attributed to individual structural elements which have not yet converged by 400 ps. In general, for each time block average the atomic fluctuations for the heavy atoms was —0.05 A greater than for the backbone atoms. To assess therelativeflexibilitiesof the helices, sheets and turns in T4L, the time development of atomic fluctuations for the various structural elements was determined. This was done by averaging the r.m.s. positional fluctuations over the backbone atoms for each residue comprising a particular secondary structural element (Post et al.. 1989). Time development of the fluctuations was monitored by varying the size of the time period
Downloaded from http://peds.oxfordjournals.org/ at University of California, San Francisco on December 5, 2014
Fig. 7. The average r.m.s. fluctuations of the backbone dihedral angles * In order to facilitate plotting of both dihedral angles on the same graph, the negative value of • was plotted.
Residues
Molecular dynamics of T4 lysozyme
Table i n . Positional fluctuations of T4L calculated over varying time blocks for the indicated structural elements and specified residues Residues Structural unit
Time period (ps)
.
A molecular dynamics simulation of bacteriophage T4 lysozyme
Gregory E.Arnold and Rick L.Ornstein' Molecular Science Research Center, Pacific Northwest Laboratory, PO Box 999, K2-18, Richland, WA 99352, USA 'To whom correspondence should be addressed
Introduction Experimental and theoretical studies have contributed to an increase in the understanding of the relationships between structure, function and dynamics in proteins and macromolecules in general. Molecular dynamics simulations provide a tool to examine both the time scale and magnitudes of atomic fluctuations. It has been suggested that fast thermal atomic fluctuations propagate slower motions and that these generate yet slower motions (Ansari et al., 1985; Dommair and Jahnig, 1989). Taken together, these architectonic levels of motion are thought
703
Downloaded from http://peds.oxfordjournals.org/ at University of California, San Francisco on December 5, 2014
An analysis of a 400 ps molecular dynamics simulation of the 164 amino acid enzyme T4 lysozyme is presented. The simulation was carried out with all hydrogen atoms modeled explicitly, the inclusion of all 152 crystallographic waters and at a temperature of 300 K. Temporal analysis of the trajectory versus energy, hydrogen bond stability, r.m.s. deviation from the starting crystal structure and radius of gyration, demonstrates that the simulation was both stable and representative of the average experimental structure. Average structural properties were calculated from the enzyme trajectory and compared with the crystal structure. The mean value of the Ca displacements of the average simulated structure from the X-ray structure was 1.1 ± 0.1 A; differences of the backbone $ and * angles between the average simulated structure and the crystal structure were also examined. Thermal-B factors were calculated from the simulation for heavy and backbone atoms and both were in good agreement with experimental values. Relationships between protein secondary structure elements and internal motions were studied by examining the positional fluctuations of individual helix, sheet and turn structures. The structural integrity in the secondary structure units was preserved throughout the simulation; however, the A helix did show some unusually high atomic fluctuations. The largest backbone atom r.m.s. fluctuations were found in nonsecondary structure regions; similar results were observed for r.m.s. fluctuations of non-secondary structure $ and ¥ angles. In general, the calculated values of r.m.s. fluctuations were quite small for the secondary structure elements. In contrast, surface loops and turns exhibited much larger values, being able to sample larger regions of conformational space. The Ca difference distance matrix and superpositioning analyses comparing the X-ray structure with the average dynamics structure suggest that a 'hinge-bending' motion occurs between the N- and C-terminal domains. Key words: computer simulation/protein dynamics/protein motions/T4 lysozyme
to give rise to lower frequency oscillations within the protein (Ansari et al., 1985). Hence, fast co-operative atomic fluctuations lead to side chain fluctuations, secondary structure fluctuations, backbone motions, domain movements and ultimately conformational changes. Thus, the collective dynamics of these regions are relevant to biological activity for proteins involved in catalysis and transport (Mao et al., 1982; Levitt, 1983; Ansari et al., 1985). Protein tertiary and secondary structures are stabilized by weak non-bonded interactions and, as a result, considerable atomic fluctuations should occur. The study and characterization of these fluctuations can lend insight into the nature of internal protein motions. X-ray crystallography and/or NMR studies primarily provide time average data but are unable to provide detailed information regarding the motional properties of a macromolecule. Assuming that the temporal behavior of proteins can be accurately modeled, molecular dynamics simulation methods can be used to investigate the collective oscillations within the molecule. Molecular dynamics simulation can provide information regarding both the equilibrium properties and the temporal characteristics of a particular molecule. We have used the method of molecular dynamics simulation to study the motional properties of bacteriophage T4 lysozyme (T4L). Other investigators have reported simulations for hen egg white lysozyme (Mao et al., 1982) but no trajectory has been reported and analyzed for T4L to date. Both enzymes are endoacetylmuramidases and both cleave the /3(1—4) glycosidic bond between A'-acetylmuramic acid and A'-acetylglucosamine, although T4L is ~ 250 times more active towards the enzymatic cleavage of Escherichia coli cell walls (Matthews et al., 1981). The two enzymes show little homology yet are thought to be evolutionarily divergent (Matthews et al., 1981). The vast amount of experimental data on the structure and function of T4L and therelativelysmall size of this enzyme make it an excellent system to study by simulation. T4L is a 164 residue, monomeric ~18.7kDa protein with no disulfide bridges. This enzyme has a globular bi-lobed structure with its active site cleft juxtaposed between the two lobes. The enzyme aids in the digestion of the bacterial cell wall facilitating release of the virus from the host and is produced late in the infection of E.coli by bacteriophage T4 (Travers, 1970). The wild-type structure of T4L has been refined to 1.7 A resolution, with R equal to 0.193, and the inclusion of isotropic thermal parameters (Weaver and Matthews, 1987). This structure serves as the starting point for our analysis. In the native enzyme structure, the catalytic cleft is occluded and thus obstructs substrate access to the active site. Consequently, conformational flexibility and dynamic fluctuations are necessary for substrate binding and enzymatic catalysis to occur (Faber and Matthews, 1990). A particularly important objective of our molecular dynamics study is to investigate the relationships of positional fluctuation relative to die movements of essentially 'rigid' structural units within the protein. In particular, we want to examine whether me methods of molecular dynamics simulation can model relative domain movements in
G.E.Arnold and R.L.Ornstein
T4L. Intra-domain movements have been reported for native T4L (Weaver and Matthews, 1987) and relative domain movements have been observed in five different crystalline allomorphs of the T4L M6I mutant (Faber and Matthews, 1990). In this paper the results of a 400 ps molecular dynamics simulation of T4L are presented. The analyses in this initial report focus on: (i) the global properties of the dynamics simulation; (ii) the relationship between atomic fluctuations and protein structural units; and (iii) the comparison of the X-ray structure to the time average dynamics structure. Of particular interest was the relative movement observed in the two lobes of the enzyme when comparing the time averaged dynamics structure with the X-ray structure.
The X-ray structure with added hydrogens was energyminimized using the method of steepest descents for 500 steps, with the positions of all heavy atoms fixed to remove any artifacts induced by the addition of explicit hydrogens. The structure was further minimized for 500 steps with only the positions of the heavy atoms of the protein fixed to allow the waters of crystallization to relax. Then, the structure was minimized for 5000 steps using the steepest descents method followed by 5000 704
Results and discussion To facilitate discussion of T4L, a brief description of the enzyme's structure is given below. T4L is a monomeric 164 amino acid, 18.7 kDa protein and is composed of two globular domains separated by an active site cleft (Figure 1) (Faber and Matthews, 1990). The N-terminal domain consists of residues 12-59 and comprises one relatively surface-exposed a-helix, helix B (39-50) and a four-stranded antiparallel /3-sheet. The first three strands of the sheet are contiguous along the sequence and are connected by /3-turns; the fourth strand is separated from the
Fig. 1. Ribbon trace of the X-ray crystal structure of T4L. Amino acid number* important to the discussion are labeled
Downloaded from http://peds.oxfordjournals.org/ at University of California, San Francisco on December 5, 2014
Materials and methods The coordinates of the reported X-ray structure of T4L at 1.7 A resolution (Weaver and Matthews, 1987), obtained from the Brookhaven Protein Data Bank (Bernstein et al., 1977), served as the initial molecular model for the molecular dynamics calculations. The system simulated comprised 164 amino acid residues and 152 crystallographic waters. All of the hydrogen atoms were explicitly modeled, bringing the total number of atoms in the system to 3099. The coordinates of the added hydrogens were generated according to idealized bond lengths and valence angles using the MOLEDT software package from Biosym Technologies (San Diego, CA). The molecular dynamics simulation was performed using the Discover simulation software package (version 2.41) from Biosym Technologies on a Cray X-MP. The default Discover model was used for the explicit (crystallographic) waters. No constraints were applied to the crystallographic waters and thus they were free to move from their crystallographically determined sites. No cross-terms were used in the energy expression and a simple harmonic potential was chosen for the bond stretching terms. All calculations were conducted with a group-based non-bonded cut-off of 10.5 A imposed over a switching distance of 1.5 A and the non-bonded pair list was updated every 20 time steps. A linear distance-dependent dielectric (equal to the interatomic separation) was used. The parameters used were those of the consistent valence force field (VFF) (Hagler, 1985; Dauber-Osguthorpe et al., 1988) except that the charges of acidic, basic, N-termini and C-termini 'functional groups' were made net neutral; although atoms were assigned net charges similar to their normally uncharged polar counterparts. This was done to compensate for the missing dielectric effect of bulk water, because only the crystallographic waters were included. This type of charge screening effect has been successfully used by others (Aqvist et al., 1985; Kruger et al., 1985; Makinen et al., 1989; Mehler and Solmajer, 1991; Paulsen et al., 1991; Solmajer and Mehler, 1991; Swaminathan et al., 1991; Braatz et al., 1992; Gu and Brady, 1992). In previous simulations of crambin, using net neutral residues with various dielectric models, a distance-dependent dielectric or a dielectric of 1 yielded nearly identical results and were superior to other dielectric models (Ornstein, 1990).
steps using the conjugate gradients method before beginning molecular dynamics. Dynamics was performed using the leapfrog algorithm with a 1 fs time step. A constant temperature for the simulation was maintained by weak coupling to a thermal bath (Berendsen et al., 1984). For the dynamics, 0.5 ps was first performed at 50 K. This was followed by elevating the temperature to 300 K, using an exponential approach (time constant of 2.0 ps), such that the target temperature was achieved after 10 ps. At 10 ps the time constant was adjusted to 0.1 ps and the dynamics was re-initialized at 300 K. The trajectory was continued for a total of 400 ps. The calculation was performed on a Cray X-MP computer and 1 c.p.u. h was required per 4.5 ps. The simulation was restarted at 100 and 200 ps with the previous structure used as the current starting structure and with new random velocities. The data set considered in this work comprised 801 structures saved at 0.5 ps intervals.
Molecular dynamics of T4 lysozyme
Time course of the simulation Several parameters have been examined that assess the stability and the convergence characteristics of the molecular dynamics trajectory. Some of the properties monitored over the time course of the simulation include the total kinetic and potential energies of the system, the potential energy of the protein and water molecules, the hydrogen bonding patterns and the radius of gyration of the protein. In Figure 2 the variations in the kinetic and potential energies of the protein plus the 152 waters of crystallization and the potential energy of just the protein have been compared. Initially the potential energy was quite low (Figure 2), indicating that the starting structure was sufficiently minimized. As the temperature of the system increased there was a concurrent and rapid rise in the potential energy, attaining a maximum value of 1598 kcal/mol at - 1 2 ps. Subsequently, there was a graded exponential decrease in the potential energy, reaching a convergent value at ~ 100 ps. The time course changes observed for the total system energy directly parallel those described for the potential energy of the system (Figure 2). In marked contrast, the system kinetic energy reached a steady state value of — 2700 kcal/mol following the initial system warming. Examination of the protein contribution to the potential energy showed little change in its potential energy after the initial system warm-up (Figure 2); however, there was a considerable decrease in the
190
200 290 Tine (put)
Fig. 2. Calculated potential energy ( • ) and kinetic energy ( • ) as a function of time for all waters and the protein in the molecular dynamics simulation. Calculated potential energy of the protein only (A). Energy units are in kcal/mol.
water contribution to the potential energy through the first 100 ps of the simulation (Figure 2). This decrease in water potential energy directly correlates with the loss in total potential energy of the system. Changes observed in the potential energy, through the first 100 ps, can be directly attributed to the changes seen in the water-water/water-protein hydrogen bonding network. For two atoms to be counted as a hydrogen bond pair they must meet the following criteria: (i) an acceptor atom hydrogen distance of 00
190
J00
190
400
Fig. 3. Number of total (A), conserved total ( • ) , water-water ( • ) , water-protein (A), backbone (O) and side chain (D) hydrogen bonds found in T4L as a function of time.
705
Downloaded from http://peds.oxfordjournals.org/ at University of California, San Francisco on December 5, 2014
others by the lone a-helix within the domain and short stretches of unordered structure on either side of the helix. The C-terminal domain consists of two disjunct regions, residues 3 — 11 and residues 81-162 (Faber and Matthews, 1990). This domain is predominantly a-helical, except for a /3-turn at position 159-162 and a 3 10 helix at position 108 — 113, helix F. The a-helices are A ( 3 - l l ) , D(82-90), E(93-106), G(115-123), H(126-134), 1(137-141) and J(143—155) (Nicholson et al., 1988). The helices are interconnected by short 2 — 3 residue sections of unclassified structure. Helices E, G, H and J pack together in a way reminiscent of an up-and-down four helix cluster (Richardson and Richardson, 1987). The two domains are joined by the 21 amino acid C-helix consisting of residues 60-80. The overall secondary structure composition of T4L is 68% a-helix, 11% /3-sheet, 10% /3-turns, 7% unclassified structure and 4% 3i0 helix.
G.E.Arnold and R.L.Ornstrin
706
TabJe I. Hydrogen bond comparisons between the time average and X-ray structure X-ray
Time
average structure Total number of hydrogen bonds Conserved hydrogen bonds Intermolecular hydrogen bonds Conserved intermolecular hydrogen bonds Water-water hydrogen bonds Water—protein hydrogen bonds Backbone-backbone hydrogen bonds Conserved backbone—backbone hydrogen bonds Side chain hydrogen bonds Conserved side chain hydrogen bonds
229 91 27 64 99 39
279 93 173 6 132 41 87 78 19 9
i y?
Fig. 4. Deviation (r.m.s.) in position away from the crystallographic coordinates for main chain atoms (—), heavy atoms ( ) and main chain atoms in secondary structure units ( ) as a function of time for T4L
Fig. 5. Time course of the calculated radius of gyration for the molecular dynamics simulation of T4L
warm-up portion of the simulation, the radius of gyration remained stable throughout the remainder of the simulation. The radius of gyration averaged over the last 350 ps simulation was 16.62 ± 0.06 A which is in good agreement with the X-ray crystal value of 16.4 A. Thus there is no decrease in packing
Downloaded from http://peds.oxfordjournals.org/ at University of California, San Francisco on December 5, 2014
of the simulation there were typically an average of 80 ± 4 backbone-backbone hydrogen bonds and an average of 73 ± 4 conserved backbone hydrogen bonds. Thus, on average, - 7 2 % of the original 102 backbone—backbone hydrogen bonds were retained through the last 350 ps of the simulation. However, if the weak hydrogen bonds that were lost during the initial 10 ps of the simulation are neglected, - 9 1 % of the average 80 backbone-backbone hydrogen bonds are conserved during the remainder of the trajectory. There is an - 6 % SD in both the mean backbone-backbone and the mean conserved backbone hydrogen bonds indicating that some backbone hydrogen bond exchange occurs during the simulation. Contrasting the backbone-backbone hydrogen bonds in Figure 3, there is only a slight discernible decrease in the number of side chain hydrogen bonds. There were 39 side chain hydrogen bonds in the initial energy-minimized structure which decreased to an average value of 32 ± 4 over the last 350 ps of the trajectory. The number of conserved side chain hydrogen bonds is markedly lower than the average number of side chain hydrogen bonds, with an average value of 11 ± 2 calculated over the last 350 ps of the simulation. This number is one-third lower than the average number of side chain hydrogen bonds and there is an 18% SD in its average. These values are indicative of a much higher hydrogen bond exchange rate associated with the protein side chains. An average dynamics structure for T4L was calculated by averaging the atomic coordinates of all atoms over the final 350 ps of the trajectory. In Table I the hydrogen bonding patterns between the time average dynamics structure and the crystal structure are compared. There is an increase in the total number of hydrogen bonds from 229 to 279 with only - 4 1 % of these hydrogen bonds being conserved between the two structures. The average dynamics structure shows a marked increase in both the number of intermolecular hydrogen bonds and the number of water-water hydrogen bonds. These results can be attributed to the clustering of water molecules into discrete groups around the protein. There is also a substantial decrease in side chain hydrogen bonds and conserved side chain hydrogen bonds which decrease by - 5 0 % and - 7 5 % respectively. The time average structure shows a slight decrease in the number of backbonebackbone hydrogen bonds, although - 7 9 % of these hydrogen bonds were conserved between the two structures. The majority of the non-conserved backbone-backbone hydrogen bonds, in the time average structure, were located, not surprisingly, in turns and regions of irregular secondary structure. In Figure 4 the time course changes in the r.m.s. deviation in a position away from the crystallographic coordinates of all non-hydrogen atoms, backbone atoms and backbone atoms in secondary structure units are examined. All three curves exhibit a rapid initial increase in r.m.s. deviation away from the starting crystal structure within the first 8.5 ps of the simulation. The heavy atom r.m.s. deviation increases in value rising from 0.96 to 1.75 A, the backbone atoms from 0.78 to 1.56 A and the backbone atoms in secondary structure units from 0.69 to 1.45 A. This is followed by several large changes in local maxima and minima values until a steady state value is reached at - 5 0 ps. The r.m.s. deviation away from the crystal structure averaged over the last 350 ps is 1.7 ± 0.1 A for the heavy atoms, 1.4 ± 0.1 for the backbone atoms and 1.2 ± 0.1 for the backbone atoms in secondary structure units. The radius of gyration was used as a general measure to compare the dynamic structures with that of the X-ray structure (Figure 5). After an initial 2.5% increase during the early
Molecular dynamics of T4 lysozyroe
In Figure 6, the reported X-ray temperature factors (Weaver and Matthews, 1987) are compared with the temperature factors determined from the simulation. Figure 6A displays the theoretical temperature factors calculated from the isotropic r.m.s. fluctuations of the heavy atoms and Figure 6B shows the isotropic
ihTiviiiiiii irwi iiiiiiiri>
Fig. 6. Thermal factor B of the heavy atoms (A) and backbone atoms (B) averaged per residue and plotted as a function of residue number. Dark trace is calculated from the dynamics simulation and the stippled trace is determined from the X-ray structure. Bars along the abscissa demarcate units of secondary structure.
r.m.s. fluctuations calculated from only the backbone atoms. Both sets of data were averaged on a per residue basis and plotted as a function of residue number. Consistent with several other molecular dynamics simulations, the experimental values are, by and large, greater than the theoretical values. The average value for the theoretical temperature factor for all heavy atoms is 16.06 A2 and for the backbone atoms is 10.25 A2, whilst the respective experimental values are 23.89 and 20.29 A2. There is, in general, good qualitative agreement between the experimental and theoretical values, demonstrated by the coincident peaks and valleys in both plots. The theoretical and experimental graphs exhibit a common pattern; however, there are some noteworthy exceptions. In the simulated molecule, Lys48 exhibits a peak in both the heavy atom and backbone atom temperature factor plots which is absent from the experimental data. The simulated backbone temperature factors are quite large from residue 47 to 57 and the majority of these residues are associated with a surface turn and unordered structure, except for residues 4 7 - 5 0 ; these residues are located at the C-terminal end of the B helix. Visual inspection of the trajectory revealed especially prominent fluctuations in Lys48 and a significant amount of fraying at the end of the B helix. The conspicuously large peaks observed for the theoretical values are all located on the protein surface. Also, there are many crystal contact points in this region; in particular residues 48, 49, 51 and 55 all make inter-protein contacts (Weaver and Matthews, 1987). Given these attributes it is not surprising that this particular area exhibits high mobility in the dynamics simulation. Interestingly, the calculated temperature factor of d u l l , in the modeled system, is quite large which is significantly different from the experimental data. This residue is the terminal amino acid in the A helix and is also situated at the junction between the N- and C-terminal domains. It has been suggested that a shift in the experimental position of Glul 1 is necessary for the enzyme to become catalytically competent (Anderson et al., 1981). Whether the high temperature factor calculated reflects some sort of an adjustment that the molecule is trying to make to become catalytically active, or whether this high value is merely an artifact of the residue being located at the edge of a secondary structure element, or whether crystal packing forces have induced a strained conformation is unclear. In general, the backbone atoms of the simulated molecule exhibit less fluctuation in well ordered secondary structure units than do the experimental backbone atoms; in the simulated structure, the largest fluctuations are found at the boundaries of highly ordered secondary structure segments. A case in point is the large C helix which connects the two domains; note the low and constant r.m.s. fluctuations through the helix proper and the spikes at either end of the helix. Contrast this to the experimental thermal values where a gradual decrease is seen from both ends of the helix, subsequendy reaching a minimum value in the middle of the helix proper. Also, note that many of the prominent peaks present in the heavy atom temperature factor plot are absent in the backbone atom temperature factor plot and that these discrepancies can be attributed to large fluctuations contributed by side chain atoms (e.g. Lysl35). 4> and * dihedral angle fluctuations were calculated in an attempt to examine the backbone motions of the protein. The dihedral angle fluctuations are plotted as a function of residue number in Figure 7. The average r.m.s. fluctuation in and * for all residues is 12.2° and the average r.m.s. fluctuation in w for all residues is 7.3°. Residues found in secondary structure elements exhibit the lowest r.m.s.fluctuations(e.g. those residues 707
Downloaded from http://peds.oxfordjournals.org/ at University of California, San Francisco on December 5, 2014
density of the protein using our molecular dynamics parameters, but rather, a slight increase upon addition of thermal energy to the system. Experimentally measured values obtained from diffuse X-ray scattering studies indicate that the radius of gyration of the molecule is ~ 18.6 A in solution. This value is 13% greater than our calculated value and the crystal coordinate value, and suggests there is a difference between the protein in solution and in the crystal (Timchenko et al., 1978). The small increase in the radius of gyration of T4L contrasts the slight decrease in the radius of gyration seen for hen egg-white lysozyme (Post et al., 1986) and carboxypeptidase A (Makinen et al., 1989); however, the slight increase in radius of gyration for T4L is consistent with other protein simulations (Swaminathan et al., 1991; Bass et al., 1992). The T4L simulation used an all-atom model, in contrast to the earlier hen egg white lysozyme and carboxypeptidase A simulations which used a united-atom approximation for non-polar hydrogens. Atomic fluctuations of the simulation
G.E.Arnold and R.L.Ornstein 60
r
-
Table II. Average r m s deviation of the secondary structure elements in T4L
mtdiM mmibtr
found in helix and /3-sheet structures have average r.m.s. fluctuations of 10.7° and 12.0°respectively).In contrast, residues which are neither in helix or /3-sheet structures have an average r.m.s. fluctuation of 18.9°. Examination of Figure 7 shows that the largest rotational r.m.s. fluctuations occur at the boundaries of secondary structure elements and all of these fluctuations occur in 3> and ¥ pairs between adjacent residues. Visual analysis of the trajectory shows that these large r.m.s. fluctuations arise from concerted anticorrelated rotations between adjacent ,¥ pairs, resulting in virtually no perturbation in the direction of the protein backbone. Similar observations have been noted in other simulations (Post et al., 1986; Paulsen et al., 1991). Secondary structure analysis In Table II the r.m.s. deviations away from the starting crystal structure, averaged over the final 350 ps of the simulation, are listed for the individual secondary structure elements. The mean r.m.s. deviation for the helical elements is 0.40 ± 0.08 A. The largest r.m.s. deviation for the helical units is found in helices A and C. There exist three distinct loci where hinge bending is thought to take place and, interestingly, both of these helices are localized at or near the putative hinge bending region of the molecule (Faber and Matthews, 1990). The first of these is in an extended section contiguous with the C-terminal end of helix A which links the C-terminal domain to the N-terminal domain. The other two sites are near the boundary regions at the two ends of the long 20 residue C helix which connects the two domains. In addition, there is a relatively large mean r.m.s. deviation of 1.12 exhibited by turn 4. This turn is located near the N-terminus of the C helix and is part of one of the purported hinge bending loci. In Figure 8 and Table III the time development of the overall atomic fluctuations in the simulated molecule is examined. The atomic fluctuations of the molecule as a whole will be used to benchmark the atomic fluctuations of individual structural elements. Average values of the r.m.s. isotropic displacements were calculated for all heavy atoms and for the backbone atoms N, C, Ca and O. Not surprisingly, the values calculated for the heavy atoms were significantly greater than those calculated for the backbone atoms because of the inclusion of side chain fluctuations. Other than the magnitudes of thefluctuations,the general characteristics of the two curves appear to be nearly identical, even though large fluctuations from side chain atoms are included in the heavy atom calculation. Initially, the increase in fluctuations are rapid and nearly linear; thefluctuationsreach - 82% 708
Structural unit
R m s deviation (A)
4-162 3-11 39-50 60-80 82-90 93-106 108-113 115-123 126-134 137-141 143-155 56-58 14-20 24-27 31-34 20-23 28-30 54-57 159-162
molecule helix A (a) helix B (a) helix C (a) helix D (a) helix E (a) helix F (3,0) heilx G (a) helix H (a) helix 1 (a) helix J (a) strand 1 strand 2 strand 3 strand 4 turn 1 turn 2 turn 3 turn 4
1 35 0 53 0 39 0.52 0 33 0 48 0 42 0 37 0.34 0.29 0.38 0.55 0.44 0.28 0.26 0 39 0 22 1 12 0 47
Fig. 8. Time development of isotropic atomic fluctuations in T4L. The data points indicate the time interval over which the fluctuations were calculated The calculations were averaged for all heavy atoms ( • ) and backbone atoms (O).
of their maximal value within the 10 ps average. Time average fluctuations after 50 ps increase slowly but constantly and the functions do not appear to be approaching an asymptotic limit. As the time period is increased, the time development of the fluctuations should increase smoothly and rapidly to a convergent value. Because an asymptotic limit is not achieved, the equilibriumfluctuationof the overall protein has not yet been reached at 400 ps. This can be attributed to individual structural elements which have not yet converged by 400 ps. In general, for each time block average the atomic fluctuations for the heavy atoms was —0.05 A greater than for the backbone atoms. To assess therelativeflexibilitiesof the helices, sheets and turns in T4L, the time development of atomic fluctuations for the various structural elements was determined. This was done by averaging the r.m.s. positional fluctuations over the backbone atoms for each residue comprising a particular secondary structural element (Post et al.. 1989). Time development of the fluctuations was monitored by varying the size of the time period
Downloaded from http://peds.oxfordjournals.org/ at University of California, San Francisco on December 5, 2014
Fig. 7. The average r.m.s. fluctuations of the backbone dihedral angles * In order to facilitate plotting of both dihedral angles on the same graph, the negative value of • was plotted.
Residues
Molecular dynamics of T4 lysozyme
Table i n . Positional fluctuations of T4L calculated over varying time blocks for the indicated structural elements and specified residues Residues Structural unit
Time period (ps)
.
