Molecular docking with VEGA ZZ and Fred

1 Introduction
2 What's you need
3 FTase download
4 Protein preparation
5 Creation of the input files for NAMD
6 Run the NAMD minimization
7 Protein refinement
8 Ligand build and conformational search
9 Definition of the docking region
10 Molecular docking
11 Analysis of the docking results
12 Tips & Ticks
 

1 Introduction

VEGA ZZ, Fred and NAMD are very efficient tools to perform a complete molecular docking analysis. This tutorial explains in details how was performed the docking between the farnesyltransferase enzyme (FTase) and a set of inhibitors as published in the paper: Bolchi C., Pallavicini M., Rusconi C., Diomede L., Ferri N., Corsini A., Fumagalli L., Pedretti A., Vistoli G., Valoti E., "Peptidomimetic inhibitors of farnesyltransferase with high in vitro activity and significant cellular potency", Bioorg. Med. Chem. Lett., 17(22), 6192-6 (2007).

2 What's you need

• VEGA ZZ release 2.2.0 or greater (click here for the software setup).
• ISIS/Draw (free downloadable for academic non-profit use at Elsevier MDL Web site) or alternatively SketchEl included in VEGA ZZ package.
• Fred (available for free for academic non-profit use at OpenEye Scientific Software).
• NAMD for Windows (click here to download it, click here to install it).
• Crystal structure of the human FTase complexed with the farnesyl diphosphate (FPP) and the peptidomimetic inhibitor L739,750. available at Protein Data Bank (PDB ID 1JCQ).

3. FTase download

You can download the FTase complex (1JCQ) structure using the PDB Web interface or the tool integrated in VEGA ZZ:

• Start VEGA ZZ and select the File PDB download menu item.
• Put 1JCQ in the PDB Id field and click Download. At the end of the download, the protein structure will be shown in the workspace (for more information, click here).
• Normalize the coordinates in order to translate the protein at the origin of the Cartesian axis (Edit Coordinates Normalize).
• Save the molecule (File Save As) with 1JCQ file name. It's strongly recommended the use of the IFF/RIFF file format because it's able to keep the maximum number of information (e.g. atom types, charges, bond orders, etc).

4 Protein preparation

The zinc ion must be unbonded to all atoms in order to assign  the correct atom types. It results bonded to the thiolate group of the L739,750, to the ASP 297 B carboxylate, to the CYS 299 B thiolate group, to the ring nitrogen NE2 of the HIS 362 B. To remove the wrong bonds, follow these steps:

• To make easier this operation, you must select the residues around the zinc ion: select the zinc (View Select Custom), clicking ZN* in Atom column and + button. Choose the Proximity tab, click on the zinc in the main window, select Residues in What box, Atom in Around box and Normal in Selection box. Put 3 as Radius and click + button. The residues in a 3 Å sphere centred on the zinc will be shown.
• In the main menu, select Edit Remove Bonds and choose One bond in What ? box. In the main window, click on the atom pairs of each bond to delete (four pairs) or directly on the bonds. Alternatively you can right-click the bond and choose Bond Remove in the context menu.
• Add the hydrogens (Edit Add Hydrogens), selecting Protein in Molecule type box to enable an extra check for atom hybridizations, Residue end in Position of hydrogens box and check Use IUPAC atom nomenclature. Finally, click Add to place the hydrogens.  Warning messages will be shown in the console to inform you that some atoms have an unusual geometry and the extra check has corrected the atom types.

Adding the hydrogens, the thiolate groups of the L739,750 and the CYS 299 B will be protonated, but they must be without hydrogens:

• Show the residues around the zinc as explained above.
• In the main menu, select Edit Remove Atom and click in the main window on the hydrogens bonded to the two sulphur. They will be removed. Alternatively, you can right-click the atom and select Atom Remove in the context menu.
•  The secondary aminic group of the L739,750 is protonated, and so click on  one hydrogen bonded to it.

The carbonyl of both C-terminal residues (HIS 367 A and GLU 424 B) are in aldehydic form and must be edited to carboxylate:

• Select View Select Custom,  choose Atoms tab, click Reset button to empty the Selection string, click HIS, 367 and A in Residue, Number and Chain columns, choose Normal in Selection box and finally click + button. The HIS 367 A will be shown.
• In the main menu, choose Edit Change Atom/residue/chain, click on the aldehydic hydrogen.
• In Element field of Edit window, change H to O and press Apply button.
• Repeat the same operation for the GLU 424 B.
• Close the Edit window.

The sucrose and the acetate are far to the binding site and so can be removed:

• Edit Remove Residue, click on ACY 3002 and press the Remove button.
• Repeat the same procedure for the SUC 3010 residue.

Select all atoms (View Select All) and remove the labels (View Label atom Off).
In order to optimize the crystal structure of the complex, it will performed an energy minimization, using the CHARMM22_LIG force field keeping fixed the protein backbone. The CHARMM22_LIG contains more parameters than the standard CHARMM22_PROT for proteins.

• Fix the atom types and the charges (Calculate Charge & Pot.), checking Force field and Charges and selecting CHARMM22_LIG and Gasteiger. Click Fix button. The total charge is -25.

Before the energy minimization, the atom constraints must be applied. We want to fix the protein backbone and the zinc ion, in order to avoid too many modifications in the experimental structure.

• Select the protein backbone (View Select Backbone).
• Select the zinc ion, showing the selection window (View Select Custom), clicking the Reset button, selecting ZN* in the Atom column, choosing Additive in the Selection box and clicking + button.
• In the main menu, choose Edit Coordinates Constraints, select Fix in Mode box and Visible atoms in Selection box. Finally, click Apply button. The fixed atoms will be painted in blue and the free atoms in green.
• Save the molecule as 1JCQ.iff in IFF format.

5 Creation of the input files for NAMD

NAMD requires some files: the PDB and the PSF file of the molecule, one or more parameter files and a command file defining the condition of the calculation.

• Save the molecule in PDB 2.2 format (File -> Save As...) with the 1JCQ.pdb file name, unchecking Connectivity and checking Constraints in the Options box. In this way the connectivity will not be saved and the atom constraints will be included in the B column of the PDB file.
• Create the topological matrix saving the molecule in the X-Plor PSF format (1JCQ.psf file name). Optionally, it's possible to check if all force field parameters are available, selecting the force field name in the Force field param. box. The CHARMM22_LIG force field was used in the atom type attribution and so you must select the same force field.
• Copy in the directory where are placed the 1JCQ.pdb and 1JCQ.psf files, the par_all22_na.inp, par_all22_prot.inp and par_all22_vega.inp parameter files that are in the ...\VEGA ZZ\Data\Parameters directory. The par_all22_na.inp and par_all22_vega.inp are required because it contains the parameters for the ligands.
• With a text editor (e.g. Notepad) create the input file with the following commands (copy & paste them):

numsteps                50000
minimization            on
dielectric              1.0

coordinates             1JCQ.pdb
outputname              1JCQ
outputEnergies          5000
binaryoutput            no
DCDFreq                 5000
restartFreq             5000

structure               1JCQ.psf
paraTypeCharmm          on
parameters              par_all22_na.inp
parameters              par_all22_prot.inp
parameters              par_all22_vega.inp
exclude                 scaled1-4
1-4scaling              1.0
switching               on
switchdist              8.0
cutoff                  12.0
pairlistdist            13.5
margin                  0.0
stepspercycle           20

fixedAtoms		on
fixedAtomsCol		B

Save the file with the 1JCQ_min.namd name. This file allows to perform a 50000 steps conjugate gradients minimization, saving the output (coordinates and restart files) every 5000 iterations. For more information about the parameters, please consult the NAMD User Guide.

6 Run the NAMD minimization

• Open the VEGA console (Start -> VEGA ZZ -> VEGA console).
• Go to inside your working directory with the cd command.
• In the console type:

namd2 1JCQ_min.namd > 1JCQ_min.out

If you have more than one CPU installed, you can speed-up  the calculation specifying the total number of CPUs:

namd2 +p2 1JCQ_min.namd > 1JCQ_min.out

In this case the PC has two CPUs and both are used for the calculation (+p2 option). The dual core (e.g. Athlon X2, Pentium D, Core 2 Duo, etc) and the Pentium 4 with hyperthreading should use the +p2 option.

At the end of the calculation two files are generated: the 1JCQ.dcd (trajectory file) and the 1JCQ.out. The first one is a binary file that can't be opened with a text editor. It contains the atom coordinates of each saved frame (10 frames, because one frame every 5000 was saved). The second one is a text files containing the output messages generated by NAMD and the energy information. We need to save the lowest energy structure:

• Open the 1JCQ.dcd file with VEGA ZZ (File -> Open). To open a DCD trajectory a molecule file is required (e.g. in PDB or IFF format) with the same name. You shouldn't have any problem because you have the 1JCQ.dcd and the 1JCQ.iff file in the same directory.
  The molecule will be shown and the Trajectory analysis dialog will be opened.
• If you saved the 1JCQ.out file, the Energy Graph button will be active. Clicking it, you can see the energy behaviour during the calculation. The energy go down as you should expect by an minimization.
• Clicking the Last button in the Trajectory analysis window, the lowest energy conformation is selected (see the workspace).
• Save the best conformation (File -> Save As) in IFF format (1JCQ_min.iff).

7 Protein refinement

The protein isn't ready for the docking because the ligands (FPP and L739,750) are already inside and should be removed to create enough space in the active site to dock the new ligands. In order to enlarge a little bit the cavity, the co-crystallization water molecules included in a sphere of 8 Å radius centred at the zinc ion coordinates, will be removed (1006, 1165, 1252, 1425, 1674, 1676, 1677,1678, 1680, 1684, 1686) because in a biological environment, a ligand could replace a water molecule and in the molecular docking is impossible to simulate it.

• Color the molecule by atom (View -> Color -> By atom)
• Select the residues inside a sphere of 8 Å radius centred on the zinc ion (see above).
• Edit -> Remove -> Residue, check Keep the window opened, double click on FPP and 739 in the residue list.
• Remove all visible water molecules clicking on them in the main window.
• Display all atoms (View -> Select -> All)
• Save the molecule (1QBQ_dock) in IFF/RIFF and PDB 2.2 formats. The PDB file is required by Fred to perform the docking.

The protein is now ready for the docking.

8 Ligand build and conformational search

In this section will be explained how one ligand (compound 4) reported in the paper was built and how was performed the conformational analysis. This step is required by the docking software  to consider the ligand flexibility because both ligand and protein are kept rigid.

• Create an empty workspace in VEGA ZZ (File -> New).
• Start ISIS/Draw, selecting Edit -> ISIS/Draw.
• Draw the following molecule:

• In the ISIS/Draw main menu, select VEGA ZZ -> Send to VEGA. The molecule will be sent to VEGA ZZ and automatically converted from 2D to 3D.
• Check both chiral centres because the conversion is unable to fix them: the former, inside the pyridodioxane ring, must be R and the latter (alpha carbon of the methionine group) must be S.
• If both centres are wrong (S, R) you can change them inverting the Z atom coordinates (Edit -> Coordinates -> Invert Z). This procedure allows to obtain the enantiomer.
• If only one chiral center is wrong, swap two bonds selecting Edit -> Change -> Swap bonds and click two atoms bonded to the chiral center.

The subsequent conformational analysis will be performed by AMMP applying the SP4 force field parameters that are dynamically assigned starting from the bond types (single, partial double, double and triple). For this reason, the bond order must be correctly assigned, but looking the ligand structure you can see that the methionine carboxylate group has the two oxygens bonded by single and double bonds instead of two partial double bonds.

• In the main menu, select Edit -> Change -> Bond type, choose One bond and Partial double in Wath ? and Bond type boxes.
• Click the carboxylate carbon and one oxygen. Repeat the operation for the second oxygen.
• Click the Done button in the Bonds window.

At this step, the ligand structure need an energy minimization in order to refine the 3D structure. The assignment of charges and potentials is not required because it was done automatically during the 3D conversion.

• Select Calculate -> Ammp -> Minimization, in the Ammp minimization window click the Default button, change Toler to 0,01 and Conjugate gradients as minimization algorithm. To start the calculation, press the Run button. The minimization will converge in few steps because the structure was previously optimized by the 3D conversion.
• Remove the atom labels (View -> Label atom -> Off).
• Save the molecule in IFF/RIFF format as ligand4.iff.

Now we perform the conformational search based on the Boltzman jump algorithm because the systematic search is too expensive in time terms due to the nine flexible torsions of the molecule.

• Open the Ammp conformational search dialog window (Calculate -> Ammp -> Conformational search) and press the Default button.
• In the Search parameters box, select Boltzman jump as Method, put 5000 Steps, 60 degrees as torsion RMSD and 2000 K as temperature (Temp.).
• Check Minimize all conformations, put 300 Steps and 0,01 Toler.

You must define the torsion angles that will be changed during the conformational analysis:

• Press the Edit torsion button in the Search parameters box: the Selection tools will be shown.
• In that dialog window, select Edit -> All flexible torsions: the flexible torsion angles are automatically detected and shown in the main window.
• Click Done to confirm the torsions.
• Check Trajectory and Energy to save all 5000 conformations in the compound4.dcd file and the energy of each conformation in the compound4.ene file.
• To start the calculation, click the Run button.

The Boltzmann jump approach generates the conformations randomly without checking if they are similar. For this reason, the cluster analysis must be performed:

• In the main menu, select Calculate -> Analysis.
• In the Trajectory analysis window, click the open button,  select and open the compound4.dcd file.
• Select the Cluster tab, choose Torsion RMSD in the Method box, check Save energy, uncheck Active atoms only and put 90 degrees in the RMSD field.
• Press the Ok button and save the new trajectory as compound4_clust in Mol2 multi model format (this format is required by Fred).
• The new trajectory file will contain about 50 conformations, the best one for each cluster.
• Go to the directory containing that file and rename if from compound4_clust.ml2 to compound4_clust.mol2.

WARNING:
Due to the stochastic nature of the Boltzmann jump approach, running two conformational analysis of the same molecule, you could obtain different results (e.g. different number of clusters).

9 Definition of the docking region

As explained above, the docking study will be performed by OpenEye's Fred. This software requires a special file called box defining the FTase space region in which the ligand will be docked. We assume that the ligand, developed as antagonist, could occupy the catalytic pocket in which the FPP and the L739,750 are placed in the crystal structure. For this reason, the box can be defined by FPP, L739,750 atoms and the zinc ions (known to play an important role in the catalysis). To increment the docking possibilities, the box will be enlarged of 8 Å in all X, Y, Z dimensions (see the next section).

• Clear the current workspace, selecting File -> New.
• Open the 1JCQ_min.iff file (File -> Open).
• Using the Atom selection tool (View -> Select -> Custom), select the zinc ion, choosing the Atoms tab, pressing the Reset button, clicking ZN* in the Atom column, choosing Normal in the Selection box and finally clicking the + button.
• Press the Reset button, select 739 in the Residue column, choose Additive in the Selection box and click the + button.
• Repeat the previous step selecting FPP in the Residue column.
• Remove the invisible atoms (Edit -> Remove -> Invisible atoms).
• Save the molecule as box in Mol2 file format.
• Go to the directory containing that file and rename if from box.ml2 to box.mol2.

10 Molecular docking

To run the docking, you must have the following files in the same directory:

WARNING:
You should have the OpenEye installation directory in the search path otherwise you must type the full path to run Fred.

• Open the command prompt.
• Change the current directory to the folder containing the docking files (cd command).
• Type this command and press return:

fred2 -pro 1JCQ_dock.pdb -box box.mol2 -addbox 8 -dbase compound4_clust.mol2 -oformat mol2
      -conftest none -prefix 1JCQ-compound4 -refine lig_mmff

where:

Option	Argument	Description
-pro	1JCQ_dock.pdb	File name of the target structure in which the ligand will be docked.
-box	box.mol2	File name of the box file defining the target region in which the ligand will be docked.
-addbox	8	Expand the box size of 8 Å for each dimension.
-dbase	compound4_clust.mol2	Database containing all structures/conformations that will be docked.
-oformat	mol2	Save the docked structures in mol2 multi model files.
-conftest	none	Disable the database check to find redundant structures. If you don't specify this switch, one conformation only is docked.
-prefix	1JCQ-compound4	Prefix used naming the output files.
-refine	mmff	Minimize each complex using the MMFF force field.

For a description of the Fred theory and for a more exhaustive description of the command parameters, consult the software documentation.

• Wait the end of the calculation.

Fred generates a *_docked.mol2 and a *_scores.txt file for each scoring function:

1JCQ-compound4_chemgauss2_docked.mol2
1JCQ-compound4_chemgauss2_scores.txt

1JCQ-compound4_chemscore_docked.mol2
1JCQ-compound4_chemscore_scores.txt

1JCQ-compound4_consensus_docked.mol2
1JCQ-compound4_consensus_scores.txt

1JCQ-compound4_plp_docked.mol2
1JCQ-compound4_plp_scores.txt

1JCQ-compound4_screenscore_docked.mol2
1JCQ-compound4_screenscore_scores.txt

1JCQ-compound4_shapegauss_docked.mol2
1JCQ-compound4_shapegauss_scores.txt

The .mol2 files contain the ligand poses with the same order of the corresponding .txt score files starting from the best to the worst complex.

11 Analysis of the docking results

The analysis of the docking results will be performed by VEGA ZZ. As first step, you need to choose the docking pose and to assembly the complex because the mol2 files contains the ligand only without the FTase structure.

• As an example, we want analyze the best ranked complex by the consensus score.
• In VEGA ZZ, open the 1JCQ_dock.iff file.
• In a new workspace open the 1JCQ-compound4_consensus_docked.mol2 (click New workspace in the Molecule placing window).
• In the Trajectory analysis window, you can select the pose moving the slider or typing the pose number in the Frame number field. The best consensus pose is the first one, and so you must select it.

To put the ligand inside the FTase structure you must copy & paste it.

• In the main menu, select Edit -> Copy.
• Change the current workspace clicking on the 1JCQ_dock label at the bottom of the 3D display area.
• Select Edit -> Paste: the ligand will be inserted in the protein structure.

The view is a little bit chaotic because all atoms are shown. To simplify it, you could select the residues closer to the ligand:

• Show the ligand alone, selecting View -> Select -> Molecule, click on the last molecule in the list of the Select molecule(s) window, press the Select button.
• Select View -> Select -> Custom, choose the Proximity tab, click on a ligand atom in the main window, select Residues in the What box, Molecule in the Around box and Normal in the Selection box. Put 4 as Radius and click the + button: the residues with a distance lower than 4 Å will be shown.
• Look and find the best interacting residues.

Another analysis possibility, is the use of the VEGA ZZ Evaluation of interaction tool, that requires the atomic charges and the CVFF force field correctly assigned:

• Fix the atom charges (Gasteiger) and the atom types (CVFF) as explained in the Protein preparation section.
• Select the ligand as explained above.
• In the main menu, select Calculate -> Interactions.
• Click an an atom ligand to put the it in the Residue/ligand field.
• Change the dielectric constant to 20 in order to non overestimate the electrostatic interactions. Please remember that the dielectric constant inside a protein is unknown and  so it may possible that you need more tests to adjust that value.
• The Threshold field allow to number of selected residues. A good value is 3 %, meaning that the residues contributing more than 3 % in the computation of the interaction energy are shown.
• In theVisualization box, check Ligand, Best residues, Worse residues and Water.
• In the Color box, check Ligand  by atom and By residue energy.
• Pressing the Evaluate button, the best (attractive, cyan) and the worse (repulsive, red) residues are shown and in the console, are reported the interaction energies for each residue (firstly the attractive ones).

12 Tips & Tricks

• You should repeat the interaction analysis more times adjusting the parameters to obtain a enough clear view.
• Remember that saving the complex in IFF/RIFF format, the atom selection, the colors are kept and the world transformations (rotation, translation and scale). This is useful to save the current view.
• Repeat the analysis for more than one pose (usually about five) because is not sure that the first ranked pose is the best one also.
• Repeat the analysis for more than one scoring function. This is required at the first docking in order to choose the more sensible function for your specific problem.
• If you have  a reasonable number of ligands), minimize the complexes with NAMD using the same procedure explained in the 5 and 6 sections of this tutorial.
• You could consider the possibility to apply atom constraints during the minimization to preserve partially the crystal structure, but if you want all atoms free to move, remember to delete the last two lines in the NAMD input file.