Intro to Molecular Modeling-Sketching, MM, MD and Docking

GETTING STARTED:

1) The instructor will provide you with student accounts and passwords. These accounts will remain active for two weeks in case you would like to try more things on your own.

2) Log onto the SGI workstations with your account and password

3) Materials for the training exercises will be available in the directory "MGMcourse".  To access this directory, type the following command at the command prompt:

· cd MGMcourse

and hit [enter]

4) For these exercises, you will be using the SYBYL suite of software.  To run SYBYL, type:

· sybyl

and hit [enter]
 

EXERCISE 1:  MOLECULE BUILDING USING SYBYL

In this exercise you will be asked to draw and save the following molecule:

 

molecule

 

1) To begin, locate the {Build/Edit} button on the top bar of the SYBYL interface and click it.

 

2) On the pull-down menu, you will see a button for {Sketch Molecule...} -- click that too.

Note:  in the future, I will use the following shorthand notation:

* {Build/Edit} >> {Sketch Molecule...}

to simplify the description of steps like 1) and 2)

 

3) SYBYL can handle as many as 10 different molecules in memory at any given time.  It does this by storing the different molecules in bins called "Molecular areas" numbered M1 through M10.  SYBYL is asking you which Molecular area you would like to use for your first molecule.  M1 is highlighted by default.  Simply click {OK} to select M1.

 

4) You now have the three widgets used for building molecules.  On the middle widget, click the {Draw} button.  You are now ready to begin drawing.

 

Here are some hints for drawing:

* In "Draw" mode, every time you click on the black screen, you will draw an atom

* If the atom is enclosed in a small yellow square, it means that the next atom you draw will form a bond to the atom with the yellow square

* If you do not want your next atom to be bound to the atom with the yellow square, click inside the yellow square and it will go away

* If you want your next atom to be bound to an atom that does not have a yellow square around it, click on the atom you want to bind to and that will give it a yellow square

* If you want to form a bond to two already existing atoms, simply click the two atoms in succession

* If you want to form a multiple bond between two atoms, click on the two atoms in succession again

 

5) Begin inputting atoms where you would like to see them in your molecule.

Note:

* When you first input your atoms, they will all be Carbon by default.  Leave them that way for now -- it will be easier to change them from Carbon to other elements later than it will be while you're sketching.

* You should not bother drawing in any atoms that will correspond to Hydrogens. Hydrogens can be added automatically later.

 

6) Sketching molecules can be tricky and many people end up with a few bonds that they did not intend to sketch.  If you would like to delete such spurious bonds:

· Go to the middle drawing widget and click {REMOVE_BOND}

· Click on the two atoms whose bond you wanted to delete

Note:

* As soon as you click {REMOVE_BOND}, you are no longer in "Draw" mode -- you are in "Delete Bond" mode and will keep on deleting bonds with every click until you get back to "Draw" mode.

* To get back to "Draw" mode, go to the middle drawing widget and click {Draw}

 

7) If you have any extra atoms that you don't want:

· Go to the middle drawing widget and click {REMOVE_ATOM}

· Click on the atoms that you would like to delete

· Don't forget to go back to "Draw" mode when you are done deleting!  Otherwise you will continue to delete atoms

 

8) Once you have placed atoms corresponding to all of the heavy atoms in flavopiridol, and have the right bonds between all of them, you are  ready to give them the right atomic type.

· Flavopiridol has 3 Oxygen atoms.  In order to change atoms from Carbon to Oxygen, click on the {O} in the bottom widget and then click on each of the atoms you would like to change.

· Do the same for the Chlorine and Nitrogen atoms.

 

9) You are now ready to add H's.  To do this, simply go to the middle widget and click on {ADDH}.

Note:

* if you have made any mistakes in specifying your bonds (e.g., drawing a single bond where you wanted a double), you will get the wrong number of Hydrogens.  Check your structure carefully.

 

10) There are numerous other molecular drawing options in the sketch tool.  These are described in detail in the online documentation that comes with our Tripos license.  To access this online documenation, minimize your SYBYL window, click on the Netscape icon, and find the Bookmark entitled {Sketcher Menu Items}. 

 

11) For now, however, let's exit the molecular sketcher by clicking the {Exit} button on the middle widget.

 

12) You now have your molecule alone on the black screen.  You can look it over carefully by manipulating it in the screen using your mouse.      Specifically:

· If you hold down the right mouse button and move your mouse, you will rotate the molecule

· If you hold down the middle mouse button and move your mouse, you will translate the molecule

· If you hold down both the middle and right buttons and move the mouse up or down, you will focus in, or pan away from the molecule

 

13) If you notice any serious chemical errors (incorrect atoms or bonds) in your molecule, go back to the sketcher as follows:

· {Build/Edit}  >>  {Sketch Molecule}  >>  select M1  >>  {OK}    and make the relevant corrections.

 

14) In step 12 you probably noticed that the molecule does not have the proper three dimensional structure.  To correct this, you are going to need to do a geometry optimization as part of the next tutorial

 

EXERCISE 2:  GEOMETRY OPTIMIZATIONS USING CLASSICAL MECHANICS

In this exercise, you will optimize the structures of two different molecules -- flavopiridol and taxol:

1) Starting the with flavopiridol molecule you sketched in Exercise 1,

   bring up the minimization widget as follows:

· {Compute}  >>  {Minimize...}

 

2) The default settings for SYBYL's classical mechanics optimizer are very good for straightening out the three dimensional structures generated by the molecule sketcher.  So, without changing any settings, go to the "minimize" widget and simply click {OK}.

 

3) The molecule will optimize fairly quickly.  Once it is done, check it over to convince yourself that the three dimensional structure looks more sensible now.

 

4) In case you need the flavopiridol molecule later, you should save it as follows:

· {File}  >>  {Save As}

· In the {File:} box, enter the name "flav.mol2"

· Select m1  >>  {Save}

 

5) Now that you have saved your flavopiridol work, you may delete the molecule from SYBYL as follows:

· {Build/Edit}  >>  {Zap (Delete) Molecule}

 

6) You are now going to optimize a more complicated molecule -- taxol.  In this case, you will not have to sketch taxol by hand because there is already a crystal structure available for the molecule.  You can load this structure as follows:

· {File}  >>  {Read...}

· Click on "1JFF.pdb" from the list of files

· Select m1 >> {OK}

· Click {NO}  (do NOT center the molecule}

 

7) What this gives you is the entire crystal structure for taxol docked in the protein tubulin.  For now you will not need tubulin, so you will need to extract taxol from the crystal structure.  Bring up the "Extract" widget as follows:

· {Build/Edit}  >>  {Extract...}

 

8) On the widget, click {Substructures}

 

9) This gives you a "Substructures" widget.  Scroll down this widget all the way to the bottom and you see the Taxol moiety.  Click on {TA1601} and click {OK}.

 

10) From the main "Atom Expression" widget do the following:

· {OK}

· Select M2 >> {OK}

 

11) At this point, you will have extracted taxol from tubulin and stored it in a new molecular area.  To get rid of the unwanted tubulin:

· {Build/Edit}  >>  {Zap (Delete) Molecule}  >>  select M1 >> {OK}

 

12) If the whole screen appears to have gone blank, don't worry!  The taxol is probably still in there, but is just outside of the viewing area.  To find it, use your mouse to pan out.  I.e., hold down both the middle and right mouse buttons and pull the mouse towards yourself.

 

13) Once you've found the taxol, try to center it with your mouse.  Once it's centered, try rotating it around.  If it keeps sliding way off to the side when you rotate, ask the instructor to "show you a trick".

 

14) The molecule may be colored monochromatically at this stage.  In order to be able to see the different atom types:

· {View}  >>  {Color}  >>  {By Atom Type}

 

15) When you bring up a molecular structure straight out of a PDB file as you have in this case, SYBYL doesn't always get all of the bonds described correctly.  Check the structure over to see if you see anything out of the ordinary.  One big problem you should see is the bonding in a six-membered ring that (if you know the taxol structure) you know should be a phenyl ring.  To correct this, bring up the "Modify Atom" widget as follows:

· {Build/Edit}  >>  {Modify Atom}  >>  select "ONLY_TYPE" >> {OK}  

 

16) You will need to specify the atoms that currently have incorrect bonding behavior.  To do this, reorient the molecule so that you have a good view of the incorrectly bound phenyl ring, click on each of the 6 Carbon atoms in the ring (they will become highlighted with little green diamonds) and then click {OK}.

 

17) You are now prompted for the type of atom you want to change each of these Carbons to.  For a description of each of the different SYBYL atom types, open up your Netscape browser and choose the {Atom Types (Tripos)} bookmark.  For now, though, you are going to want to change all of these carbon atoms to 'aromatic carbons'.  To do so, scroll down on the widget, select "C.ar" and hit {OK}.

 

18) At this point, you will have changed the first of the six C atoms.  The widget is now asking you is you want to change the second atom to "C.ar" as well (note where the green diamond is on the molecule to see which atom you're about to change).  You do wish to change this one to "C.ar" as well so click {OK}.  In a similar fashion, the widget will cycle through all of the atoms you've asked to change and will disappear when you're done.

 

19) Notice that the one N in the structure also has an incorrect bond. Change the N to type "N.am" and the carboxyl C next to it to "C.2".​

20) You are going to need to hydrogenate your structure in order to optimize it properly.  To do so:

· {Build/Edit}  >>  {Add}  >>  {Hydrogens}

 

21) You are now ready to begin an optimization, so first save your structure:

· {File}  >>  {Save As}  >>  choose "m2"  >>  call the molecule "taxol_start.mol2"  >>  {Save}

and then bring up the Minimize widget as follows:

· {Compute}  >>  {Minimize...}

 

22) This time, we are not going to use just the default settings.  We would like to try to reproduce some interesting experimental results that suggest that there is an observable attraction between one of the two phenyl rings close to each other in the taxol crystal structure with the third phenyl ring on the other side of the taxol molecule.  To try to reproduce this, we are going to:

· Change the maximum number of optimization iterations (in the widget's "Max Iterations:" text box) from 100 to 1200

 

23) We are also going to make the interaction forces more realistic as follows:

· Click the {Modify...} button next to "Energy Setup" to give you the Energy Setup widget

· On this new widget, click the button next to "Charges" that currently says {None}.  Change it to "Gasteiger-Huckel".

· Since the two phenyl rings of interest are currently more than 8.0 Angstroms away from each other, change the "NB Cutoff:" text box from "8.0000" to "15.0000".

· Click {OK}

 

24) Submit the job by clicking {OK} on the minimize widget.

 

25) When the job is done (in a few minutes), compare it to your starting structure by loading the starting structure back into SYBYL:

· {File}  >>  {Read...}  >>  select "taxol_start.mol2" from the "Files:" list and select "m1" as the molecular area  >>  click {OK}

 

26) To make it easier to compare the two structures, you might wish to recolor the starting structure:

· {View}  >>  {Color}  >>  {Atoms...}  >>  select M1 and click {All}  >>  {OK}

· Choose a color from the list (purple is often good)  >>  {OK}

 

27) At this point you should see that the minimization has moved two of the rings somewhat closer together, but you also note that they are prevented from getting too close because they're oriented perpendicular to each other. NMR studies, on the other hand, suggest that these two rings should really be quite close though, and should be oriented parallel to each other. What has gone wrong? You will find out in the next exercise. For now make a note of the final total energy for this structure for future reference (in the text window at the bottom of your screen, scroll up until you find the "Total Energy :" line.
 

EXERCISE 3:  MOLECULAR DYNAMICS USING CLASSICAL MECHANICS

In this exercise, you will look at how temperature and time-dependence affects molecular structures and will use these effects to try to solve the taxol structure problem.

 

1) Delete the molecule you just optimized:

· {Build/Edit}  >>  {Zap (Delete) Molecule}  >>  select "M2"  >>  {OK}

· Do not delete the taxol starting structure stored in M1.

 

2) Bring up the molecular dynamics widget:

· {Compute}  >>  {Dynamics}  >>  {Setup Dynamics...}

 

3) Click on the {Modify...} button on the TOP RIGHT of the widget to bring up the Interval Setting widget.

 

4) Change the simulation's "Length:" from "1000" fs to "2500" fs.

 

5) We are going to run this simulation at constant pressure, so change the "Ensemble:" setting from "NTV" to "NTP" by clicking on the {NTV}.

 

6) Change the Pressure from "0" atm to "1" atm.

 

7) Note that the temperature is set at 300 K, which is fairly close to body temperature -- leave the temperature at this value.

 

8) Change the "Initial Velocity:" setting from "Previous" to "Boltzmann" by clicking on the {Previous}

 

9) Change the "Seed:" value to "32759". Note that this number is supposed to be randomly selected, and it is often fine to let the program choose it for you, but I have selected this seed value so that you can see an interesting dynamics result in a relatively short period of time -- for other seed values, you would eventually see similar behavior, but might have to wait longer.^m>

 

10) Click {OK} to bring yourself back to the main Dynamics widget.

 

11) Note that the "Energy Setup:" settings that you specified in Ex. 2 have been carried over for you to this exercise (i.e., Charge = Gasteiger-Huckel and NB Cutoff = 15.0 Angstroms).  As a result you will not have to change any of the energetics settings this time.

 

12) Change the "Job Name:" to "taxol_md" and click {OK}.

 

13) Once the dynamics job has completed, replay the trajectory:

· {Analyze}  >>  {Dynamics}  >>  {Replay Trajectory...}  >>  select "taxol_md.his"  >>  {OK}

What do you notice about the relative positions of the phenyl rings?  You can re-replay the trajectory from different angles (by using the mouse to reorient the molecule) if you want a better view of what's going on.

In the dynamics simulation, you should see that all three phenyl rings actually come fairly close together. Also, for a little while jut before the end, the three rings seem to orient themselves roughly parallel to each other, unlike the optimized structure obtained in the previous exercise. This suggests that there may be one local minimum structure with the two closest rings largely perpendicular to each other), but there may also be also another much different local minimum structure that likes to have the rings closer to parallel. In a straight geometry optimization, the simulation algorithm only finds the first minimized structure, but sometimes when you do molecular dynamics at a warm enough temperature (300 K), you give the molecule enough energy to hop from one potential well to another.

 

14) Let's see if we can figure out which of these two local minima is lower in energy:

· {Analyze}   >>  {Dynamics}  >>  {Analyze Dynamics...}  >>  select "taxolmd.his"  >>  {OK}

· Select "M1" >> {OK}

· Click {OK} for a data reduction interval of "1"

· click {OK} to start examining data from time "0"

· change the "Data End" valie to "2500" and click {OK}

 

15) This gives you a table of energetic properties for every fifth geometry step in the entire dynamics simulation.  Unfortunately, it is hard to tell from such a large table what the lowest energy conformations are.  Fortunately, from the table widget, you can sort your results:

· {View}  >>  {Sort}  >>  select "TOT_ENG" as the primary sorting criterion and leave the others unchanged  >>  click {Sort}

You now see that the table has been rearranged according to increasing energy.  The lowest energy structures near the top of the table are all right around the 2000 fs time range -- about four fifths of the way through the simulation.  Do you recall approximately what sort of configuration the molecule was in around 4/5 of the way through the simulation when you replayed the trajectory?

 

16) If you can't remember what things looked like at that point, you can find out exactly what the molecular conformation is at a given point:

· On the table widget, click in the leftmost column of any row you would like to see the corresponding structure of

· Click the {Show RowSel} button on the table widget

In particular, look at the structure corresponding to the lowest energy (403: R2010.000) -- what can you say about the phenyl rings in this structure?

 

17) Even though the point (403: R2010.000) of the dynamics simulation is the lowest energy conformation, it is not necessarily the global minimum energy structure for the molecule.  This is because you carried out the dynamics simulation at 300 K, and at this temperature the molecule has too much energy to stay in the absolute lowest energy conformation.  To find that conformation, you should now take the structure for (403:R2010.000) and do a straight geometry optimization:

· {Compute}  >>  {Minimize...}

· Note that all optimization settings have been carried over from your last minimization job.  As a result, you don't need to change anything and can just click {OK}

What happens to the two phenyl rings during the optimization?

How does your total energy compare with what you got at the end of Ex. 2?

 

EXERCISE 4:  MOLECULAR DOCKING USING FLEXIDOCK

This exercise will look at the docking of a taxol analog in tubulin.

 

1) Delete the optimized taxol molecule from Ex. 3:

· {Build/Edit}  >>  {Zap (Delete) Molecule}

 

2) Re-read the original taxol starting structure:

· {File}  >>  {Read...}  >>  select "taxol_start.mol2"  >>  {OK}

 

3) Label all of the atoms:

· {View}  >>  {Label}  >>  {Atom Name...}  >>  select "M1" and click {All}  >>  {OK}

 

4) We are going to change an alcohol into an ether by removing H28 and  replacing it with an N-Propyl group.  Find H28, memorize where it is, and:

· {Build/Edit}  >>  {Sketch Molecule...} >> "M1:" >> {OK}

· Click {REMOVE_ATOM}  and then click on the atom that had been labeled "H28"

· Click {GROUP}  >>  click {N-PROPYL}  and then click on the "O" atom from which you just removed "H28"

· Click {ADDH}

· Click {EXIT}

 

5) Read in the tubulin crystal structure again:

· {File}  >>  {Read...}  >>  select "1JFF.pdb" and "m2:"  >>  {OK}

· {No}  --  do not center the molecule

 

6) Molecular docking typically doesn't pay attention to atoms much farther than 10 angstroms outside the active site, so in order to simplify our task, we're going to trim away unnecessary parts of the molecule:

· {Build/Edit}  >>  {Delete}  >>  {Atom...}

· Select "M2:"

We are actually going to start by telling SYBYL what part of the molecule to NOT delete:

· Click {Substructures...}  >>  choose "TA1601" (very bottom)  >>  {OK}

· Click {Sets...}  >>  enter "10.0" in the textbox labeled "R =" and click {OK}

So we have now selected all atoms within 10.0 Angstroms of the active site.  In order to retain those atoms and delete everything else:

· Click {Invert}  >>  {OK}

    

7) Bring up the FlexiDock widgets:

· {Compute}  >>  FlexiDock  >>  {Create an Input File...}

· Select "M2:"  >>  click {OK}

· Read the "No Pocket Defined" message and click {OK}

 

8) First you will need to define your active site:

· On the left widget click the {Define Pocket} button -- this brings up a new pocket-defining widget

· Select "M2:" since your active site will be part of the tubulin

· Click {Substructures...}  and once again choose "TA1601" to define the center of your active site  >>  {OK}

The "No Pocket Defined" message in 7) told us that FlexiDock automatically extends the active site model to include anything within 4 Angstroms of your defined center.  This isn't very big, so lets extend it a bit:

· Click {Sets...}  >>  choose R = {3.0}  >>  {OK}

which means that your active site will now extend 3.0 + 4.0 = 7.0 Angstroms.

· Click {OK} to complete the pocket selection.

 

9) The active site, in which we're going to dock the taxol analog defined in steps 2-4), still contains the original taxol structure however.  To extract it:

· On the right widget, click {Extract Ligand}  >>  select "M2:"  >>  click {Substructures} 

· In the Substructures widget, select "TA1601"  >>  {OK}

· Back in the extraction widget, click {OK}

· Dump it into "M3:"  >>  {OK}

 

10) To avoid confusion, you should probably get rid of this extracted taxol entirely.  From back in the main SYBYL menu:

· {Build/Edit}  >>  {Zap (Delete) Molecule}  >>  select "M3:"  >>  {OK}

 

11) You will now need to revise the ligand location:

· In the right widget, change the "Ligand Location" button from {M3} to {M1}.

 

12) Click the {Remove} button next to the "Water has been removed" checklist item.

 

13) You have already checked the atom types in Ex. 2 , so click the small check button next to "Atom types have been checked"

 

14) Click {Add} next to "Hydrogens have been added"

 

15) Click {Add} next to "Charges have been added"

 

16) Define your ligand's rotatable bonds:

· Click {Choose} next to "Rotatable bonds have been set"

· Click {Add} below "Bonds to Flex in Ligand"

· Select "M1:"

We now have to click on all bonds in the ligand that we want to be able to flex.  If we hit the {All} button, we could have the entire molecule be flexible, but in the interests of time, we're going to limit this a bit, by choosing specific bonds from the molecule.  This is done by clicking on each of the two atoms between which the bonds is formed (make sure the bond turns green).   Do not bother selecting bonds that:

* Are part of a ring structure (the rings in taxol are quite rigid)

* Include a hydrogen (H flexion doesn't contribute much  to docking potentials)

* Include a terminal methyl (quite spherical and nonreactive)

 

17) For a high quality docking simulation we would probably want to allow part of the active site to flex as well.  In this case, we are just doing a simple screening study, though, and will leave it rigid.  Click {OK} on the "Pick Rotatable Bonds" widget.

 

18) You can specify expected H-bonds in the ligand-receptor complex by pairing up expected H-bond donors with H-bond acceptors.  For a simple screening job it is fine to go without explicit H-bond input, so just click the check box next to "H-bond sites are marked".

 

19) Because our taxol analog has been created by modifying a taxol cut directly out of the crystal structure, and because we have been careful not to recenter molecules when we loaded them in, the ligand is already prepositioned close to the cavity and we can just check off the small check box next to "Ligand is pre-positioned in cavity".  Otherwise we would have had to hand dock a guess structure.

 

20) Specify the job name as "taxol+tubulin"  (written as one word, and entered in the "Name to Use ...." text box)

 

21) Click {Write Out FlexiDock Input File}

 

22) Submit the job:

· Click {FlexiDock It}

· Click {OK} for the file and parameter file settings

· Click {OK} for the random seed

· Click {OK} to leave the "Maximum number of generations to allow" (genetic algorithm setting) unchanged

The job may take upwards of five minutes to run.

     

23) Often times your docking run will yield multiple candidate structures, so SYBYL collates the results into a molecular database.  In this case, the program probably only found one or two good structures, but let's bring it up in database format anyway:

· {Build/Edit}  >>  {Zap (Delete) Molecule}  >>  {All}  >>  {OK}

· {File}  >>  {Molecular Spread Sheet}  >> select "New"   >>  select "Database"  >>  {OK}

· Select "taxol+tubulin.mdb"  >> {Open}

 

24) You have only two rows, of which the first is your starting structure and the second is your successfully docked structure.  If you want to display one of the structures:

· On the table widget, click in the leftmost column of the row you want to display a structure for

· Click the {Show RowSel} button on the table widget

If you want to load both of the two structures into your window at the same time:

· Main SYBYL Window  >>  {File}  >>  {Read...}  >>  click the "taxol+tubulin.mdb"  "Sub-Directory"  >>  select the appropriate *.mol2 file  >>  {OK}

· Use the {View}  >>  {Color}  >>  {Color Atoms....}  options if you need help distinguishing the different structures


CONTACT
David Johnson
emaildkjohnson@ku.edu
Access to more than 40 computational chemistry software programs and databases
High-performance computational tools accelerate drug discovery and minimize costs
Analyze complex, multidimensional data sets to rapidly generate biological insights
KU Today