This exercise is intended to provide some experiences with the following operations:

(a) Construction of a small but biologically relevant molecule,

(b) Study of the conformational dynamics of that molecule in vacuo,

(c) Construction of a solvent shell around the molecule,

(d) Analysis of the same conformational dynamics within the solvent shell.

A small molecule is being used for the exercise so that the computation time needed is reasonably short. The same approaches can be used with large molecules, but the dynamics runs could become multi-day affairs on the SGI workstations in the lab.

The molecule to be studied is adenosine 3',5' cyclic monophosphate (cAMP). This species acts as a "second messenger" within cells; its formation is catalyzed by adenylate cyclase and is highly regulated by hormones and other signaling molecules. The conformational processes to be explored in this exercise are the "flutter" of the ribose ring pucker (the interconversion of exo-endo conformations of the ring) and rotation about the C-N single bond (the glycosidic bond) that links the base to the ribose ring. The endo conformation has the out-of-plane ring atom on the same side of the ring as C_5' while the exo conformation features the out-of-plane atom on the side opposite C_5'. Rotation of the base relative to the ribose ring interconverts the syn and anti conformations. By convention, the position of the base relative to the ribose ring is specified by a conformational angle c which is defined at the O4' - C1' - N1 - C2 torsion angle.

Some suggestions for accomplishing the tasks listed above using SYBYL are given below. If you use different software for the exercise, your report should indicate this.

1. Get the coordinates for cAMP

This structure can be built de novo, constructed by modifying an adenosine residue built with the Biopolymer >>> Build feature, or extracted from a CSD or PDB structure file. Be sure to load the charges for the structure (COMPUTE >> CHARGES, MMFF94 charges). Minimize your cAMP structure using the MMFF94 force field until the forces on all atoms are less than 1 kcal mol^-1 A^-1. Expand the drawing of the molecule so that it fills about 2/3 of the screen. Save a file containing the structure under an appropriate name, using the .mol2 format.

2. Run a dynamics simulation for the in vacuo system

This simulation will in essence be for the molecule in the gas phase since no other molecules will be present to interact with the cAMP. Select Compute >>> Dynamics >>> Setup Dynamics. Set the parameters for a NTP (constant temperature and pressure) simulation at a pressure of 1 atm and a temperature of 300^o, using 10000 1 fs integration steps. Take a snapshot every 5 fs. (Every 5 fs the Cartesian coordinates of the molecule at that time will be added to a data file.) Set the coupling times to 50 fs. Set the force field to MMFF94. Enter a descriptive name for the job name and click OK. Observe the screen image, which is updated periodically as the dynamics simulation proceeds. (Quiz Question: Why does the CH₂ group of the molecule move so much more vigorously that the PO₂ group?) It will take about a couple of minutes for the simulation to run.

You can play back the dynamics trajectory by Analyze >>> Dynamics >>> Replay Trajectory. The atom colorings for this replay can be rather psychedelic!

3. Analyze conformational dynamics

SYBYL uses a database strategy to provide a format for analysis of dynamics simulations ("trajectories"). To get out the conformational angles of interest as a function of time click Analyze >>> Dynamics >>> Analyze Dynamics… and select the file that was saved during the dynamics run (specified earlier). Select a name (table name) for the database when asked and choose Data Reduction interval equal to 1. Accept 0 as the data starting value and accept 100000 as the data ending value. (Click on OK to select the standard values.) A window with a database will come up. It will show time, temperature and energies as a function of the simulation step. Move the database window and the structure drawing so that both are observable at the same time. Click on column 6 of the data base--it will turn black when selected. Then click AutoFill and select Torsion from the choices offered. Enter a name for column 6 that indicates it will contain the c (chi) torsion angle. When asked to select atoms, click on C2', C1', N1 and C2 of the cAMP drawing. The program will go away for awhile as it goes through the saved coordinate sets and computes this angle for each set. When completed, you can scroll down the database to see how this angle changes with time. Now highlight column 7 and choose to define another torsion angle. This time select C4', C3', C2', and C1' as the angle to be calculated.

The data of the database could now be exported to a file which could be read by SigmaPlot or some other graphing program. To prepare a plot using SYBYL, click Graph then Scatter. Select for the X-Axis time, for the Y-Axis Chi, for the Z-axis Omit Axis, for Color Uniform. Click Layout…, observe the parameters definable there and click on Connecting Lines. Select a display area and click Create. SYBYL creates a plot of chi vs. time in the selected drawing area. This can be rotated or translated like a molecule drawing and printed as any SYBYL graphical object. This plot is regarded as a "background" and can be erased by Delete All Backgrounds.

Look at plots of temperature vs. time and the ribose ring pucker angle vs. time. The starting structure was not thermally equilibrated, with the result that energy and structure change during the first part of the simulation until the kinetic energy and potential energy achieve the appropriate values for the specified simulation temperature. Since the temperature was changing during the course of the simulation any conformational conclusions drawn from this first dynamics trajectory are therefore suspect. More reliable conclusions can be reached by running a a second simulation for 10000 1 fs time steps using the results of this first one as the starting point. Save the results of this second calculation in a separate file from the first one for analysis. Analyze it using the SYBYL database functions described above.

Make hardcopy plots of the conformational angles vs. simulation time. Analyze changes of the torsion angles to estimate the rates at which the syn-anti and exo-endo interconversions take place. (Recall that a rate is number of events per unit time.) Turn in these estimated rates along with the hardcopies. Note: It is likely that each person will get slightly different answers and plots, so comparisons to what someone else may have obtained should be done gingerly.

4. Simulation of cAMP in water

Clear all molecules (ZAP). Load the data for the cAMP structure (from the previously written out file) into molecule area M1. Then create a water molecule in molecule area M2. (There is a constructed water molecule at the bottom of the fragment list that comes up when Get Fragment… is selected under Build/Edit.) Assign MMFF94 charges to both molecules. There are many options for construction of solvent shells in SYBYL. For this exercise, a process called Molecular Silverware will be used to layer solvent molecules around the solute. Under Build/Edit select Solvent, then Solvate… Select M1 (the cAMP molecule) to solvate and M2 as the solvent. Choose Droplet, and accept Tripos as the radius. Choose to create 1 solvent layer and observe the result. This process can be repeated until the desired thickness of solvent molecules surrounds the solute. For this exercise, create a solvent "droplet" of about 10 A radius with the cAMP molecule at the center of the droplet.

Even though the algorithm used to position the solvent molecules is fairly sophisticated there will be bad constants between some solvent molecules and between solvent molecules and the solute. Minimize the energy of the collection of solute + solvent molecules created to remove these contacts before starting a dynamics calculation.

Run a dynamics run in the same way that was used for cAMP in vacuo. Analyze conformational motions in solvated cAMP as was done previously. Report the same results as was done for the in vacuo simulation and comment on the differences in the results from the two types of calculations.

5. General Comments

[1] You may wish to work through the SYBYL tutorial on introductory dynamics. This can be brought up in the web browser by going to the TriposBookshelf >>> Force Field >>> Force Field Tutorials >>> Introduction to Dynamics Tutorial. A useful tutorial on use of the "Molecular Spreadsheet" (data base) for analysis of dynamics simulations can be found at TriposBookshelf >>> Force Field >>> Force Field Tutorials >>> Advanced Dynamics Tutorial, section 2.5.4.

[2] There are several ways of creating solvated ensembles in both SYBYL. The solvent molecules need not be water molecules.

[3] The number of water molecules in the simulations done for this exercise is rather small. Typically, one wants to have a "layer" of water molecules that extend significantly beyond the edges of the solute molecule. For proteins, the number of solvent molecules can therefore become rather large--in the thousands to tens of thousands, with the calculations running correspondingly slower.

[4] The cAMP-solvent simulations are not representative of what would normally be done when doing calculations that include solvent. The number of solvent molecules is small and no attention has been given to the errors introduced by the truncation of force field terms at the edges of the model. These latter effects are typically handled by the periodic boundary conditions (PBC) that will be discussed in class.

[5] Notice the behavior of water molecules in the cAMP-solvent simulations. If this calculation were being done more rigorously, the "evaporative" tendencies of these molecules would have to be reduced. To prevent these excursions, a weak restoring potential can be applied to the oxygen atoms of the solvent molecules near the edge of the model. The potential has to be chosen carefully since it can cause the model to collapse in on the solute molecule if it is chosen to be too strong. A second approach to keeping solvent molecules in place is to create a model in which there is a layer of mobile solvent around the solute which is then enclosed by another layer of solvent molecules that is kept rigid (not included in the dynamics). Collisions of mobile solvent molecules with the rigid solvent molecules keep things in place without using artificial constraints or potentials.

[6] When estimating rates, each dynamics simulation produces a single observation of the behavior of an angle or interactomic distance. To get good statistics, the experiment/simulation would have to be repeated several times and the results of each run compared. And to get reliable rate, the sensitivity of the calculated rates to the assumptions of the model would have to be investigated.

[7] Check on what SYBYL is adding to the contents of your directories. Dynamics simulations run with SYBYL create a number of data files, some of which can become very large. Remember that there are time and size limitations on the files that can be stored on the lab computer system; good housekeeping is essential. Translation: Remove files that are not needed as soon as possible.

6. A report of your results is due February 20.