Introduction


The conformational energy of a molecule with respect to molecular mechanics is defined as the sum of the bond stretching, bond angle bending, harmonic dihedral bending, sinusoidal dihedral torsions and both van der Walls and Coulombic non-bonded interactions. These individual terms form the basis of an empirical potential energy function which can be minimized with respect to the atomic coordinates of the molecule in order to determine stable conformations. Equilibrium in a biochemical system is governed by the Gibbs-Helmholz equation and the tendency of the system to move toward states of lower free energy:

DG = DE + PDV - TDS

In a well-ordered system, DE may dominate this equation and the stable conformations of the molecule may be determined by minimizing internal energy, which is the sum of the potential plus kinetic energies. The kinetic energy however, is a function of temperature only. Therefore, minimization of the potential energy will determine stable molecular conformations if conformational entropy can be ignored.

The techniques of molecular mechanics have been applied recently to modeling enzyme-substrate complexes, particularly where no crystal structure is available and direct experimentation is not feasible. In one recent application, molecular mechanics minimization has been used to obtain intermolecular geometrys for use in subsequent ab initio MO calculations on the enzyme active site.(1) These theoretical investigations are potentially useful for reinforcing proposed reaction mechanisms for enzyme catalysis which have not been experimentally proven.

There are limitations in the technique however. While the mathematical form of the empirical potential energy function varies little from one force field to the next, it is a simplification for the sake of computational speed. For example, the electronic charge distribution of the molecule is modeled as simple atom centered point charges. This neglects the anisotropic distribution of charge in an atom such as the carbonyl oxygen fo the peptide backbone. Perhaps more important however, are the adjustable parameters of the potential energy functions, which assume different values among the various force fields. While the derivation of these parameters from spectroscopic studies, x-ray crystal structures and thermodynamic data suggest adequate comparison with experiment, the ability of various force fields to each predict the stable conformations of a molecule and relative energies thereof is a more relevant comparison.

Barone (2) compared f-y maps for the model peptide N-acetylcyclopropylglycine N-methylamide generated using the ECEPP (3) and AMBER (4-6) force fields. Good agreement was found between potentials with the only notable exception being the absence of a sharp minimum in the helical region of the f-y map generated using the AMBER potential, although this region was still accessible in the AMBER potential.

Roterman (7) compared the energy minimization of the tandemly repeated peptide (Asn-Ala-Asn-Pro)-9 using the CHARMM (8), AMBER and ECEPP potentials. It was concluded that when starting from the same conformation but using any two different potentials, the final conformational resemblance to each other varied from acceptable to highly unsatisfactory and that the ordering of the final conformations and energy differences between them varied greatly for all three potentials. It was found that if the sample of starting conformations was small, predictions based on the three potentials would usually diverge.

While molecular mechanics provides a static snapshot of the molecular system, the dynamic motions of molecules are of fundamental importance in biochemical systems. The dominant role of x-ray crystallography over the past several decades has given rise to the classic view of biomolecules with each atom rigidly defined in space. Of course, the x-ray structure provides only the average atomic coordinates and at room temperature, these atoms exhibit considerable motion about their average position. This might be the rattling motion of an atom group encaged by solvent or in an aqueous environment, a protein may exhibit gross conformational changes lasting several hundred seconds. For proteins, these fluctuations about the average position can play a significant role in protein function. In the oxygen transport protein myoglobin for example, local side chain motions seem to be essential for the entrance and exit of ligands. (9) In the classic case of hemoglobin, allosteric ligand binding induces conformational changes which affect protein activity. (10)

The fundamental principle of molecular dynamics (MD) simulation is the numerical integration of Newton's equations of motion for a system of interactiing particles over a specified period of time. The force acting on each atom is the negative gradient of the potential energy function with respect to atomic positon:

Fi(t) = (-d/dri) V(r1,r2,...,r3Natm)

The acceleration of each atom with mass mi is:

ai(t) = Fi(t)/mi

and the atomic positions are then found by:

d2ri/dt2(t) = Fi(t)/mi

For each atom i there are three scalar equations describing the force each atom j is exerting on atom i:

miai(t) = Sjfij(t)

This equation is solved by replacing it with difference equations which are solved successively for small time steps. For the x coordinate of atom i, the following two Taylor expansions apply:

xi(t + Dt) = xi(t) + vi(t)Dt + ai(t)Dt2/2

xi(t - Dt) = xi(t) - vi(t)Dt + ai(t)Dt2/2

The addition of these two equations gives:

xi(t + Dt) = 2xi(t) - xi(t - Dt) + ai(t)Dt2

and substitution using miai(t) = Sj fij(t) produces:

xi(t + Dt) = 2xi(t) - xi(t - Dt) + Sjfij(t)Dt2/mi

From this equation, known as the Verlet (11) algorithm, the position of atom i at time t + Dt can be predicted if the positions at times t and t - Dt are known.

If the two difference equations are subtracted, the velocity of atom i at time t can be determined:

vi(t) = xi(t + Dt) - xi(t - Dt)/2Dt

The major limitation of the MD technique is that it is approximate. The force acting on each atom is the negative gradient of the potential energy function and this function and its parameters are determined empirically.

In order for these theoretical methods to be useful, they must be evaluated using experimental methods. One such experimental technique is circular dichroism (CD) spectroscopy. CD is an extremely sensitive probe of molecular conformation. This is due to its inherent dependence on the relative spatial orientation of the electronic and mangnetic transition dipole moments of a molecule. As a result, CD spectroscopy can provide useful information regarding the secondary and tertiary structure of proteins in solution.

Model compound studies have been carried out previously in order to correlate chiroptical proerties of peptides in solution with molecular conformation. (12) The primary goal of these investigations was to test the theoretical formalisms used to calculate the chiroptical properites as a function of conformation by comparison with experimental properties. The ultimate goal is to use the techniques developed therein to probe the conformation of biological macromolecules in solution.

A useful experimental system in which to formulate a comparison of force fields is found in the investigation of the CD spectra of cyclic dipeptides. L-Alanyl-L-alanine diketopiperazine (AADKP) is a conformationally restricted, optically active compound with an experimentally definable energy minimum. CD spectroscopy is particularly well-suited to a comparison of this kind given the technique's sensitivity to molecular conformation. While a cyclic dipeptide in an apparently simple system, all the mechanistic possibilities for the generation of optical rotatory power are present (13), even when the problem is limited to the two lowest absorption bands. As the CD spectra of these compounds are readily correlated to molecular conformation, and further, as the spectra in certain cases seem to be the result of a solvent dependent distribution of conformational isomers (12), the ability of a force field to accurately represent this provides one criterion for evaluating the field.

The present research has two complimentary goals. The first is to provide an interpretation of the CD spectra of the model compound L-AADKP on the basis of its conformation and the composition of the solvent system used to obtain the spectra. This will be achieved using the techniques of molecular modeling in conjunction with theoretical CD calculations. To the extent that this is successful, the second goal of the research is to provide for a comparison of the AMBER and CHARMM force fields.

Next, some Background