Dihedral Distributions

Analysis of energy fluctuations indicates that the NEIMO method accurately solves its equations of motion for molecular systems, but does not verify that NEIMO produces molecular motions similar to Cartesian dynamics, which solve the equations of motion for the full 3 - 6 degrees of freedom. In order to determine how well NEIMO simulations represent the dynamics of the Met-enkephalin peptides, the distribution of dihedral angles during these simulations was determined. Cartesian dynamics using a 1fs timestep and NEIMO dynamics using a variety of timesteps were run for 5.0 ps at a simulation temperature of 300 K, during which the dihedral angles were output every 0.1 ps. Figure shows the resulting distributions from simulations of Met-enkephalin. The numbering of the dihedral angles is shown in Figure and further identified in Table . The top graph in Figure shows the distribution from Cartesian dynamics and the bottom shows the distribution from a NEIMO dynamics simulation; both simulations used a 1 fs timestep. The distributions from the two simulations are very similar, with the backbone dihedrals (6, 9, 12, and 17) showing the least flexibility, as would be expected, and the methionine sidechain dihedrals showing the greatest variation during the simulation. The average values for each dihedral, , can be calculated from such distributions. Because dihedral angles have a periodicity of 2 (360), the average cannot be calculated directly, but is derived from the average cosine and sine[49]:

Once is known, the standard deviations can be calculated easily from simulations:

where

Because of the periodicity of dihedral angles, Equation () can always be enforced by appropriate additions or subtractions of .

The average values, , and standard deviations, , for the distributions in Figure are shown in Figure . The average values are also shown in Table , and are compared to the initial conformation. The NEIMO results are very similar to the results from the Cartesian simulations, indicating that the reduction in the number of degrees of freedom does not, in general, affect the torsional flexibility of the molecules. There are two exceptions to this here: of Met 5 undergoes a transition from roughly 30 to -60 (300) in the Cartesian simulation, but remains near 45 in the NEIMO simulation. Secondly, the angle of Gly 2 is rotated from -60 to 60 in the Cartesian simulation, but remains near -60 during the NEIMO calculation. It is possible that fixing the angle terms increases the barriers to rotation enough to prevent these transitions during 5 ps NEIMO simulation at 300 K. The rotational transition of Met 5 did occur after approximately 40 ps of a 50 ps NEIMO simulation using 5 fs timesteps. A 600 K NEIMO simulation using 5 fs timesteps saw both transitions occur by 20 ps, but the temperature was high enough that further transitions continued in both directions over these barriers. It is important to note that the NEIMO formalism explicitly includes the capacity for bond stretches and angle bends between clusters, but the current implementation uses only the dihedral degrees of freedom.

A slightly different type of plot shows the average dihedrals from 10 different NEIMO simulations in Figure . The simulations were identical except for the timestep, which ranged from 1 fs to 10 fs. As explained above, timesteps were chosen to give an integer number of dynamics steps per 0.100 ps. It is clear that the results are extremely consistent for timesteps up to 10 fs. Only the two outer sidechain dihedrals and of Met 5 have significantly different distributions at different timesteps. has for 8, 9, and 10 fs timestep simulations, but for the smaller timesteps. It is possible that the larger timesteps occasionally enable the molecule to jump over rotational energy barriers which cannot be cleared by simulations using smaller timesteps which, in effect, calculate energies and forces at more points along the trajectory.

In order to quantify the dihedral distributions, we represented each distribution by a gaussian, using the average, , and standard deviation, , from the 50 datapoints:

These gaussians were normalized to give a probability function,

which is a normal distribution. Note that the constant in Equation () is derived for nonperiodic variables, and the total probability in Equation () does not exactly equal 1.0 if is so large that the probability is non-zero for every value of . This is not the case for any of the distributions we have analyzed.

Distributions from two different simulations can be compared by calculating the overlap, , of the functions and :

If and are defined as

where and , and and are defined as in Equation () for and , then the product of these functions is also a gaussian:

Here, is defined as usual from , where

The constant in Equation () is

Inserting Equation () into Equation () gives a formula for the overlap:

equals 1 if the two distribution functions are identical and equals 0 if there is no overlap.

The overlaps from the 10 NEIMO simulations plotted in Figure are shown in Figure . Each line represents the overlap between the 1 fs timestep simulation and one of the simulations with a larger timestep. A second figure, Figure , specifically shows the overlaps between the 1 fs simulation and the 2, 5, and 10 fs simulations at a higer resolution. As expected, there is almost 100%overlap among the NEIMO simulations, which indicates clearly that the molecular dynamics are very consistent across a range of timesteps of 1 fs to 10 fs. The only exceptions to these are dihedrals of Met 5 and the of Gly 2. The relatively small overlap of the latter is actually due to the very small value of (0.3) for the 1 fs NEIMO simulation. The discrepancy in the methionine sidechains is also due primarily to differences in rather than for the smaller timestep simulations. At larger timesteps, however, both and differ.

Overlaps between the dihedral distribution from the Cartesian simulation, and those from the NEIMO simulations, are shown in Figure . Here, the overlap is quite small for of Met 5 and the backbone dihedral of Gly 2, as indicated by the large differences in note above. A third very-low overlap is seen for the of Gly 2. This difference is completely hidden in Figure since it is due entirely to the extremely low value of in the NEIMO simulations. The value is so low, in fact, that it does not appear in Figure for the 1 fs NEIMO simulation. The overlaps are greater than 65%for 19 of the 22 dihedrals for every NEIMO timestep. Excluding the methionine residue, overlaps are greater than 90%for 13 of the 16 dihedrals.