9034
. I Phys. . Chem. 1995, 99, 9034-9044
Conformational Flexibility of Manxane Revealed by Adiabatic Mapping, Normal Mode Analysis, and Molecular Dynamics Richard B. Sessions,*Tt David J. Osguthorpe,’ and Pnina Dauber-Osguthorpe’ Molecular Recognition Centre, Medical School, University of Bristol, University Walk, Bristol BS8 ITD, U.K., and Molecular Graphics Unit & School of Chemistry, University of Bath, Claverton Down, Bath BA2 7AY, U.K. Received: October 14, 1994; In Final Form: March 13, 199.5@
We have investigated the relative stability of the conformational minima and the barriers and modes for interconversion of manxane (bicyclo[3.3.3]undecane) using a molecular mechanics approach. The study involved adiabatic mapping of the energy surface as a function of the torsion angles, normal mode analysis, and molecular dynamics simulations followed by digital signal analysis. Energy mapping revealed two C, transition state, in which the conformational minima, C 3 h and C,. The energy surface has a single C 3 h plane of symmetry is maintained. There are two types of C, C, transitions. One of these (11) involves a symmetric transition state, and the other (11’) has two nonsymmetric transition states. The calculated overall c 3 h interconversion agreed with the observed value. Partitioning of the energy of energy barrier for c 3 h the minimum energy conformations and the various transition states revealed the underlying causes for relative stability. The low-frequency modes of the two conformations, as obtained by normal mode (NM) analysis and by extracting characteristic modes of motion from the molecular dynamics (MD) trajectory, were very similar, with some differences due to the inclusion of anhannonicity in the MD treatment. The evolution of modes with time was monitored using sliding frequency distributions (SFD).The results highlighted the importance of energy transfer between modes. Channeling of energy into one of the low-frequency, torsionally active, modes is a prerequisite for the occurrence of a conformational transition. The characteristic features of the active mode in the period preceding a conformational transition can determine the nature of the transition.
-
Introduction Molecular mechanics calculations provide a tool for understanding the structural and energetic preferences of molecules. Using an analytical representation of the energy of a molecule as a function of internal coordinates and interatomic distances, it is possible to reproduce or even predict experimental properties such as molecular structure and conformational energy differences. In addition, these calculations can reveal the underlying factors determining the observed properties by partitioning the energies into components, down to atomic level, thus providing insights unavailable from experiments. Similar interest exists in obtaining detailed information about the flexibility and dynamic behavior of molecules. The theoretical study of molecular motion has involved three approaches to date: normal mode (NM) analysis, adiabatic mapping, and molecular dynamics (MD) simulations. The NM approach has long been used in the interpretation of vibrational spectra.’-3 This analysis partitions the overall motion of the molecule into characteristic motions and provides a “pictorial” description of these motions in the form of normal mode vectors. However, it is only applicable for molecules undergoing small fluctuations around an equilibrium structure. In order to deal with large conformational changes, including conformational transitions, the method of energy (adiabatic) mapping has been used. This involves calculating the energy of the molecule at grid points along the coordinates involved in the conformational changes and reveals aH minimum energy conformations and transition states. Identifying optimal reaction paths connecting different conformations (or reactants and University of Bristol. University of Bath. Abstract published in Advance ACS Absrrucrs, May 1, 1995
0022-365419512099-9034$09.00/0
-
-
products in chemical reactions) is an essential part of transition state theory. Procedures to locate transition states have been developed first for small molecules, in the field of quantum ~ h e m i s t r y . ~Similar -~ ideas have been followed in developing algorithms for locating transition states for conformational changes in larger molecule^.^-^^ These algorithms start from one minimum energy conformation and move “uphill” along the direction of an appropriate normal mode coordinate to a transition state, or drive a set of variables (reaction coordinate) which connect two conformations. In principle, this provides an efficient way of finding transition states without scanning all available degrees of freedom. However, for large systems it may not always be obvious which normal mode to follow or how to choose the reaction coordinate. The third method for studying flexibility and dynamic behavior of molecules is the MD simulation, in which the equations of motion for the atoms in the system are solved n ~ m e r i c a l l y . ’ ~The , ’ ~ resultant trajectory represents a realistic and comprehensive description of molecular motion, taking into account the full flexibility of the molecule, including conformational transitions. Although the laws underlying the motion are quite simple, the overall motion can appear very complex, resembling noise. Fourier transform (FT) techniques have often been used as an aid in understanding the trajectories resulting from MD simulations. For example, for assemblies of atoms or small molecules, the FT of velocities or coordinates has been used to calculate the frequency distribution in the Others have used the FT of the dipole moment,Is the kinetic energy,I9 or specific structural properties of i n t e r e ~ t ~ Oto- reveal ~~ the associated frequencies. We have developed a set of methods for characterizing the various components of motion in MD simulations utilizing digital signal processing techniques. The characteristic frequencies of motion are revealed by Fourier 0 1995 American Chemical Society
Conformational Flexibility of Manxane transforming the atomic coordinates. Subsequently,it is possible to eliminate some modes of motion and focus on motions of interest by using filtering techniques24 (e.g. eliminating highfrequency bond stretching and angle bending and retaining the low-frequency conformational m ~ t i o n * ~ - * This ~ ) . approach has been extended recently to enable the extraction of vectors defining the characteristic motion for each frequency in the MD simulation, analogous to those obtained from normal mode analysis. If only small oscillations about a single conformational minimum occur during the MD simulation, extracting modes from MD simulations leads to results similar to those obtained from normal mode analysis. However, intermixing of modes or even the appearance of new modes may occur in parts of the simulations which involve large conformational changesz8 In order to gain further understanding about the relative importance of the various modes, we introduce here the “sliding frequency distribution” (SDF), which involves monitoring the changes in frequency (and associated modes) during the simulation. Whereas the methods for finding transition states mentioned above require prior knowledge of a suitable normal mode (or reaction coordinate), the identity of the modes which are important for the transitions will emerge from the analysis of the MD simulation. The conformational analysis of cycloalkanes has attracted attention since the early 1 9 6 0 ~ There . ~ ~ has recently been an upsurge of interest in understanding the detailed conformational behavior of common and medium ring cycloalkanes, spurred by current increases in computing power. Examples include the use of adiabatic mapping to explore interconversion barriers in c y c l ~ h e x a n e , a~number ~ ~ ~ ’ of four- and five-membered-ring compounds,32 and a comprehensive study of the four- to 12membered ring cycloalkane series.33 Molecular dynamics has been used to investigate pseudorotation in ~ y c l o p e n t a n eand ~~ a substituted ~ycloheptane.~~ Medium ring bicyclic compounds are an intriguing class, as the large degree of strain energy and restricted conformational mobility can confer highly unusual physical and chemical proper tie^.^^.^^ We have chosen to investigate the conformational behavior of the medium ring bicyclic hydrocarbon manxane (bicyclo[3.3.3]undecane) since its strained geometry38 and inherent chemical symmetry lead to a limited and well-defined number of possible conformations and interconversions. The molecule can exist in just two conformations. The first conformation has C3h symmetry and is characterized by having each trimethylene bridge oriented in the same sense (clockwise or anticlockwise) when viewed along the bridgehead-bridgehead axis. The second conformation has C, symmetry corresponding to one of these bridges being flipped in the opposite sense. The possible conformational transitions and altemative paths for interconversion are shown in Figure 1. The relative orientations of the bridges are conveniently described by the pairs of torsions 4 and ly (Figure l), for each bridge, as these interchange their signs on flipping. The bridge flips also involve a change in the closely coupled pair 8 and t, from -40, -40 to ~ 9 5 -95. , The purpose of this work is to investigate the flexibility of manxane in terms of the relative stability, transition barriers, and dynamic behavior of the possible conformations. Of particular interest is the structural, energetic, and dynamic behavior of the molecule just before, during, and just after a conformational transition.
Methods The molecular mechanics calculations in this study, including minimizations, normal mode analysis, adiabatic mapping, and molecular dynamics, are based on a valence force field (VFF) representation of the potential energy of the molecule. The potential energy is expressed as an analytical function of all
J. Phys. Chem., Vol. 99, No. 22, 1995 9035
Figure 1. Interconversion pathways. The
C3h and C, conformations are shown in black and gray lines, respectively. The three bridges are labeled for one of the C3h and one of the C, conformations. The a and /3 methylenes are identified for bridge A. The main torsion angles, defining the conformation of the molecule, I$ (Cbh - C i - C{ Ct), )I (C; - C{ - C’‘ - Cbr), 6 ( C i - C,, - C; - C{), and T (Cg - C; - Cbh,- C,)2 are shown for bridge B.
the internal coordinates and the interatomic distances.39 The parameters for this force field were derived by fitting experimental properties of model compounds. The observables used in the fitting included equilibrium structures, vibrational frequencies, rotational barriers, energy differences between conformers, sublimation energies, and crystal packing. A significant part of the experimental data (most molecular structures) relate to the minima of the potential energy surface. Other observables (e.g. strained structures and vibrational frequencies) define the shape of the surface, and some (e.g. rotational barriers and energy differences between conformers) give information about higher energy regions. The validity of the force field for the manxane molecule specifically will be evaluated by comparing the corresponding calculated and experimental structural and thermodynamic properties. Minimization and Adiabatic Mapping. The initial structures of the two possible conformations of manxane were built using molecular graphics techniques. The minimum energy conformations were obtained by minimizing the energies of the initial structures with respect to the coordinates. Strict convergence criteria were applied as required prior to NM analysis (the highest derivative of the energy with respect to Cartesian coordinates was smaller than kcal/&. The adiabatic energy maps were calculated using the VFF expression for the potential energy, with two additional terms to force 4, ly to adopt the values 40,ly0, (or 8,t to adopt values 80,TO) while minimizing all other degrees of freedom. The forcing terms were of the form Ef = k(4 - 4 0 ) ~ where k = 500.0 kcal/rad2. The values the torsions were forced to adopt were incremented by steps of lo”, to produce the complete maps. The range over which the torsions were varied was chosen so that it included the minimum energy and conformational transition regions. Normal Mode (NM) Analysis. In NM analysis the Lagrange equations of motion are solved by assuming small oscillations in the vicinity of a local minimum. Under these conditions the potential energy is approximately quadratic (or harmonic). The second-derivatives matrix is diagonalized to yield a set of 3n solutions: q i= 1, cos(2nv,t
+ Ek)
k = 1, 2, 3, ..., 3n
(1)
qi is the mass-weighted Cartesian displacement coordinate, the lik are the normalized eigen vectors that define the normal modes of motion, the Yk = i1’I2/2nare the frequencies, as obtained from the eigen values, A, and Ek are the corresponding phases. (Six of these solutions correspond to rigid body translation and rotation and have zero frequencies.)
9036 J. Phys. Chem., Vol. 99, No. 22, 1995
Sessions et al.
Molecular Dynamics (MD). In MD studies no assumptions are made about harmonicity. Instead, the equations of motion are solved numerically.’* The minimized structure of the C3h conformation was used to provide initial coordinates for the simulation. Random initial velocities consistent with the required temperature were assigned. A leap-frog algorithm with a time step of 1 fs was used. Temperature scaling was used only in the initial equilibration period and tumed off during the rest of the simulation. A principal tool used in the analysis of the MD trajectories was the Fourier transform The FT identifies the frequencies and amplitudes of the different sinusoids that combine to form an arbitrary function of time. For a discrete function of time, h(tm),the FT,H(va),is defined as
TABLE l: Energy and Structure of Manxane
Conformations and Transition Statesa C 3 h (exptl)b
@A *A
4B VB @C *C
@A
SA
OB tB
@c TC
N in= 1
(3)
CQ
cE
CQ CE Ci-Cbh-C: C:-Cbh-C: Ci-Cbh-C: H-Cbh-Ca H-C-H
total
(5)
where H,(va) is the FT of the mass-weight Cartesian displacement coordinates. The relative direction is given by sign(lia) = sign(Ai(va)) for Ai(va)> Bi(va)
(6)
sign(l,) = sign(Bi(va)) for Ai(va)< Bi(va) For further details see ref 28. A dot product, S k K , between vectors representing two modes, k and k‘, was used as an indicator of the similarity between these modes. (7) i
In principle, each independent mode should correspond to a “peak” in the frequency distribution. Fluctuations in conforma-
-54 54 54 -54 0 0
-34 65 65 -34 20 20
42 (36) -42(-36) 42 (36) -42(-36) 42 (36) -42(-36)
95 -95 48 -48 43 -43
71 -71 45 -45 40 -40
98 -98 70 -70 44 -44
99 -99 45 -45 73 -73
78 -109 30 -63 42 -100
122 120 122 122 119 117 119 114 115
125 128 119 118 119 118 116 115 117
119 116 125 128 119 117 116 115 117
121 119 121 119 126 127 121 115 115
152 (150)
Filtering MD Trajectories. The filtering technique is used in order to extract the conformational motion and energetics bond from the overall myriad of motions in the MD ~ i m u l a t i o n . ~ ~ - * ~ angle The technique involves three steps: Fourier transforming (eq torsion nonbond 2 ) , applying a filter, F(v), and inverse Fourier transforming.
11;al = IHj(va)I = [Ai(vJ*+ Bi(va)21’’”/2
-65 65 0 0 65 -65
119(120) 119(118) 119(120) 119(118) 119(120) 119(118) 116(114) 116(114) 116(114) 102 (105) 103 (102)
C-H
Where hi(tm) can be a structural property (e.g. Cartesian or intemal coordinate) or an energy component. Extracting Characteristic Modes. As discussed above (eq l), a characteristic mode is described by a typical frequency and a set of atomic amplitudes for all atoms. An analogous description can be obtained from MD simulations using Fourier transform techniques. A sample mode, a,is defined by the set of amplitudes, lis, of all coordinates at a specific frequency, va. The absolute value of li, is given by
62(67) -62(-67) 62 (67) -62(-67) 62 (67) -62(-67)
H-Cbh-Ca-Co
c-c
a
TS II’b
Ca
c!
where A and B are the real and imaginary parts of the FT. Frequency Distribution. The frequency distribution function, g(v),gives the number of characteristic motions with a frequency v. Using the fluctuations in the mass-weighted Cartesian displacement coordinates, qi(tm),as the time function (hi(tm))in the FT definition (eq 2 ) , the frequency distribution is given by26
TS I1 TS II’a
Structure -53 0 53 0 50 63 -50 -63 64 65 -64 -65
TS I
C,
121, 119 119 119, 121 119 124 125 119 113, 116 116, 113
1.55 (1.54) 1.11 (1.09)
2.8 20.4 6.9 13.9 44.0
Energies 3.1 25.6 7.1 13.9 49.7
2.6 29.7 10.8 12.9 56.0
2.8 29.2 10.3 13.2 55.5
2.4 36.0 11.3 12.3 62.0
3.2 28.2 13.7 15.9 61.0
a Angles are in degrees, distances in angstroms, and energies in kcaU mol. The electron diffraction experimental structure gives average C-C and C-H bond lengths, H-C-H and H-Cbh-CQ bond angles, and H-Cbh-C,-Cp torsion angles, which are compared with the corresponding calculated values of the C 3 h conformer. Coordinates of the calculated structures are available as supplementary material.
tion during the molecular dynamics and the discrete and finite nature of the trajectories cause line broadening in the frequency distribution. Thus, in order to obtain the overall contribution to a characteristic motion, we sum the set of “sample modes” obtained for the discrete frequencies in the frequency range Av which have a large Skk‘. The set of characteristic modes will not always be identical to the set of modes obtained by normal mode analysis. Since the simulation includes anharmonic effects, the extracted modes may correspond to intermixing of individual normal modes, or even reveal additional modes, which may occur due to conformational events such as intramolecular free bond rotations.28 Software. The minimizations NM analysis and MD simulations were carried out with DISCOVER.41 Correlation matrices were calculated and plotted using CMPMOD,42 and MD trajectories were analyzed with the package FOCUS.43.44
Results Structure. Only one of the two possible manxane conformations, the C3h, has been observed experimentally. The structure of manxane itself has been determined by gas phase electron d i f f r a ~ t i o n .A ~ ~comparison of the observed and calculated C3/r structure of manxane is given in Table 1. Torsion angles lie within 6” and valence angles within 2”. Both the observed and calculated structures indicate that the molecule is strained. In particular, the valence angles deviate from tetrahedral values
J. Phys. Chem., Vol. 99, No. 22, 1995 9037
Conformational Flexibility of Manxane by 5-11'. The X-ray structures of two manxane derivatives are also available: the isoelectronic derivative, manxine hyd r ~ c h l o r i d e ,and ~ ~ the 1,5-dihydroxy The rms deviation of the Cartesian coordinates of the ring heavy atoms between the calculated structure and both X-ray structures is 0.05 A, which is similar to the rms deviation between the heavy atoms of both X-ray structures of 0.03 A. Thus, the calculations reproduce the available structural information well. The fact that the force field was parametrized using a wide range of cyclic and acyclic hydrocarbon^:^ and in particular that it reproduces this strained manxane structure, gives confidence in the force field in regions away from conformational minima where intemals are necessarily away from minimum energy geometries. Energy Surfaces. The energy surfaces corresponding to the three types of possible transitions I, 11, and 11' (defined in Figure 1) are shown in Figure 2a(i)-c(i), respectively, in terms of 4 and q. The corresponding energy surfaces in terms of the torsions 8 and t are shown in Figure 2a(ii)-c(ii). Each panel shows two minima interconnected by one or two saddle points corresponding to transition states in the conformational interconversions; these are characterized by having a single imaginary frequency. The energies of the minima and the barriers to interconversion are summarized in Table 1. Contributions to the total energy from the different components of the force field equation39are included. In Figure 2a(i) and (ii) there are two local minima, corresponding to the C 3 h and C, conformations. The energy difference between the two minimized conformations is 5.7 kcal/ mol. Qualitatively, the C 3 h c, transition involves motion of a single bridge which can be described as a reflection of the /3 methylene in a plane defined by the two bridgehead atoms and C, transition the two a carbons of that bridge. If the C 3 h involved only such a reflection of the3!, methylene group (A), a clash (1.44 A) would result between the p hydrogen atoms on the pair of bridges which point toward each other (A and B). This clash would manifest itself in a high nonbonded repulsion contribution to the overall energy. The extent of this clash is relieved by adjustments to the torsion and valence angles. The 4, q torsions of bridges A and B of the C, conformation deviate from a staggered conformation by 7-10'. These deviations increase the torsional energy only slightly, since the corresponding torsional potential is not very steep, particularly at the minimum energy region. In addition, the valence angles in the eight-membered ring containing bridges A and B of this conformation have opened up by x3', at a cost of %5 kcal/mol. The combined effect of the changes in torsion and valence angle was to flatten this eight-membered ring, and the actual distance between the two p hydrogens in the minimized C, conformation was increased to 2.1 A. The increase in valence and torsion angle energy is more than compensated for by the relief of the repulsion energy. The energy surfaces, as function of 4, q or 8, t,are shallower around the C, minimum than around the C 3 h minimum. This is due to the fact that changes in torsion angles in the C 3 h conformation increase both angle strain and nonbonded energies, whereas in the C, conformation nonbond energies are counteracted by angle strain; thus the overall energy hardly changes. The transition state for the C 3 h C, conversion, TS-I, has one bridge in a planar conformation. This involves two torsions in an eclipsed conformation and opening up of the valence angles cbhcacp and CaCpCa from 119' to 125' and 128', respectively. Thus, the main contribution to the 12.0 kcal/mol transition barrier comes from torsion ( ~kcdmol) 4 and valence angle (%9 kcdmol) strain. The nonbond interactions are more favorable ( ~ kcdmol) 1 in the transition state than in the
-
-
-
TABLE 2: Low-Frequency Normal Modes of Mamane
c3h
CS
sym
no.
NM(MD) NM(MD) frea sym no. frea
A"
1
208(208)
A"
1
73(78)
E'
2+3
259(241)
E"
4+5
271(267)
A'
6
348(358)
A' A' A" A" A'
2 5 3 4 6
214(202) 301(313) 239(248) 266(261) 338(326)
A'
7
387(339)
A'
7
379(397)
descriDtionofmotion ~~
(i) twist, around the bridgeheads axis, a methylene motion (ii) wag, p methylene motion (iii) rocking, motion of two bridges (iv) breathing, a-p-a valence angle bends (v) wag, a a n d p methylenes
a Frequencies are in cm-]. Normal mode frequencies of the minima are available as supplementary material.
minimum energy conformation due to the relief of the close interactions of the Ha-HB atoms of adjacent bridges. In Figure 2, parts b(i),(ii) and c(i),(ii), both energy minima correspond to C, conformations. The transition state for the C, C, conversion type 11, TS-11, (Figure 2b) is very similar to the transition state of the C 3 h C, conversion, TS-I. Both have one planar bridge, and the other two bridges are not pointing to each other. Thus, the total energy and the energy components of these two transition states are very similar. (See Table 1.) On the other hand, transition type 11' (Figure 2c) involves flipping the third bridge while the two other bridges are still pointing to each other. Thus, on top of the angle strain of the eight-membered ring containing these two bridges, extra strain is required to flatten the third bridge. Therefore, a symmetric transition state, with 4 = = 0, requires %15 kcal/mol of valence angle strain (TS-II'a in Table 1). As can be seen in Figure 2, the lowest energy path for transitions I and 11maintains symmetry with respect to the mirror plane passing through the three p carbons, i.e. 4 = q = 0 at the transition state. However, the surface corresponding to transition type 11' has a double saddle located off the diagonal connecting the two minima; hence, the lowest energy pathway involves bridgehead-bridgehead twisting. This twisted transition state, TS-II'b, has higher torsion and nonbonded strain than the symmetric transition state, TS-II'a, but significantly lower valence angle strain energy. The asymmetric transition state is lower in energy than the symmetric one by about 1 kcal/mol and is characterized by having one imaginary frequency (whereas the symmetric TS has two). Thus, the lowest energy barriers for a C, C, interconversion are M . 8 and 11.3 k c d mol for type I1 and type 11', respectively. The rate-determining step for the complete c 3 h C3h interconversion via path 11' is the second step, whereas for path 11all the transition states have similar energies. Thus, the lowest energy C3h c 3 h interconversion involves path 11and an overall energy barrier of 12.0 kcal/mol, which is in agreement with experiment ( ~ 1 kcal/mol). 1 (See ref 48 and references therein). The agreement between the observed and calculated conversion energy lends additional support to the validity of the force field for this molecule in regions far from the conformational minima. Normal Mode Analysis. In order to investigate the types of motion the molecule undergoes, we have carried out normal mode analysis of the C 3 h and the C, conformations. The lowfrequency characteristic modes of motion for the conformations are described in Table 2 and shown in Figures 3a and 4a. The common characteristic of these modes is that they describe conformational motion; that is, they involve mainly oscillations in torsion angles. Modes of symmetry A' and E' are symmetric with respect to the mirror plane passing through the three p carbons, while modes of symmetry A" and E" are asymmetric with respect to this plane. Examination of the main torsional
-
-
-
-
-
9038 J. Phys. Chem., Vol. 99, No. 22, 1995 60
0
-60
Sessions et al.
-120
-60
TAU
(ii) 60
0
-60
PHI /L\
0
-hO
60
v) 0
-120
-60
w
3:
0
0
-120
-60
TAU
(ii)
(1) 0
-60
60
PHI 0
-611
60
H v)
n.
-120
-60
TAU
(ii) 0
b0
60
PHI
Figure 2. Adiabatic energy surfaces representing the three possible conformational transitions (contour level 1 kcaYmol). (i) The energy as a
-
torsions (deg). (ii) The energy as a function of the 6' and t torsions (deg). Selected sections of the filtered torsion function of the C$ and trajectories corresponding to periods with conformational transitions are superimposed on the respective energy maps. (a) C 3 h C, energy surface, type I transition (9-1 1.8 ps). (b) C, C, energy surface, type I1 transition (20-24 ps). (c) C, C, energy surface, type 11' transition (12-15 PSI.
-
-
J. Phys. Chem., Vol. 99, No. 22, 1995 9039
Conformational Flexibility of Manxane
CD
1
1
208
2&3
/ 259
4
a a
2
24 I
3
267
307
339
Figure 3. Comparison of the normal modes of the C 3 h conformer with modes extracted from the 5-10 ps section of the trajectory. (Two views are given for each mode: down the bridgehead-bridgehead axis and perpendicular to this axis.) (a) Normal mode coordinate displacement vectors. (b) Extracted mode coordinate displacement vectors.
4
5
6
A@
7 379
components of these modes shows that modes of symmetry A‘ and E’ correspond to motion along the diagonal connecting the two minimum energy conformations, while modes of symmetry A” and E” correspond to motion perpendicular to this diagonal. The lowest frequency mode, type (i) symmetry A”, is a twisting mode in which the a methylenes on each side of the mirror plane are rotating in opposite directions. For the CY,conformation, this type of motion includes 4, v distortions along the steep energy gradient (perpendicular to the interconformation diagonal). However, all bridges move in-phase, and thus the clashes between bridges are maintained at a minimum and the overall vibration is of low energy. The corresponding motion for the C, conformation is much lower in frequency. Only bridges A and C have significant motion for this mode. As seen from the maps, the regions of the C, energy minima are overall a lot shallower then the C3h minimum energy conformation. In addition, motion along the off-diagonal direction for bridges A and C is associated with particularly shallow energy gradients. (Top right-hand comer of Figure 2, parts a and c). As discussed in the previous section, bridgehead-bridgehead twisting of this conformation hardly changes the total energy, and thus motion in this direction requires very little energy. The next to lowest mode of motion, type (ii) symmetry E’ or A’, can be described as a “wag”. It mainly involves motion of the /3 methylene groups. In terms of changes in the main torsional angles, the mode describes motion in a direction along the interconformation diagonal on the energy maps. Unlike the twist mode, discussed above, the amplitude and the direction of rotational motion of the three bridges are not equal. The third mode, type (iii) symmetry E” or A“, is a “rock”. This mode involves a rocking motion of two bridges about an axis passing through the center of the molecule and the /3 position of the third bridge. This mode, like the first, involves motion in a direction perpendicular to the interconformation diagonal. The next vibration, type (iv), involves a breathing type of motion along the bridgehead axis which mainly consists of
397
Figure 4. Comparison of the normal modes of the C, conformer with modes extracted from the 25.0-30.0 ps section of the trajectory. (Two views are given for each mode: down the bridgehead-bridgehead axis and perpendicular to this axis.) (a) Normal mode coordinate displacement vectors. (b) Extracted mode coordinate displacement vectors.
angle bending rather than torsional motion. The next mode, type (v), is a wag which involves in-phase motion of the a and /3 methylenes. This mode also has relatively large components of angle bending. Thus, conformationally significant modes will not be seen at higher frequencies than those mentioned already. Molecular Dynamics Trajectories. While normal mode analysis can give valuable information about molecular motion in the vicinity of a minimum, the dynamic behavior of the molecule in other regions (including conformational transitions) has to be investigated by other methods. Classical (Newtonian) dynamics simulations provide a suitable probe. A number of short MD simulations at different temperatures and initial random velocities were carried out. Several conformational transitions occurred in a 20 ps trajectory at 600 K, so this trajectory was extended to 2 ns. As discussed above, we expect the conformational modes of motion to be of frequencies lower than 300 cm-I. Thus, in order to focus on conformational motion the 4, v and 8, z torsion trajectories were filtered to retain only frequencies below 300 cm-I. The changes in conformation, as described by the pairs of filtered torsion angles 4, and 8, z are monitored in parts a and b of Figure 5, respectively. The first 50 ps of the trajectory is characterized by a large number of conformational transitions. This reflects the fact that the random initial velocities did not correspond to equipartition of energy into the various modes, with torsional modes getting a higher proportion of the energy. The excess energy in conformational motion at this period manifested itself also in the ability to overcome higher barriers. All the examples of the high energy C, C, $J
-
9040 J. Phys. Chem., Vol. 99, No. 22, 1995
bridge
Sessions et al.
A
bridge B
A
bridge C
,
-
.
.
.
.
I
.
.
_
500
0
,
-
-
-
.
1000
I
-
1500
-
-
2000
time ( p s ) bridge
J
-
A
500
0
1000
1500
2000
time (ps)
3, j
potential
-
kinetic l
0
'
"
~
l
500
'
"
'
I
1000
"
"
1
1500
2000
time ( p s ) Figure 5. Conformation and energetics during the MD simulation. Torsions are in degrees, energies in kcal/mol, and time in picoseconds (ps). (a) Filtered trajectories of the 4 (black) and (gray) torsion angles of each bridge. (b) Filtered trajectories of the 8 (black) and t (gray) torsion angles of each bridge. (c) Potential, kinetic, and total energy trajectories. The original and filtered energies are shown in gray and black, respectively. (Temperature scaling was carried out during the first 50 ps.)
type 11' transition (three) occur in this period, whereas all other C., C.,transitions in the simulation are of the lower energy type 11. A full qh-. qh(clockwise to anticlockwise) interconversion occurs at 39.7 ps. The reverse conversion, Gh-. qh,back to the initial conformation occurs at 45.4 ps. No further transitions occur until nearly 1 ns has elapsed, whereupon sequential flips of bridges A, B, and C (921-923 ps) rapidly lead to a Gh qi interconversion. Soon after
-
-
(1028.5 ps) bridge A flips, rapidly followed by bridge C, and the reverse sequence occurs after about 8 ps to regenerate the qh conformation (1036.6 ps). A further, more protracted, set of events occur next (1559- 1580 ps) causing a Cz, -. Gh conversion. The final event in the trajectory is a flip and rapid flip back (1762-1764 ps) of bridge A. The underlying link between the energy surface of the molecule and the trajectory it follows during the MD simulation
Conformational Flexibility of Manxane
J. Phys. Chem., Vol. 99, No. 22, I995 9041
is demonstrated in Figure 2. The filtered trajectories of the 4, q and 8, z torsions, of selected time periods in which conformational transitions occurred, were overlaid on the corresponding energy surfaces. As can be seen from these figures, during most of the time the molecule undergoes torsional fluctuations in the regions within -4-5 kcavmol of the local minima. The conformational transitions follow fairly closely the path of lowest energy between the the two minima. In the case of transitions type I and II (Figure 2a,b) this path is very close to the diagonal line connecting the two minima; that is, the symmetry with respect to the mirror plane is retained. For the type 11’ transition, an off-diagonal route is taken, which passes close to one of the two asymmetric transition states. Using the filtering technique, it is also possible to reveal the relation between the energies obtained by minimizations or adiabatic mapping and those calculated from MD simulations. The original energy trajectories obtained during the simulations (Figure 5c) include large contributions from high-frequency bond and valence angle fluctuations which obscure any changes in conformational energy. The energy trajectories were therefore filtered to retain only frequencies lower than 5 cm-I. Figure 5c depicts, in addition to the original energy trajectories, the filtered trajectories of the potential, kinetic, and total energy, which are associated with the changes in conformation during the simulation. There is no drift in the total energy of the system over the 2 ns, indicating the stability of the numerical integration. The underlying minimum energies of the two conformations can be obtained from the filtered trajectory using the relation
where P and K are the filtered potential and kinetic energies, Pminis the minimized potential energy, and PeXCeSS and Kexcess are the excess conformational potential and kinetic energy. Assuming the overall excess energy is equally partitioned between potential and kinetic energy, the minimum energy is defined by
Pmin= P - K The filtered potential and kinetic energies for the periods corresponding to the C, conformation are 101 and 51 kcal/mol, respectively, resulting in a minimum energy of 50 kcavmol. The corresponding filtered potential, kinetic, and minimum energy for the C3h conformation are 98, 54, and 44 kcavmol, respectively, in good agreement with the values obtained by minimization or mapping (see Table l),justifying the assumption of equipartitioning of energy into kinetic and potential terms. The transition barriers can also be estimated from filtered energy trajectories. For this purpose we need to retain the fluctuations due to torsional modes which can lead to transitions. As discussed before, the important modes are of frequencies below 300 cm-I. The corresponding energy fluctuations will be in the range up to 600 ~ m - ’ . ’ ~The barrier to the C 3 h C, transition is given by the energy difference between the potential energy peak at the transition point (4 and q change signs) and the lowest potential energy in the period corresponding to the C 3 h conformation. This barrier is thus calculated to be ~ 1 0 . 5 kcavmol. Extracted Modes. The nature of the oscillations for each of the conformations can be elucidated by applying digital signal techniques to the corresponding periods during the dynamics. These techniques include extracting sets of vectors which represent the characteristic modes of motion (eqs 5, 6) and filtering the trajectories to generate trajectories which include only selected types of motion, to be monitored on a graphics
-
display system (eq 4). During the time periods 5-10 and 2530 ps of the simulation the molecule oscillates in a single conformation (C3h and C,, respectively) with no transitions (Figure 5). The typical low-frequency motions occurring during each of these periods were examined by extracting the sets of sample modes with frequencies up to 500 cm-I. The correlation between each pair of sample modes is shown in the bottom panel of parts a and b of Figure 6 for the c 3 h and C, conformations, respectively. The correlation between these modes and the modes obtained from the NM analysis is depicted in the top panel. The corresponding frequency distributions are given above the correlation matrices and indicate which modes are active. The high correlations along the diagonal in the autocorrelation matrices indicate that there are bands of frequencies each with a typical motion. The high correlation along the diagonal in the cross correlation matrices indicates that the same types of motion, with similar frequency ranking, occur in the MD simulation and in the NM analysis. Groups of sample modes in a frequency band which have a high autocorrelation were added up to yield a set of characteristic modes for this time period. These extracted modes are shown in Figures 3b and 4b, along with the corresponding NMs (Figure 3a and 4a). There is considerable similarity between the normal and extracted modes, although there are some variations in frequency (see, for example, modes 6 and 7 C 3 h ) and in motion. These differences are due to the fact that during an MD simulation the local environment of each atom is constantly changing, and thus the “effective” force on the atom is different from the corresponding force at the minimum energy conformation which is used to calculate the NM. The MD simulation takes explicit account of anharmonicity, whereas the NM analysis uses the harmonic approximation around the particular minimum. Evolution of Modes with Time. The NM description of molecular motion implies equal and constant distribution of energy into the different characteristic modes. MD simulations provide a more realistic description of energy distribution among modes of motion; while equal distribution occurs over the long term, short-term fluctuations occur. As demonstrated for trajectories of chemical reactions of small molecules, the energy distribution between normal mode coordinates changes with time.49 Energy transfer between modes plays a major role in initiating conformational transitions as well. We monitored the evolution of the low-frequency modes during the MD trajectory by calculating a sliding frequency distribution (SFD) (Le. a set of g(v)’s corresponding to time windows of 5 ps and a starting time which is shifted by 1 ps at a time). Figure 7 shows SFDs during periods of the trajectory when conformational transitions occurred; frequencies and characteristic motions were observed similar to those of the extracted modes in Figure 3b. Inspection of these SFD plots shows that before each c 3 h C, transition there is a growth in intensity of one or more of the lowfrequency torsionally active modes. There are examples of each of the three lowest frequency modes, the twist, wag, and rock, leading to a bridge flip. Thus, it is the rocking mode (iii) which grows before the 11.6 ps transition, and the twist mode (i) grows before the 921.5 ps transition. Both the wag and rock modes (ii and iii) grow before the 1028.5 ps flip, and the wag mode (ii) grows before the 1559.5 and 1762.7 ps events. Hence, each of the lowest three modes can lead to a conformational transition provided enough energy is channeled into it. However, most commonly transitions will be initiated by the wagging mode (by itself or in combination with another mode), since this mode involves torsional changes along the lowest energy pathway between the C3h and C, conformers, as discussed above. The mode which is active just before the transition also determines which one of the three bridges will flip. The mode which grows
-
Sessions et al.
9042 J. Phys. Chem., Vol. 99, No. 22, 1995
(a)
(b)
Figure 6. Autocorrelations between sample modes extracted from the MD simulation (bottom panel) and cross correlations between sample modes extracted from the MD and modes obtained from NM analysis (top panel). (See eq 7.) Correlations larger than 0.5 are shown in black, and correlations between 0.3-0.5 are shown in gray. Correlations of modes of the C.V,and of the C.rconformations are shown in parts a and b, respectively. The corresponding frequency distributions, g(v), are shown above the correlation matrices.
in intensity during the 5 ps preceding each of the conformational transitions was extracted and analyzed. In each case the bridge which flips is the one which displays the largest motion along the interconformation diagonal in the adiabatic map. The C,, C.,transitions could not be analyzed in similar detail since there are very few examples in which a C., conformation exists for a long enough period before a transition to generate SFD plots. However, some insight into the link between the active modes of motion and the types of transitions occurring can be gained by comparing g(v) for the time period before the type 11’ transition at 13.4ps and the time period before the type I1 transition at 21.0 ps. The most active low-frequency mode before the transition type 11’ is the twist mode (i). As mentioned above, this mode involves torsional changes which are asymmetric with respect to the mirror plane, and thus it induces a transition which has an asymmetric nature. In the period before the type I1 transition the wag mode (ii), which had negligible intensity before, gained significant intensity. This mode involves motion along the lowest energy pathway and thus leads to a transition of a symmetric nature. An examination of the two lowest frequency modes (Figure 4) reveals that the relative amplitudes of motion of the-three bridges is very different. The twisting mode (i) has large (nearly equal) amplitudes for bridges
-
A and C and a very small amplitude for bridge B, whereas the wagging mode (ii) has the largest amplitude of motion for bridge B, somewhat smaller for bridge A, and very small for bridge C. As shown in Figure 1, transitions type I1 and 11’ involve flips of bridges B and C, respectively (a flip of bridge A would C3h transition). Thus, the wag mode will be lead to a C, effective in inducing transitions of type 11, whereas the twist mode will induce transitions of type 11’. The combined application of SFD and characteristic mode extraction techniques provides a systematic and general method for investigating energy transfer in MD simulations. As we have shown, the frequency of modes which are enhanced prior to a conformational transition is apparent from the SFDs. By extracting the modes corresponding to these frequencies from the appropriate time period section of the trajectory, a precise description of the mode, in terms of atomic amplitudes, is obtained as well. It should be noted that although in this study the modes near the equilibrium conformationswere known from NM analysis, this knowledge is not a prerequisite for characterizing the importance of modes during the MD simulation. Thus, it is possible to apply this analysis to larger and more complex systems for which NM analysis is either not appropriate due to
-
(b)
(c)
925-930
-
924-929
1564 1569
&
13-18
923-928
1 2 - 17n -,-
922-927
-
1562- 1567
11-16-
921-926
-
1561-1566
920-925 919-924 918-923 917-922 916-921 5-10
2
4-9 3-8
1-6
915-920
~ u
-A
-h
1760-1765
-,-,-h/ 1\ 757-1762
1554-1559
912-917
a
1552-1557
:-yo
Ah 1 7 6 1- 1 7 66
1557-1562
1555-1560
1551-1556
4
h
~-,-hn-~,/1 7 5 8 - 1 7 6 3
1553-1558
911-916
--hn----
, . 1 7A 63 - 1 7 6 8
1 5 5 8 - 1 5 63
1556-1561
k
1764- 1769
1759-1764
A
p u
-,-.-.---.,-&
1559-1564
913-918
910-915
0-5
A -
1765-1770
1563 1568
1560-1565
914-919
2-7
(dl
1565.1570 , -
b
rz
h
h 1
h
1550-1555
1756-1761 1755-1760 1754-1759 1753-1758 1752-1757
-r
A~
A
h .
h 2
n I A
1751-1756
A
J
1750-1755
frequency
frequency
frequency
Figure 7. Sliding frequency distributions (SFD) corresponding to time periods in which conformational transitions occurred: (a) 0.0-20.0 ps; (b) 910.0-930.0 ps; (c) 1550.0- 1570.0 ps; (d) 1750.0- 1770.0 ps. Frequencies are in cm-I. The frequency of the mode that is being enhanced for each SFD set is indicated by an arrow.
9044 J. Phys. Chem., Vol. 99, No. 22, 1995
Sessions et al.
anharmonicity of the system or too prohibitive in terms of computer resources.
Conclusions We have used an array of molecular mechanics techniques to investigate the structural, energetic, and dynamic properties of manxane. The minimum energy structure and the adiabatic mapping of the energy surface reproduce well the available experimental information. In addition, the ability to partition the energy into its components provided insight into the origins of the relative stability of specific minimum energy conformations and transition states. For example, it turns out that the main source of the higher energy of the C,conformation, relative to the C3h conformation, is the valence angle strain, while the nonbond and torsion energies are nearly identical. Similarly, it was revealed that the major component in determining the barrier height for transitions is the strain involved in forcing one of the bridges to a planar geometry. For transitions of type I and I1 this strain is counteracted by the absence of two bridges pointing to each other. The higher barrier of transition 11’ is accounted for by the existence of two bridges pointing to each other. Partial relief of strain is obtained for this transition state by bridgehead-bridgehead twisting. The dynamic properties of manxane were investigated by NM analysis and MD simulations. Both methods proved adequate in describing the oscillations around an equilibrium conformation. A comparison of the low-frequency motions of the molecule by NM analysis and motions from the MD trajectory, analyzed by mode extraction, gave similar results but with differences due to the inclusion of anharmonicity in the MD treatment. However, the MD simulations, followed by digital signal analysis, provided further insight with regard to the links between molecular oscillations and conformational transitions. It revealed that energy transfer between modes is crucial for initiating a conformational transition. Channeling of energy into one of the low-frequency, torsionally active, modes is a prerequisite for the occurrence of a conformational transition. Activation of one of the three lowest frequency modes of the C3h conformer can lead to conversion to C,. The wag mode is particularly efficient in inducing a transition since it involves torsional fluctuations along the interconformation diagonal on the energy surface. Similarly, the nature of the active mode in the C, conformation determines the type of a subsequent C, C, transition. The twisting mode induces an asymmetric (type 11’) transition, whereas the wagging mode induces a symmetric (type 11) transition. The difference in potential energy between the two manxane conformers in the MD trajectory was revealed by filtering out high-frequency components and agreed with the results from adiabatic mapping. Thus, by applying the filtering and mode extraction techniques to MD simulations, it was possible to characterize the motion of the molecule at minimum energy regions as well as point out features leading to conformational transitions.
+
Supplementary Material Available: Minima and transition states: coordinates, pictures, and normal mode frequencies (1 1 pages). Ordering information is given on any current masthead page. References and Notes (1) Herzberg, G. Molecular Spectra and Molecular Structure 11. Infrared and Raman Spectra of Polyatomic Molecules; D. van Nostrand Co. New York, 1945. (2) Wilson. E. B.: Decius, J. C.: Cross, P. C. Molecular Vibrations; McGraw Hill: Ne% York, 1955. (3) Hehre, W. J.; Radom, L.; Schleyer, P. v. R.; Pople, J. A. Ab Initio Molecular Orbital Theory: J. Wiley & Sons: New York, 1986.
(4) McIver, J. W.; Komornicki, A. J. Am. Chem. Soc. 1972, 94,26252633. ( 5 ) McIver, J. W.; Komornicki, A. J . Am. Chem. Soc. 1973.95,45124517. (6) Halgren, T. A.; Lipscomb, W. N. Chem. Phys. Lett. 1977,49,225232. (7) Crippen, G. M.; Scheraga, H. A. Arch. Biochem. Biophys. 1971, 144, 462. (8) Nguyen, D. T.; Case, D. A. J . Phys. Chem. 1985,89,4020-4026. (9) Czerminski, R.; Elber, R. J . Chem. Phys. 1990, 92, 5580-5601. (10) Wales, D. J. J . Chem. Soc. Faraday Trans. 1993, 89, 1305-1313. (11) Smart, 0. S. Chem. Phys. Lett. 1994, 222, 503-512. (12) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon: Oxford, 1987. (13) McCammon, J. A,; Harvey, S. C. Dynamics of Proteins and Nucleic Acids; Cambridge University Press: Cambridge, 1987. (14) Kim, E.; Lee, Y. H. Phys. Rev. B: Condens. Matter 1994,49, 17431749. (15) Edvardsson, S.; Ojamae, L.; Thomas. J. 0. J . Phys. Condens. Matter 1994, 6, 1319-1332. (16) Ezra, G. S.; Martens, C. C.; Fried, L. E. J . Phys. Chem. 1987, 91, 3721 -3730. (17) Fredkin, D. R.; Komomicki, A.; White. S. R.; Wilson, K. R. J . Chem. Phys. 1983, 78, 7077-7092. (18) Tiller, A. R. Macromolecules 1992, 25. 4605-4611. (19) Rauhut, G.; Clark, T. J . Chem. Soc., Faraday Trans. 1994, 90, 1783-1788. (20) McCammon, J. A.; Gelin, B. R.; Karplus, M. Nature 1977, 267, 585. (21) Cardini, G.; Schettino, V. Chem. Phys. 1990, 146, 147-153. (22) Axelsen, P. H.; Prendergast, F. G. Biophys. J. 1989, 56, 43-66. (23) Aaqvist, J.; Leijonmarck, M.; Tapia, 0. Eur. Biophys. J. 1989. 16, 327-339. (24) Oppenheim, A. V.; Schafer, R. W. Digital Signal Procesring Prentice-Hall Intemational: London, 1975. (25) Sessions, R. B.; Dauber-Osguthorpe, P.; Osguthorpe, D. J. J. Mol. B i d . 1989, 210, 617-634. (26) Dauber-Osguthorpe, P.; Osguthorpe, D. J. J . Am. Chem. Soc. 1990, 112, 7921-7935. (27) Dauber-Osguthorpe, P.; Osguthorpe, D. J. Biochemistry 1990, 29, 8223-8228. (28) Dauber-Osguthorpe, P.; Osguthorpe, D. J. J . Comput. Chem. 1993, 14, 1259-1271. (29) Hendrickson, J. B. J . Am. Chem. Soc. 1961, 83, 4537-4547. (30) Dixon, D. A.; Komornicki, A. J . Phys. Chem. 1990, 94, 56305636. (31) Laane, I.; Choo, J. J . Am. Chem. Soc. 1994, 116, 3889-3891. (32) Rosas, R. L.; Cooper, C.; Laane, J. J . Phys. Chem. 1990,94, 18301836. (33) Kolossvary, I.; Guida, W. C. J . Am. Chem. Soc. 1993. 115, 210721 19. (34) Cui, W. L.; Li, F. B.; Allinger, N. L. J . Am. Chem. SOC. 1993,115, 2943-2951. (35) Li, F. B.; Cui, W. L.; Allinger, N. L. J . Comput. Chem. 1994, 15, 769-781. (36) Alder, R. W. Tetrahedron 1990, 46, 683-713. (37) Alder, R. W.; Heilbronner, E. E.; Honegger, E.; McEwen, A. B.; Moss, R. E.; Olefirowicz, E.; Petillo, P. A.; Sessions, R. B.; Weisman, G. R.; White, J. M.; Yang, Z.-2. J . Am. Chem. Soc. 1993, 115, 6580-6591. (38) Gundersen, G.; Murray-Rust, P.; Rankin, D. W. H.; Seip, R.; Watt. C. I. F. Acta Chem. Scand. 1983, A37, 823-831. (39) Dauber-Osguthorpe, P.; Roberts, V. A,; Osguthorpe, D. J.; Wolff, J.; Genest, M.; Hagler, A. T. Proteins: Struct. Funct. Genet. 1988, 4 , 3147. (40) Brigham, E. 0. The Fast Fourier Transform and its Applications; Prentice-Hall: Englewood Cliffs, 1988. (41) DISCOVER, Biosym Technologies: San-Diego, CA. (42) Dauber-Osguthorpe, P.; Osguthorpe, D. J. CMPMOD, Molecular Modelling Unit; University of Bath: Bath, U.K., 1991. (43) Osguthorpe, D. J.; Dauber-Osguthorpe, P. J . Mol. Graph. 1992, 10, 178-184. (44) Distributed by Oxford Molecular Ltd., Magdalen Centre, Oxford Science Park, Sanford-on-Thames, Oxford OX4 4GA, U.K. (45) Wang, A. H.-J.; Missavage, R. J.; Bym, S. R.; Paul, I. C. J . Am. Chem. Soc. 1972, 94. 7100. (46) Murray-Rust, P.; Murray-Rust. J.; Watt, C. I. F. Tetrahedron 1980, 36, 2799-2806. (47) Ermer, 0.;Lifson, S. J . Am. Chem. Soc. 1973, 95, 4121. (48) Olah, G. A,; Liang, G.; Schleyer, P. v. R.; Parker, W.; Watt, C. I. F. J . Am. Chem. Soc. 1977, 99, 966-968. (49) McDonald, J. D.; Marcus, R. A. J . Chem. Phys. 1976, 65, 2180-2192. JP942757G
