Section 20.1.3. Results

In Fig. 20.1.3.1, the non-bonded contributions to the total potential energy are shown. The non-bonded interactions comprise Lennard–Jones and electrostatic interactions. Solvent–solvent, solute–solute and solute–solvent interaction energies are shown separately. All of these appear converged after approximately 100 ps. The solvent–solvent energy remains close to its initial value during the whole simulation, the water molecules having relaxed during the pre-equilibration period, while the protein was restrained. The protein internal energy increases during the first few hundred picoseconds, but this is compensated by a decrease in the protein–solvent energy as the protein adapts to the force field and the pre-relaxed solvent environment. This effect becomes negligible after about 200 ps, from which time point the system can be viewed as equilibrated with respect to the energies. The distribution of kinetic versus potential energy and the total (bonded and non-bonded) energy of the system relaxes even faster (results not shown).

Figure 20.1.3.1| top | pdf | Non-bonded energies (in kJ mol−1) of the simulated system as a function of time.

Not all properties converge as fast as the energies. Fig. 20.1.3.2 shows the root-mean-square atom-position deviation (RMSD) from the X-ray structure for each of the four individual chains for both Cα atoms and all atoms. Here, several relaxation periods can be distinguished. After the initial increase, which occurs during the first 50 ps of the simulation, a plateau is reached, and the system is apparently stable until 300 ps. The values reached are 0.12 nm for the Cα atoms and 0.20 nm if all atoms are considered. These numbers are comparable with results obtained in crystal simulations of other proteins of equivalent length reported in the literature (van Gunsteren et al., 1983; Berendsen et al., 1986; Shi et al., 1988; Heiner et al., 1992; Levitt et al., 1995). After 300 ps, however, the values increase slowly again. For the Cα atoms, there is apparently a second plateau from 300 to 850 ps, but during this period the RMSD for all atoms continues to increase monotonically. After 850 ps, a final plateau is reached. During the second nanosecond of the simulation (1000–2000 ps), the RMSDs are 0.21 nm for the Cα atoms and 0.29 nm for all atoms. The RMSD of chain 1 is an exception. There is a strong increase after 1200 ps due to a movement of a particular part of the chain which will be addressed later. To ensure that the RMSD values have converged, longer runs would be required.

Figure 20.1.3.2| top | pdf | Root-mean-square atom-positional deviations (RMSD) in nm from the X-ray structure of the four different protein molecules in the unit cell as a function of time. Rotational and translational fitting was applied using the Cα atoms of residues 1–72. The upper and lower graphs show the deviations for the Cα atoms and for all atoms, respectively.

Although the RMSD values shown in Fig. 20.1.3.2 grow larger than those usually observed in the course of short simulations, the hydrogen-bonding pattern and thus the secondary-structure elements observed in the X-ray structure are reproduced well (Table 20.1.3.1). Most of the hydrogen bonds reported (Vijay-Kumar et al., 1987) show high occupancies during the whole simulation, especially the ones within secondary-structure elements. Only six out of the 44 hydrogen bonds present in the X-ray structure are disrupted during the simulation. Hydrogen bonds in the α-helix (residues 23–34) show very high occupancies, ranging from 75% at the N-terminus to well over 90% inside the helix. Only its C-terminus shows signs of instability, with the α-helix deforming towards a 3₁₀-helix. The β-sheet pattern is, apart from chain 1 in the region of residues 49–64, as stable as the α-helix. Occupancies range from 55% up to 95%. The six hydrogen bonds not reproduced (five 3₁₀ and one α-helical) can be rationalized as follows. The bridges 10–7 and 65–62 are part of the most mobile regions of the protein. These regions involve residues 7–10, 51–54 and 62–65 (Vijay-Kumar et al., 1987). The hydrogen bond at the C-terminal end of the α-helix (35–31) is lost, and the preceding hydrogen bond (34–30) is partly changed, indicating that the end of the α-helix is deformed towards a 3₁₀-helix. The donor of the bond 40–37 is replaced by residue 41, and the four-residue 3₁₀-helix that was stabilized by hydrogen bonds 58–55 and 59–56 is replaced by an α-helical hydrogen bond 59–55. The high occupancy of this particular bond and the complete absence of the two observed experimentally indicate an early rearrangement in this part of the structure (before the analysis period) which is stable during the rest of the simulation.

Backbone–side-chain hydrogen bonds are less well reproduced than backbone–backbone interactions. While some are present 80–90% of the time, others are present less than 50% of the time. Two out of the seven hydrogen bonds in which a backbone atom is the donor are not observed in the simulation; both involve the OG1 atom of Thr7 as an acceptor. The hydrogen-donor atoms are the backbone nitrogen atoms of residues Thr9 and Lys11, both of which have high experimental B factors (18.32 Ų for Thr9 and 13.56 Ų for Lys11). The mean of the experimental B factors is 10.73 Ų for the backbone atoms and 13.41 Ų for all protein atoms. Where a side-chain atom is the donor, three out of the five hydrogen bonds present in the X-ray structure are not found in the simulation. All of these involve the side-chain nitrogen atom of a lysine residue as the donor, the experimental B factors of which range from 23.92 Ų for the NZ atom of Lys48 up to 30.06 Ų for the NZ atom of Lys33. Of the four side-chain–side-chain hydrogen bonds, not one is observed as in the crystal. The 54–58 hydrogen bond is, however, replaced by a 55–58 hydrogen bond. All the others involve very mobile atoms with large experimental B factors as donors and acceptors.

There is one intermolecular hydrogen bond (Table 20.1.3.2) in the starting structure which is not seen in the simulation. The donor is the side-chain nitrogen atom of Lys6, which has an experimental B factor of 20.55 Ų, and the acceptors are the side-chain oxygen atoms of Glu51, with B factors of 32.13 and 33.44 Ų. Most of the hydrogen bonds not reproduced in the simulation contain at least one mobile atom. Although these atoms do not remain fixed at their equilibrium positions, they may still stabilize the structure on average.

In Fig. 20.1.3.3, the deviation of the Cα atoms of the different chains from the X-ray structure and from the mean MD structure are presented together with the deviation of the mean MD structure from the X-ray structure. Overall, the individual protein molecules remain close to the experimental structure; however, parts of the structure do deviate substantially. The region involving residues 7–11, which experimentally has high B factors (implying high mobility), has, in three out of the four cases, an RMSD for Cα atoms of 0.3 nm or greater. Chain 3, in contrast, is close to the X-ray structure, and the mean MD structure is closer to the X-ray structure than are any of the individual chains. This suggests that the simulation does not deviate systematically from the X-ray structure, but rather that different regions of conformational space are sampled by the different molecules. The same holds for the very mobile region between residues 47 and 64, where chain 1 deviates dramatically from both the mean MD structure and the X-ray structure. In the other chains, this part of the protein remains close to the X-ray structure. Overall, the deviation from the X-ray structure is largest in the C-terminal region. This part is also ill-defined in the experiment, with occupancies of 0.45 for residues 73 and 74, and 0.25 for the terminal two glycines. Other parts of the protein, especially the stable secondary-structure elements, stay close to the X-ray structure. In the average structure, the α-helix, including its C-terminal part which was deformed to a 3₁₀-helix, deviates by a maximum of 0.08 nm from the X-ray structure, although, as seen before, the individual chains may deviate more. The β-sheet regions also stay close to the X-ray structure. As with the helix, residues 1–7, 40–45 and 64–72 stay within 0.1 nm RMSD from the X-ray structure. The β-strands formed by residues 10–17 and 48–50 are not as similar to the experiment, since they lie close to mobile regions and are thus influenced by neighbouring mobile residues. For the strand formed by residues 10–17, from residue 12 onwards the same structural similarity is reached as for all other secondary-structure elements, and residues 48–50 are, again, strongly influenced by the moving part of chain 1.

Figure 20.1.3.3| top | pdf | Root-mean-square Cα-atom-position deviation (RMSD) in nm from a reference structure as a function of the residue number using the final 1.6 ns of the simulation. RMSDs of the four protein molecules in the unit cell from the mean molecular-dynamics structure (dashed line) and from the X-ray structure (solid line) are shown in the first four graphs. The bottom graph shows the RMSD of the mean (over the four molecules) MD structure from the X-ray structure.

In Fig. 20.1.3.4, the impact of different fitting protocols on atomic mean-square position fluctuations (RMSFs) is examined. B factors are related to mean-square position fluctuations according to where the angle brackets indicate a time or a combined time and ensemble average. Molecule 4 was selected because it is the most stable. RMSFs of the Cα atoms were calculated directly from the simulation trajectory (Fig. 20.1.3.4a) after applying a translational fit using the Cα atoms of residues 1–72, which are well defined in the X-ray structure (Fig. 20.1.3.4b), and after applying both a rotational and a translational fit on residues 1–72 (Fig. 20.1.3.4c). In Fig. 20.1.3.4(d), a translational and rotational fit was applied to all Cα atoms (residues 1–76). Removal of the overall translational component of motion reduces the positional fluctuations by 0.04 nm on average. Only the RMSFs in the proximity of the end of the large α-helix formed by residues 23–34 are not affected by the introduction of translational fitting. In contrast, it is exactly this region where the fluctuations are substantially lowered by introducing an additional rotational fit. The regions before residue 27 and after residue 42 are only slightly affected by the removal of overall rotation. These findings suggest that the entire protein translates by about 0.04 nm, while the α-helix region remains close to its initial position, thus rotating relative to the rest of the protein. Inclusion of the four C-terminal residues in the fitting procedure only affects the RMSFs of these residues and residues in the rotating part of the molecule, indicating that these four residues move together with the rest of the molecule. The atom-position fluctuations obtained by applying a full (rotational and translational) fit are determined by internal motion only. The largest RMSFs for residues 1–72 are 0.12 nm. RMSFs of the two C-terminal glycines are 0.26 nm if the last four residues are excluded from the fitting and 0.22 nm otherwise.

Figure 20.1.3.4| top | pdf | Root-mean-square Cα-atom-position fluctuations (RMSFs) in nm are shown for molecule 4 as a function of residue number. RMSFs are averaged over the final 1.6 ns. In (a), no fitting was applied; in (b), translational fitting was applied using the Cα atoms of residues 1–72; and in (c), the rotational component of motion was also removed. In (d), translational and rotational fitting was applied over all Cα atoms (1–76).

If the same properties are examined but averaged over all the chains, similar trends can be observed (Fig. 20.1.3.5). If no fitting is applied, the RMSFs of 0.24 nm, on average, indicate that the different molecules show relative translation and rotation. After translational fitting is applied, the mean RMSFs drop to 0.18 nm. Thus, the molecules translate within the unit cell. If, in addition, the rotational component of overall motion is removed, the whole helix region is much less mobile than before, and the mean RMSFs drop to 0.14 nm. The same holds for region 47–64, dominated by the rotation of part of chain 1. Fluctuations are generally much larger than before when only chain 4 was observed, again indicating that the distinct chains behave in an uncorrelated way. The size of the peaks of the RMSFs averaged over all chains is around 0.22 nm, compared with 0.27 nm when overall rotation is still present. Thus, in addition to internal rotations, relative rotations of the different molecules occur. If the fit is not applied only to the well defined Cα atoms of residues 1–72, the RMSF value becomes slightly higher – apart from the C-terminal region – but this effect is small, with the mean RMSF staying at 0.14 nm. However, the relative heights of the peaks differ, which shows that it is crucial to define a standard fitting protocol that must not be changed during the course of the analysis.

Figure 20.1.3.5| top | pdf | Root-mean-square Cα-atom-position fluctuations (RMSFs) in nm are shown using the same fitting protocols as in Fig. 20.1.3.4, but averaged over all four protein molecules in the unit cell.

In Fig. 20.1.3.6, we concentrate again on molecule 4. Comparing different averaging periods of 400 and 800 ps with different starting points, it can be seen that, in general, the later the simulation, the less motion observed. During the period 800–1200 ps, only the small region between residues 9 and 12 shows more mobility than that between 400 and 800 ps. In the stable region, the rest of the molecule shows the same mobility as in the earlier time period. In the parts that are most mobile between 400 and 800 ps, the motions decrease significantly after the latter time point, indicating that equilibrium is reached. Focusing on the longer averaging periods, 400–1200 ps versus 1200–2000 ps, we see that over the whole chain mobility is comparable, indicating clear equilibrium as far as internal motions are concerned. The fluctuations during the 400–1200 ps period are of the same size as those of the shorter subperiods, 400–800 ps and 800–1200 ps. They are thus determined by movements on a timescale shorter than 400 ps.

Figure 20.1.3.6| top | pdf | Root-mean-square Cα-atom-position fluctuations (RMSFs) in nm are shown for molecule 4, with full translational and rotational fitting over the Cα atoms of residues 1–72. Different averaging periods are compared. (a) shows the RMSFs for the period 400–800 ps, (b) shows those for the period 800–1200 ps. In (c), the results for the previous two periods are averaged (400–1200 ps), and in (d), the results for the period 1200–2000 ps are shown.

In Fig. 20.1.3.7, the effect of different averaging periods is again considered, but taking the whole unit cell into account. Comparing RMSFs between 400 and 800 ps with those of the following 400 ps period, no significant difference is seen. In fact, contrary to what was observed when chain 4 was examined, somewhat more mobility is evident in the later period compared with the earlier one. This difference shows that although the configurations of the single chains converge rapidly, the different chains visit different parts of phase space. The fact that the fluctuations in the period 400–1200 ps are between those of the two shorter analysis periods is another indication that the single-chain movements are converged on these short timescales. In the last 800 ps of the simulation, the RMSFs are substantially higher than in the 800 ps window before that. All of the peaks can be traced back to one of the single chains. If only one of the four molecules differs strongly from the other three, this one determines the magnitude of the fluctuations of the average. The peak at residue 10 comes from chain 2, the ones around residue 20 and the whole region 47–64 are determined by chain 1. The peak at residue 33 originates in chain 3, which at this point differs substantially from the mean MD structure (Fig. 20.1.3.3).

Figure 20.1.3.7| top | pdf | Root-mean-square Cα-atom-position fluctuations (RMSFs) in nm are shown using the same averaging periods as in Fig. 20.1.3.6, but averaged over all four protein molecules in the unit cell.

Fig. 20.1.3.8 displays the atomic root-mean-square position fluctuations for the Cα atoms of the four protein molecules during the whole analysis period, together with corresponding values obtained using equation (20.1.3.1) and the crystallographic B factors. Rotational and translational fitting was applied using the Cα atoms of residues 1–72, and the fluctuations were averaged over the final 1.6 ns. The mobility of the stable secondary-structure elements in the simulation is comparable with that inferred from the experiment. There is a correlation between the more mobile parts of the proteins in the simulation and large B factors in the X-ray structure, but the magnitude of the fluctuations is overestimated in the simulation. The movements of the single chains can be rationalized as follows. In chain 1, the whole region from Gly47 onwards rotates around a stable axis formed by residues 41–46. This part lies, as do all the flexible regions, on the exterior of the protein. Residues 19 and 20, which are stable in all but this single chain, are in contact with this moving part. This rotation, which tends to compact the protein, occurs during the 200 ps period between 1350 and 1550 ps after the start of the simulation, in which the atom-position RMSD from the X-ray structure increases significantly (Fig. 20.1.3.2). Overall, chain 2 is more stable than chain 1. Nevertheless, the end of the unwinding helix shows large fluctuations. In the course of this deformation, the side-chain nitrogen atom of Lys11 moves from close to the OE atom of Glu34 towards the backbone oxygen atom of Lys33, which is associated with a change in the position of Gly10. A similar but smaller motion occurs in chain 4. Both lysines, Lys33 and Lys63, are fully exposed to the solvent and have no intramolecular contacts. In chain 3, the flexible residues are also not part of secondary-structure elements and are located on the outside of the protein. The backbone oxygen atom of Gln62 that moves in all the four chains has, in addition, the closest contact to another heavy atom: the OG1 atom of Ser65 is only 2.51 Å away, and the van der Waals repulsion of these atoms causes them to move further away from one another. The mobile residues in chain 4 are again in contact with the solvent, Gly35, Gly47, Gln62, the end of the helix and Gly10. The terminal residues of all the protein molecules are very mobile, as observed experimentally in the crystal.

Figure 20.1.3.8| top | pdf | Root-mean-square Cα-atom-position fluctuations (RMSFs) in nm for the four protein molecules in the unit cell as a function of the residue number. Full translational and rotational fitting was applied to the Cα atoms of residues 1–72 using the final 1.6 ns of the simulation [(a)–(d)]. (e) shows the corresponding values defined by equation (20.1.3.1), obtained from experimental B factors.

Backbone dihedral-angle fluctuations and transitions are examined in Tables 20.1.3.3 and 20.1.3.4 using different analysis periods. After the first 400 ps of analysis, the dihedral-angle fluctuations differ only slightly, but if longer averaging times are chosen, the different protein molecules show larger differences from one another. These fluctuations also increase for longer analysis times, indicating that they are not yet converged after 2 ns. In the period from 800 to 1200 ps, chain 3 shows a large increase in mean-square dihedral-angle fluctuations, whereas the Cα-atom-position RMSDs with respect to the X-ray structure during the same time fluctuate around a plateau. Thus, there is a lot of flexibility without the simulation structure diverging from the experimental one. Protein molecule 3, for example, shows the largest fluctuations of all the four molecules, and it shows the lowest atom-position RMSDs of Cα atoms from the X-ray structure at the end of the simulation (Fig. 20.1.3.2), indicating that it explores phase space around the equilibrium structure. If, in contrast, the Cα-atom-position RMSDs, with respect to the X-ray structure, increase significantly, larger dihedral-angle fluctuations are also observed, for example, in molecule 1 after 1200 ps.

Table 20.1.3.3| top | pdf | Root-mean-square fluctuations of polypeptide backbone and ψ dihedral angles (°) for the different molecules using different time-averaging periods The bottom row shows the averages over all four protein molecules. Molecule 400–800 ps 400–1200 ps 400–2000 ps 1 18.4/22.9 19.5/23.7 31.0/33.6 2 17.2/17.0 18.6/18.7 23.6/26.8 3 18.5/20.4 25.6/26.3 35.3/37.5 4 19.7/18.8 19.4/20.3 21.6/28.8 All 26.1/26.2 28.0/28.6 35.2/38.1

Concerning relaxation, observations similar to those made before can be made when analysing dihedral-angle transitions (Table 20.1.3.4). The number of transitions for the different chains differs by about 20%. Within a single chain, however, the number of transitions increases in proportion to the observation time. Again, the protein molecules showing the most transitions do not have the largest Cα-atom-position RMSDs from the X-ray structure. Thus, only certain dihedral-angle flips induce large movements that determine the RMSD value.

In Fig. 20.1.3.9, the number of water oxygen atoms with a given atomic root-mean-square position fluctuation (RMSF) are plotted. The time evolution and the shapes of these curves are similar to those obtained for bulk water, a Gaussian distribution with the maximum at larger RMSF values and larger standard deviations when using longer averaging times. Using a diffusion constant of bulk SPC water at 300 K of 3.9 × 10⁻³ nm² ps⁻¹ (Smith & van Gunsteren, 1995), the root-mean-square position fluctuation for an average water molecule would be 1.25 nm for a 400 ps period, 1.77 nm for a 800 ps period, and 2.5 nm for a 1600 ps period. Comparison of these values with the distributions in Fig. 20.1.3.9 indicates that the motion of most of the crystal water molecules is restricted by the crystalline environment. A number of water molecules adopt well defined positions for periods of the order of some 100 ps. Several water molecules were also observed to visit the same location repeatedly. This indicates that, although not remaining fixed in space, they stay close to a crystallographically determined site which they revisit periodically, alternating with other water molecules.

Figure 20.1.3.9| top | pdf | The number of water molecules with a given root-mean-square oxygen-position fluctuation (RMSF) in nm are shown for different averaging periods: 400–800 ps (solid line), 400–1200 ps (short-dashed line), 400–2000 ps (long-dashed line).