BioSci 1540

Introduction to Biomolecular Simulation

Analysis

You should now have a long trajectory in your $MDT/production directory called dyn00.dcd. In this session we will use VMD to analyze this simulation.

Visualization of production run

Start off by loading your original .psf and .pdb. The command below assumes that you are in your $MDT/model_files directory:

Tk Console
1
2
mol new gramicidin_model.psf
mol addfile gramicidin_final_model.pdb

Now change directories so that you are in $MDT/production and load the trajectory:

Tk Console
1
mol addfile dyn00.dcd waitfor all

Take some time and using different representations and selections, play the trajectory (the controls to play the loaded trajectory are on the VMD Main window) and see what is going on with the water, the ions and protein. You can select the ions by using the ion keyword in the selection or using their resname (resname POT for potassium or resname CLA for chloride). It is often easiest to visualize the ions if you set the Drawing Method to VDW in the Graphical Representations > Draw Style tab.

QUESTION 1: Write down a qualitative description of what you see, and if there are any interesting features of the simulation that you notice (max. 1 paragraph)

Water String Dipole

If you closely examine the string of water in the pore, you will notice that the water molecules are not randomly oriented. Instead the oxygen of one water points toward the hydrogens of the next. This directs the negatively charged oxygen in one direction with the positively charged hydrogens in the other.

QUESTION 2: Predict what happens to the orientation of the water as positively charged ion like potassium moves through the pore. When the ion is at the center of the pore, will the water above and below it be oriented in different directions?

We can quantify this orientation by measuring the net dipole of the string of waters. The following tcl script will calculate the z-component of the net dipole for each snapshot in your trajectory and write the results to a file:

Tk Console
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
set pwat [atomselect top "water and same resid as {water and within 4 of protein and z > -8 and z < 8}"]
set nframes [molinfo top get numframes]

set fout [open dipole.dat w]

for {set k 0} {$k < $nframes} {incr k} {
$pwat frame $k
$pwat update

set dipole [measure dipole $pwat -debye]
set dz [lindex $dipole 2]
puts $fout $dz
}

close $fout

Plot the data by launching the application Grace which is located in the Applications folder.

  • Open dipole.dat by selecting Data > Import > ASCII.. and then selecting dipole.dat as the file. You may have to navigate to $MDT/analysis to locate it.
  • Select the options:
    • Set type: XY
    • Load as: Single Set
    • Click OK to plot the data

QUESTION 3: Does the net dipole change during the simulation? If so, does it change smoothly in time or does it make a sudden jump? Give a qualitative description of what you think is going on in this plot.

Save an .jpeg image of the plot and put it in your analysis directory. To save the image, in Grace, under the File > Print setup.. dropdown, select Device: JPEG, and set the path in the File name element. Hit Apply then Accept and then Ctrl+P. You should now see the .jpeg file in the directory you specified.

Ion Position

The previous section shows how you can write tcl code to loop over all of the frames in the trajectory to measure the dipole moment of the water string as a function of time.

The basic steps involve:

  • Opening a file to write data to: set fout [open <filename> w]
  • Making an atom selection involving the atoms you want to track: set sel [atomselect top "SELECTION"]
  • Getting the number of frames in the trajectory: set nframes [molinfo top get numframes]
  • Writing a loop to iterate over the frames, inside of which you:
    • Update the frame: $sel frame $k
    • Extract some data
    • Write the data to file: puts $fout $data
  • After finishing the loop, closing the file: close $fout

Write some tcl code using the above ingredients to get the z-position of the ion that was originally placed in the middle of the channel. You will have to define the atom selection (we did this in a previous section when we were moving the potassium into the channel), and then get its z-position, which can be done using something like: set zpos [$sel get z]

Write down the code you wrote in the same file where you have been recording the answers to earlier questions, and also plot the data using Grace as we did in the previous section.

Addendum

If you are interested in a more detailed and careful analysis of this system, I recommend reading Energetics of ion conduction through the gramicidin channel by Allen, Andersen and Roux.