Simulating Met-Enkephalin

 

Contents

 

A) General

B) Building the structure from the sequence.

C) Preparing files for Amber

D) Preparing files for Amber: Alternate Route

E) Getting going in the Gas Phase:  Minimization

F) Getting going in the Gas Phase: Dynamics 1

G) Getting going in the Gas Phase: Dynamics 2

H) Getting going in the Gas Phase: Simulated Annealing

I) Getting going in the Condensed Phase: Implicit Solvation

J) Implicit Solvation: Trajectory Analysis

K) Getting going in the Condensed Phase: Explicit Solvation

L) Explicit Solvation: Trajectory Analysis

 

A) General

 

This tutorial is set up to run on a unix machine where the programs Amber7, Ribosome, Rasmol, Grace, and VMD are installed and available to all users. If you have difficulty, take a look at the relevant links.  

 

Commands are to the bash and are in italics and are preceded by “>

Commands to a program running in the bash are also in italics and are preceded by “>>

 

All necessary files should be in the tutorial directory that contains this htm file.

 

It is useful to place this directory in your home directory and change into it

 

>cd ~/tutorial

 

Back to Top

 

B) Building the structure from the sequence.

 

We will use the program Ribosome for this purpose

 

Prepare an input file like mete.ribo

 

Now execute ribosome

 

>ribosome mete.ribo mete.pdb $Ribo/src/res.zmat                        # should give you a pdb file.

 

You can take a quick look at it with Rasmol

 

>rasmol mete.pdb

 

Back to Top

 

C) Preparing files for Amber

 

We will use the program Amber7 to do some calculations with our newly prepared structure.

 

>cd parm                                                                                        # we are going to use a subdirectory

>cp ../mete.pdb .                                                                             # get the pdb file we just created

>xleap                                                                                            # starts the amber utilitly xleap

 

>>source leaprc.ff99                                                                     # chooses which force field to use. See manual for more info.

 

>>METE= loadpdb mete.pdb                                                       # creates a “unit” called MET containing the information in the pdb file created in step B.

 

You will notice from the messages that the C-terminal oxygen and all H-atoms that were not present in the pdb file have been added. You will also notice that there was one atom that was unrecognized. We can see what happened by:

 

>>edit METE

 

The molecule control in this mode is not very intuitive. The middle mouse button should allow you to rotate the molecule. The middle and right mouse buttons together should allow you to scale the molecule and the right button alone should allow you to translate it.

 

At the N-terminal end, you should see four protons. This is due to a conflict in naming conventions between Ribosome and Amber. One solution is to choose “Erase” and delete the unbonded proton, which you should now do. Alternatively, we could have named H from the N-terminal tyrosine H1, H2, or H3 in the pdb file before creating the unit (see section D below).

 

You can now choose “Select” and select all or some atoms (with the left mouse button) and then choose the “edit selected atoms” option from the “Edit” menu. Here you can see atom names, types, and charges assigned by amber. Now exit the table and the “Unit editor” to return to Xleap.

 

Some useful commands

 

>>check METE                                                                              # A useful command

>>charge METE                                                                            # another useful command

>>desc METE                                                                                # another useful command

>>desc METE.2                                                                             # another useful command

>>desc METE.2.O                                                                         # another useful command

 

We would now like to create parameter/toplogy and coordinate files for later use.

 

>>saveamberparm METE mete.top mete.crd                               # saves parameter and coordinate files corresponding to the “unit” METE

 

>>saveoff METE mete.lib                                                             # saves the work that we did in xleap

>>quit                                                                                            # prepares coffee :)

 

Back to Top

 

D) Preparing files for Amber: Alternate Route

 

After you gain some experience, you will find it is easier to write scripts to leap. You can do this in either the text version of leap (tleap) or the X-Windows version (xleap). The first option (mete-fix.leap) uses mete-fix.pdb. This is a file with the name of the N-terminal H atom changed to H1.

 

>cd scripts

>tleap

>>source mete-fix.leap

 

The second option (mete-seq.leap) shows you how to build a structure from a sequence within leap as opposed to using ribosome.

 

>xleap

>>source mete-seq.leap

 

You can take at quick look at it

 

>>edit METE

>>quit                                                                                            # time for coffee again

 

Back to Top

 

E) Getting going in the Gas Phase:  Minimization

 

>cd ../../gas

>cp ../parm/mete.top ../parm/mete.crd .                                        # get our parmtop and crd files from step D

 

We will now use the script gas-min to run the input file min.in

 

>./gas-min &

 

Check the output

 

>less min.inf

>less min.out

 

To really see what happened we should look at the structures.

 

We will use the program VMD in order to do this

 

>vmd

 

In the “VMD Main” window choose “File” - “New Molecule”

Click “Browse” and choose mete.top

Set the “File Type” to parm7 and click “Load”

Click “Browse” again and choose mete.crd

Set the “File Type” to rst7 and click “Load”

Click “Browse” again and choose min.rst

Set the “File Type” to rst7 and click “Load”

Close the “Molecule File Browser” and use the “VMD Main” window to look at the structure of Met-Enkephalin before and after minimization.

Quit

 

You should notice that not very much happened to the structure. This is usually the case with minimization. If we wish to see more dramatic changes we should try dynamics.

 

 

Back to Top

 

F) Getting going in the Gas Phase: Dynamics 1

 

We will now use the script gas-md1 to run the input md1.in

 

>./gas-md1 &

 

Check the output

 

>less md1.inf

>less md1.out

 

To really see what happened we should look at the structures again.

 

>vmd

 

In the “VMD Main” window choose “File” - “New Molecule”

Click “Browse” and choose mete.top

Set the “File Type” to parm7 and click “Load”

Click “Browse” again and choose md1.crd

Set the “File Type” to crd and click “Load”

 

You will see the trajectory we just calculated appearing as snapshots every 200fs. Overall you should have 250 structures.

 

Now from the “Mouse” menu choose “Label” - “Bonds” and then click on the N-terminal Nitrogen and then on the C-terminal Carbon.

 

From the “Graphics” menu choose “Labels”. Highlight both visible atoms and click “Hide”. Now change the “Atoms” box to show “Bonds” and select the visible bond. Click “Graph” tab and then the “Graph” button and if you are lucky you will see a graph of the selected bond-length over the trajectory (you need to have Grace working in order to do this).

 

Close the Graph and the “Labels” box. You can now use the movie controls to watch the trajectory with the bond length shown.

 

Stop the movie and go to the “Extensions” menu and click “ramaplot” to open the Ramachandran plot tool. Select the current molecule (mete.top) from the “Molecule” selection box. You will see 5 dots corresponding to the 5 amino acids and the backbone dihedral angles. You can watch the movie again to see how this plot changes over time. Similarly, you can click an individual dot to see a distribution of its values over the simulation.

 

Quit VMD.

 

We will now analyze the output a bit using Grace (called using the command xmgrace).

 

>mkdir md1

>cd md1

>cp ../md1.out .

>~/tutorial/process_mdout.perl md1.out

>xmgrace

 

In xmgrace, choose “Data” - “Import” - “Ascii”. Change the filter to “*.*” and select “summary.TEMP”. This should show you a profile of the temperature over the dynamics run.

 

Now choose “Edit” - “Arrange Graphs” and make 2 rows and 2 columns. Click on the upper right graph and click back on the “Read sets” window (or choose “Data” - “Import” - “Ascii” again).  Selecting “summary.ETOT” will give you a profile of the total energy over the run. To see the Kinetic and Potential components, load “summary.EPTOT” and “summary.EKTOT into the remaining two graphs.

 

When you are satisfied, quit grace

 

Back to Top

 

G) Getting going in the Gas Phase: Dynamics 2

 

In the above procedure we simply set the starting temperature at 300K and let the system go. However, this does not always work. An alternative is gentle heating:

 

We will now use the script gas-md2 to run the input md2.in

 

>cd ~/tutorial/gas

>./gas-md2 &

 

after the job has finished

 

>vmd

 

In the “VMD Main” window choose “File” - “New Molecule”

Click “Browse” and choose mete.top

Set the “File Type” to parm7 and click “Load”

Click “Browse” again and choose md2.crd

Set the “File Type” to crd and click “Load”

 

What differences do you notice between md1.crd and md2.crd ?

Quit VMD

 

>mkdir md2

>cd md2

>cp ../md2.out .

>~/tutorial/process_mdout.perl md2.out

>xmgrace

 

 

Load in summary.TEMP, summary.ETOT, summary.EPTOT, and summary.EKTOT as before.

 

What differences do you notice between md1.out and md2.out ?

 

When you are satisfied, quit grace

 

Back to Top

 

H) Getting going in the Gas Phase: Simulated Annealing

 

As a variation on the heating theme and as a method for finding low energy conformers of the system, we can perform simulated annealing:

 

We will use the script gas-md3 to run the input md3.in

 

>./gas-md3 &

 

after the job has finished

 

>vmd

 

Load in md3.crd in the same way as before and you can see what we did.

 

Quit VMD

 

>mkdir md3

>cd md3

>cp ../md3.out .

>~/tutorial/process_mdout.perl md3.out

>xmgrace

 

Load in summary.TEMP, summary.ETOT, summary.EPTOT, and summary.EKTOT as before to see what the simulated annealing is all about.

 

Back to Top

 

I) Getting going in the Condensed Phase: Implicit Solvation

 

Now we will try to take some account of the solvent. We will do this with one incarnation of the Generalized Born model including a term dependent on the surface area. Although this model treats the solvent as a dielectric continuum, it can be quite useful for many things.

We will begin with simulated annealing:

 

>cd ~/tutorial/gb

 

We will now use the script gb-md1 to run the input md1.in

 

You may have noticed (in the script) that we will begin from our previously minimized structure. We will need this structure, and the topology, to proceed.

 

>cp ../gas/min.rst .

>cp ../gas/mete.top .

 

Run the script

 

>./gb-md1 &

 

While this job runs we can take a look at the next step:

 

Now we want to run some dynamics on our annealed structure. To do this we will reheat it slowly to get stable dynamics

 

We will use gb-md2 to run md2.in

 

If gb-md1 is finished we can run this heating.

 

>./gb-md2 &

 

After gb-md2 is finished we can do some analysis.

 

>mkdir md12

>cd md12

>cp ../md1.out ../md2.out  .

>~/tutorial/process_mdout.perl md1.out md2.out

>xmgrace

 

Load in summary.TEMP, summary.ETOT, summary.EPTOT, and summary.EKTOT to see what the two simulations just did.

 

Quit xmgrace

 

>cd ..

>vmd

 

In the “VMD Main” window choose “File” - “New Molecule”

Click “Browse” and choose mete.top

Set the “File Type” to parm7 and click “Load”

Click “Browse” again and choose md1.crd

Set the “File Type” to crd and click “Load”

Click “Browse” again and choose md2.crd

Set the “File Type” to crd and click “Load”

 

What differences do you see between the gas phase and the Generalized Born model?

 

We are now in a position to run some production dynamics. That is, we have a steady state and we will use it to begin 10ns of dynamics at constant temperature.

 

The script (gb-mdp) and input (mdp.in) for this run are included but as this can take a while to run we will take a look at the already completed output.

You can run the script later to see how long this job takes on your machine.

 

We will first take a quick look at the output:

 

>mkdir mdp

>cd mdp

>cp ../mdp.out  .

>~/tutorial/process_mdout.perl mdp.out

>xmgrace

 

Load in summary.TEMP, summary.ETOT, summary.EPTOT, and summary.EKTOT to see what the production run did.

 

Back to Top

 

J) Implicit Solvation: Trajectory Analysis

 

Now we will analyze the trajectory from the production run.

 

>cd ../traj

 

You should see the trajectory file mdp.crd. As we wrote a snapshot every 2ps for 10ns (mdp.in), we have 5000 frames. This trajectory takes a bit too long to “watch” so we will “downsize” it a bit using the Amber utility “ptraj” and the control script.

 

>ptraj mete.top script 

 

You can now load the resulting “downsized” trajectory (mdp-s.crd) and the associated topology (mete.top) into VMD and watch what happened.

 

Although VMD has some good analysis tools, we will explore the capabilities of ptraj.

 

For example, script1 allows us to monitor the Root Mean Squared Deviation (RMSD), from our starting structure, over the course of the trajectory.

 

>ptraj mete.top script1 

 

will give you 2 data files that you can load into Grace

 

You can see that the deviation has settled down a bit in the second half of the trajectory. We will use this part to calculate an average structure with script2

 

>ptraj mete.top script2

 

If you look at the average structure in either vmd or rasmol, you will see it appears a bit messy. This is normal as, for example, some groups move and rotate more than others and so their average coordinates simply reflect this. As an exercise you can try to minimize this average structure at a later time. For now, we will take a look at the RMSD of our trajectory from the average structure using script3.

 

>ptraj mete.top script3

 

will give you 2 data files that you can load into Grace and compare to the previous 2.

 

Ptraj can also calculate various other properties such as:

 

 

>ptraj mete.top script4

 

 

>ptraj mete.top script5

 

 

>ptraj mete.top script6

 

 

Back to Top

 

K) Getting going in the Condensed Phase: Explicit Solvation

 

In this section, we will see what happens if we include solvent explicitly. In order to do this we will need topology and coordinate files corresponding to a solvated system.

 

We will begin with our annealed GB structure from section I). We will therefore need a pdb file of this structure. There are a few ways to do this but we will use ptraj.

 

>cd ~/tutorial/parm/wat

>cp ~/tutorial/gb/md1.rst .

>cp ~/tutorial/gb/mete.top .

 

As well as being a good analysis tool, ptraj is quite a good converter between various file formats. E.g see ptraj1.in:

 

>ptraj mete.top ptraj1.in

 

gives us a pdb file

 

We will now build our system in leap using the script mete-wat.leap

 

>cp md1.pdb.1 ../scripts

>cd ../scripts

>xleap

>>source mete-wat.leap

>>edit METE

 

You will notice that this gives a box of water. We will apply periodic boundary conditions to this box to simulate a dilute solution.

 

>>quit

>cp mete-wat.crd ~/tutorial/wat

>cp mete-wat.top ~/tutorial/wat

>cd ~/tutorial/wat

 

We will begin by minimizing the water molecules with positional restraints on the solute. This is mainly to fix any bad contacts or holes formed during the water placement:

 

We will do this using wat-min1 to run min1.in

 

>./wat-min1 &

 

Depending on the speed and load of your machine, this may take a little while.

 

When the minimization is finished we will then run some gentle heating dynamics with constant volume (NVT), also using restraints on the solute in order for the water to adjust properly to the presence of the solute and the box boundaries.

 

We will do this with wat-md1 and md1.in

 

>./wat-md1 &

 

You will notice that we only ran this for 10ps.Ordinarily you would probably do 50ps as we did for the gas phase example.

 

In order for the restrained solute to “feel” the equilibrated water gradually, it is common to perform successive minimizations with reduced positional restraints.

 

Once wat-md1 is finished we will use wat-min2 to run min2.in and min3.in

 

>./wat-min2 &

 

Again, one would ordinarily do more steps in each minimization (e.g. Until DRMS<0.1) and more minimizations.

 

Now we will heat the whole system (without restraints) to 300K with constant volume (NVT)

 

Once wat-min2 is finished we will use wat-md2 to run md2.in

 

>./wat-md2 &

 

We will now read in the coordinates AND velocities of the constant-temperature, constant-volume (NVT) run in order to obtain a realistic density. That is, by running with constant temperature and pressure (NPT), the volume of the simulation box is able to adjust and correct the density.

 

Once wat-md2 is finished we will use wat-md3 to run md3.in

 

>./wat-md3 &

 

While wat-md3 is running, we can have a quick look at the results of wat-md2 (The NVT run).

 

>mkdir md2

>cd md2

>cp ../md2.out  .

>~/tutorial/process_mdout.perl md2.out

>xmgrace

 

Load in summary.TEMP, summary.ETOT, summary.EPTOT, and summary.EKTOT to see the unrestrained heating phase with explicit solvent.

 

Once wat-md3 is finished we can take a look at the results of the NPT run.

 

>cd ..

>mkdir md3

>cd md3

>cp ../md3.out  .

>~/tutorial/process_mdout.perl md3.out

>xmgrace

 

On this occasion you can make 3 rows and 2 columns and read in summary.TEMP, summary.PRES, summary.VOLUME, summary.DENSITY, summary.EPTOT, and summary.ETOT.

 

You should notice a drop in Volume as the pressure coupling starts to work. As the mass of the system is constant, it is hardly surprising to see the associated increase in density. This is a common feature of equilibration with explicit solvent. That is, the solvent placement algorithm almost always results in too low a density. The only way to reliably fix this is to run with NPT until the density stabilizes.  

 

We are now in a position to run some production dynamics. That is, we have a steady state and we will use it to begin 10ns of dynamics at constant temperature and pressure.

 

The script (gb-mdp) and input (mdp.in) for this run are included but as this can take a while to run we will take a look at the already completed output.

You can run the script later to see how long this job takes on your machine.

 

We will first take a quick look at the output:

 

>mkdir mdp

>cd mdp

>cp ../mdp.out  .

>~/tutorial/process_mdout.perl mdp.out

>xmgrace

 

Load in summary summary.TEMP, summary.PRES, summary.VOLUME, summary.DENSITY, summary.EPTOT, and summary.ETOT to see what the NPT production run did.

 

Back to Top

 

L) Explicit Solvation: Trajectory Analysis

 

 Now we will analyze the trajectory from the production run.

 

>cd ../traj

 

You should see the trajectory file mdp.crd. This file has already been “downsized” and contains 500 snapshots covering 10ns.  

 

You can load this trajectory (mdp.crd) and the associated topology (mete-wat.top) into VMD. However, as this was a periodic NPT run, the trajectory file contains the box size associated with every snapshot. Therefore, you need to set the “File Type” of mdp.crd to crdbox to see something sensible.

 

You will notice that the system appears to be spread over quite a lot of space. This is perfectly normal as over time the water molecules tend to drift away from the central box. This is not a problem as long as all the coordinates “image” back into the central box.

 

We can use ptraj and scr1 to accomplish this

 

>ptraj mete-wat.top scr1 

 

You can load the imaged trajectory (mdp-image.crd) into VMD remembering to use the crdbox “File Type”.

 

Although this is useful to understand something about periodic boundary conditions, in order to analyze our polypeptide, it is helpful to remove the water from the trajectory with scr2.

 

>ptraj mete-wat.top scr2 

 

You can load the stripped trajectory (mdp-strip.crd) into VMD if you wish but as we no longer have water molecules we need to use mete.top instead of mete-wat.top. We also no longer have box information so the “File Type” for the trajectory should be just crd and not crdbox.

 

We can now perform the same kind of analysis that we did for the implicit solvent in section J). As with vmd, ptraj needs a topology that corresponds to the trajectory so for analysis of mdp-strip.crd we need to use mete.top.

 

We can now take a look at the RMSD of our trajectory from the starting structure using scr3

 

>ptraj mete.top scr3 

 

will give you 2 data files as before.

 

We can also combine the calculation of the average structure from the second half of the trajectory, the positional fluctuations, the geometric properties, and the hydrogen bond characteristics into one control script (scr4). As an exercise, you can try to minimize your average structure.

 

>ptraj mete.top scr4 

 

and get the RMSD of our trajectory from the average structure with scr5

 

>ptraj mete.top scr5

 

Try to compare all of these properties from the explicit solvent run with those of the implicit solvent run we calculated in section J).  In particular, how do the average structures on the second halves of the respective trajectories compare to each other?

 

Finally, we can also analyze some aspects of the explicit solvation of particular groups in our peptide.

 

scr6 calculates the radial distribution function (rdf) of water oxygens around the hydroxyl proton of the Tyrosine residue (output as hh-rdf_standard_xmgr). In addition, by declaring the bond involving this proton as a potential “H-bond acceptor”, we can analyze the H-bonding tendencies between this bond and the water molecules.

 

>ptraj mete-wat.top scr6

 

It is interesting to compare this to the same analysis for the para-proton on the Phenylalanine residue. scr7 calculates the rdf of the water oxygens around this proton (hz-rdf_standard_xmgr ) as well as counting possible H-bond interactions.

 

>ptraj mete-wat.top scr7

 

Try to load hh-rdf_standard_xmgr and hz-rdf_standard_xmgr into the same graph using xmgrace. Although these curves are a bit rough (due to too few snapshots) you should see a very clear difference between them.

 

 

Back to Top