BioSci 1540

Introduction to Biomolecular Simulation

Model Building

Some initial notes

  • Throughout the tutorial, terminal refers to the application in /Applications/Utilities.
  • The Tk Console is a terminal-like window that you can access in VMD via the drop down menu Extensions in the VMD Main window.
  • The VMD atomselect command requires that you either specify a molid or top (e.g. atomselect top "SELECTION" or atomselect 3 "SELECTION"). If you are using top as is used below, make sure that the T appears next to the correct model.
  • Anytime you use the atomselect command in the Tk Console, it is recommended to create a new representation using the Create Rep button in Representations window of VMD, and enter the SELECTION in the Selected Atoms field to visualize the subset of atoms to make sure it is selecting what you think it should be.

Download files

Begin by downloading the archive of files from lab github repository, which can be found at: https://github.com/synapticarbors/BioSci1540/zipball/master

Move the .zip archive from where it was downloaded (most likely the Downloads directory) to a new folder that you will use for the remainder of the labs. Preferably, this will be in a location that is backed-up to the class server so that you can access it on other machines in the lab. Once you have moved the archive, which will be named something like synapticarbors-BioSci1540-XXXXX.zip, you can extract the files by either double-clicking on it and then rename the folder md_tutorial, or opening a terminal window and running the following commands:

Terminal
1
2
3
cd <Location of .zip file>
unzip synapticarbors-BioSci1540-XXXXX.zip
mv synapticarbors-BioSci1540-XXXXX md_tutorial

Be sure to replace the XXXXX in synapticarbors-BioSci1540-XXXXX by whatever the names of the files are for your specific case. For the rest of the lab, the path to the md_tutorial directory will be written in shorthand as $MDT.

Construct initial model for gramicidin in membrane and solvent

If we were setting up a system entirely from scratch, we would likely begin with a set of coordinates from the Protein Data Bank. The experimental structure might be missing the coordinates for a subset of the residues, and would be lacking solvent and lipids and probably would not contain the positions of any hydrogens in the protein. Due to the time constraints of the lab, you will begin with a system in which gramicidin has already been embedded in a DMPC lipid bilayer with water and ions, although hydrogen atoms have not been added yet. Take a moment to look at the page for 1JNO structure of Gramicidin A on the PDB website, whose coordinates provide the initial structure of the protein. http://www.pdb.org/pdb/explore/explore.do?structureId=1JNO

QUESTION 1: What method was used to determine the structure of 1JNO?

Create topology file and add hydrogens

While the .pdb file tells NAMD the initial coordinates to begin the simulation with, we also need to provide it with a topology file that tells it which atoms are bonded, which are connected via an angle spring, etc. The topology file carries a .psf extension and can be created using the psfgen tool that comes with VMD. This tool can also be used to add hydrogen atoms to the system.

In the Tk Console, first change directories to $MDT and then load the starting coordinates:

Tk Console
1
2
cd $MDT/model_files
mol new gramicidin_initial_struct.pdb

Take a moment to look at the structure, change the representation to highlight the protein, lipid, water, and ions.

The first thing we have to is split the structure into segments and create a new .pdb file for each. We will use the VMD command atomselect to select a subset of atoms and then write that subset to a .pdb file

Tk Console
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
set wat [atomselect top "water"]
$wat writepdb wat.pdb

set lip [atomselect top "lipid"]
$lip writepdb lipids.pdb

set ga1 [atomselect top "segname GA1"]
$ga1 writepdb  ga1.pdb

set ga2 [atomselect top "segname GA2"]
$ga2 writepdb  ga2.pdb

set cl [atomselect top "resname CLA"]
$cl writepdb  cl.pdb

set pot [atomselect top "resname POT"]
$pot writepdb  pot.pdb

Now load psfgen

Tk Console
1
2
package require psfgen
resetpsf

Load the topology file that shows how the atoms in different residues and groups should be built, and define some aliases:

Tk Console
1
2
3
4
5
topology ../toppar/top_all27_prot_lipid_gram.inp

pdbalias residue DLE LEU
pdbalias residue DVA VAL
pdbalias atom ETA O OG

Define the segments:

Tk Console
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
segment GA1 {
first NONE
last NONE
pdb ga1.pdb
}

segment GA2 {
first NONE
last NONE
pdb ga2.pdb
}

segment TIP3 {
first NONE
last NONE
pdb wat.pdb
}

segment MEMB {
first NONE
last NONE
pdb lipids.pdb
}

segment CLA {
first NONE
last NONE
pdb cl.pdb
}

segment POT {
first NONE
last NONE
pdb pot.pdb
}

Load in the coordinates from the .pdb files we previously created:

Tk Console
1
2
3
4
5
6
coordpdb ga1.pdb GA1
coordpdb ga2.pdb GA2
coordpdb wat.pdb TIP3
coordpdb lipids.pdb MEMB
coordpdb cl.pdb CLA
coordpdb pot.pdb POT

Guess the coordinates of any atoms specified by the topology file, but missing in the coordinate files (e.g. hydrogen atoms):

Tk Console
1
guesscoord

and finally write out the new coordinates and topology file:

Tk Console
1
2
writepsf gramicidin_model.psf
writepdb gramicidin_model.pdb

Now load the new model into VMD to visualize it and make sure everything looks good:

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

QUESTION 2: Did psfgen do a good job of selecting the positions of the hydrogen atoms for the water? Why or why not? (HINT: It may be easier to visualize the water if it is shown in the VDW representation).

Place a potassium atom in the channel

In order to investigate the behavior of a potassium ion (K+) in the channel, each group in the lab will place an ion at a different position in the channel. We will do so by swapping one of the K+ atoms in the bulk solvent for a water molecule in the channel.

First determine the extent of the protein by selecting the backbone atoms and then measure their minimum and maximum values along the z-direction:

Tk Console
1
2
set prot [atomselect top "protein and backbone"]
measure minmax $prot

This will produce the minimimum and maximum values of the positions in x,y,z as two lists {xmin,ymin,zmin} and {xmax,ymax,zmax}.

Use the zmin and zmax values to create a new representation in the Graphical Representations window showing just the water molecules in the channel. First try the selection water and z > zmin and z < zmax, where zmin and zmax have been replaced with the actual values from the last step. You should see that there are several water molecules that lay outside of the channel. To focus on the water in the channel, take a look at the description of the same and within keywords in the VMD manual, and modify selection to show only the water molecules close to the long axis of the channel.

QUESTION 3: Record your final atomselection.

There should be nine water molecules that are in the channel in a single file line along the z-axis. Each group will be assigned a different water molecule to swap with a K+ ion, so that we can analyze the stability of K+ in the channel at different positions. After you have been assigned a water molecule, determine its resid by pressing 1 and then use the mouse to select and label the oxygen atom of the water molecule. Open up the Graphics > Labels window to see more information about the selected atom and use the ResID to then make selection in the Tk Console:

Tk Console
1
2
set watsel [atomselect top "resame TIP3 and resid XXX"]
set potsel [atomselect top "resname POT and resid 1"]

Now determine the center-of-mass position of the water and K+ ion:

Tk Console
1
2
set watpos [measure center $watsel]
set potpos [measure center $potsel]

Then determine the vector between them:

Tk Console
1
set r [vecsub $watpos $potpos]

and swap their positions:

Tk Console
1
2
$potsel moveby $r
$watsel moveby [vecscale $r -1]

If everything looks correct, then save a new pdb file that contains the coordinates with the K+ in the channel:

Tk Console
1
2
set all [atomselect top "all"]
$all writepdb gramicidin_final_model.pdb

Since the topology has not changed, we can continue to use the .psf file we created in the previous step.

Create restraint files

The final step is to create a file that restrains the positions of the protein and K+ ion in the channel during the minization and heating of the system, so that these steps do not artifically perturb them. This file is just a .pdb file in which we use the column that specifies the beta factor to indicate which atoms are restrained. First delete all of the molecules that have been previously loaded into VMD and then:

Tk Console
1
2
3
4
5
6
7
8
9
10
mol new gramicidin_final_model.pdb

set all [atomselect top "all"]
$all set beta 0.0

set prot [atomselect top "not {water or ion or lipid} and noh"]
$prot set beta 1.0

set pot [atomselect top "resname POT and resid 1"]
$pot set beta 1.0

Now take a look at your system using a representation where the Coloring Method has been set to Beta. You should see the atoms that you have selected to restrain in one color and all of the other atoms in another. If this looks correct save the file:

Tk Console
1
$all writepdb gramicidin_res.pdb

We will also create a restraint file to use during our production run, which will apply a weak restraint to the CA atoms in gramicidin’s backbone. This is meant to keep the protein from diffusing in the plane of the membrane. This will make our analysis easier since we will not have to correct for the periodic boundary condition if one of the monomers cross the boundary when the other has not. The weak restraint should have minimimal effect on our results, although this strategy is not commonly used in research (we are only doing this for the purposes of this lab).

Again using the structure from gramicidin_final_model.pdb, we will reset the beta values and select a different set of atoms:

Tk Console
1
2
3
4
5
6
7
set all [atomselect top "all"]
$all set beta 0.0

set ca [atomselect top "protein and name CA"]
$ca set beta 1.0

$all writepdb gramicidin_res_prod.pdb

We are now ready to move onto Minimization and Equilibration.