Molecular Modeling Team

Details of the sietraj.tcl script

The sietraj.tcl script can be treated as a black box for the SIE calculation. However, for those interested, the entire sietraj.tcl script is shown below with some comments giving a general overview of the various sections. (The source code for this script can be found in $BRIMM_ROOT/apps/sietraj directory.) The script makes heavy use of functions defined in the loaded Tcl libraries. The arguments passed to the script are documented in the code. The start and end atom numbers for the target and ligand can be found by generating a pdb file for the complex and looking for the TER records. For most cases this script does not need to be edited.
#!/usr/local/bin/tclsh # Tcl script for calculating solvated interaction energies (SIE) for an # Amber-generated MD trajectory. # The arguments to this script are: # * prmtop : prmtop file # * traj : trajectory file # * istart : first snapshot to include in the SIE calculation # * iend : last snapshot < iend # * inc : Every inc snapshots starting with istart will be included # * targ_range : Range of atom numbers comprising the target. # Must be no spaces around the -, e.g., 1-1947. # If the numbering is not contiguous, enter all ranges separated # by commas with no spaces: 1-1947,2000-2034,2100-2234. # Same thing applies to lig_range. # * lig_range : Range of atom numbers comprising the ligand. # * outfile : Output file containing SIE energy terms # * dbflag : If set to 1, output to an sqlite3 database instead. # * dbtag : Tag the results in the database as belonging to one data set. # Allows the database to be used for several data sets. # * mut2ala : If not zero, indicates the residue number to be mutated to Ala. #$Id: sietraj.tcl 350 2012-03-20 16:11:15Z Enrico Purisima $ scan $argv %s%s%d%d%d%s%s%s%d%d%d prmtop traj istart iend inc targ_range lig_range \ outfile dbflag dbtag mut2ala set targ_range [split $targ_range ,-] set lig_range [split $lig_range ,-]
The following loads the runtime library as well as various supporting Tcl functions.
source $env(BRIMM_ROOT)/scripts/Britcl_startup.tcl source $env(BRIMM_ROOT)/scripts/amber_util.tcl set amb2syb $env(BRIMM_ROOT)/scripts/amb2syb.tcl
Next we set the parameters for the solvated interaction energy, which is calculated as

α * (EvdW + Ecoul + ERF + γ * Δ SA) + const

ERF is the reaction field energy calculated using BRI BEM (Purisima, 1998; Purisima & Nilar, 1995). The atomic radii (AMBER vdW radii) will be scaled by 1.1. These parameters were obtained by fitting to a database of 99 protein-ligand complexes (Naim et al., 2007). The solute and solvent dielectrics are 2.25 and 78.5, respectively.
# Optimized parameters for Dout = 78.5 # If you change these, make sure to update averages.awk as well. set alpha 0.1047577 set gamma 0.01289386 set const -2.89 set rscale 1.1 bemvars configure -Din 2.25 bemvars configure -Dout 78.5 bemvars configure -normalize 1
The molecular surface is generated using the marching tetrahedra algorithm (Chan & Purisima, 1998) and variable solvent probe radii (Bhat & Purisima,2006). solvRadAnionic and solvRadCoul are for carboxylates and polar groups, respectively. solvRadVdw is for all other atom types. For very large systems, you might get an error message that MarTet_MAXN2, MarTet_MAXN3 or MarTet_TMP_SIZE is too small. Increase these values as necessary. The values here should work for most systems.
bemvars configure -solvRadVdw 1.6 bemvars configure -solvRadCoul 1.29 bemvars configure -solvRadAnionic 1.08 bemvars configure -minMarTetArea 0.000001 bemvars configure -marTetPrintDetail 0 bemvars configure -bemSurfPrintDetail 0 set MarTet_MAXN2 40 set MarTet_MAXN3 100 set MarTet_TMP_SIZE 200000 bemvars configure -useCurrentSurface 0 bemvars configure -usePrevEnormal 0 bemvars configure -useExplicitBoxCenter 0
The molecule is generated from the prmtop file and the atom types and force field parameters are assigned. The TARGET and LIGAND sets are defined. The SIE is calculated from the rigid infinite separation of these two sets of atoms.
# Create the molecule, define the TARGET and LIGAND sets and assign force field # parameters. Molecule m1 m1 InitFromPrmtop $prmtop m1 DefineSetByRange 1 "TARGET" [lindex $targ_range 0] [lindex $targ_range 1] if {[llength $targ_range] > 2} { for {set i 2} {$i < [llength $targ_range]} {incr i 2} { m1 AppendRangeToSet "TARGET" [lindex $targ_range $i] [lindex $targ_range [expr $i + 1]] } } m1 DefineSetByRange 1 "LIGAND" [lindex $lig_range 0] [lindex $lig_range 1] if {[llength $lig_range] > 2} { for {set i 2} {$i < [llength $lig_range]} {incr i 2} { m1 AppendRangeToSet "LIGAND" [lindex $lig_range $i] [lindex $lig_range [expr $i + 1]] } } SetSybtypeFromPrmtop m1 $prmtop $amb2syb ReadChargesFromPrmtop m1 $prmtop SetEmbeddedAmberAtypeFromPrmtop m1 $prmtop amber AssignBemRadiusFromAmbtype m1 amber $rscale
If a virtual alanine mutation is desired, a second molecule (m2) is created that will contain the complex with the alanine mutation.
# Mutate a residue to Ala? if {$mut2ala > 0} { array set alaCharges {N -0.4157 H 0.2719 CA 0.0337 HA 0.0823 CB -0.1825 HB1 0.0603 HB2 0.0603 HB3 0.0603 C 0.5973 O -0.5679} array set alaSybtypes {N N.am H H CA C.3 HA H CB C.3 HB1 H HB2 H HB3 H C C.2 O O.2} SetAtomNameFromPrmtop m1 $prmtop SetAmberAtypeNameFromPrmtop m1 $prmtop amber atypeNameId set HCType $atypeNameId(HC) set HType $atypeNameId(H) # Flag atoms to discard/change in mutated residue set natoms [m1 cget -natoms] set atomflag [NewInt [expr $natoms + 1]] set atom_map [NewInt [expr $natoms + 1]] # Create modified molecule for SIE calculation Molecule m2 m2 InitAtoms $natoms puts "Res # $mut2ala ([FixCgetString [[m1 GetSubstructure $mut2ala] cget -name]]) will be mutated to Ala" AlaMutForSIE m1 m2 $mut2ala atomflag atom_map $HCType $HType alaCharges alaSybtypes AssignBemRadiusFromAmbtype m2 amber $rscale UpdateTargLigSet m1 m2 $atomflag $atom_map }
The starting snapshot frame is located. The output file is opened. There is also the option to write to an sqlite database instead of a regular file.
if {[GotoFrameInTrj m1 $traj $istart fp_trj] != 0} { puts "Error in setting frame position" exit -1 } if {$dbflag == 0} { set fp_out [open $outfile "w"] } else { if {$tcl_platform(wordSize) == 8} { load $env(BRIMM_ROOT)/lib/tclsqlite-64.so } else { load $env(BRIMM_ROOT)/lib/tclsqlite-32.so } proc db_busy {a} { set rnd [expr 1000 * rand()] set ms [format "%0.f" $rnd] #puts "Database busy ... Will try again in $ms ms." after [format "%0.f" $rnd] ; return 0 } sqlite3 db $outfile db busy db_busy db eval {create table if not exists sie (id integer primary key autoincrement, batch integer, iframe integer, deltaG double, E_rf double,rfc double, rft double, rfl double, vdw double, elec double, darea double)} }
The SIE is calculated for the chosen frames. If this is a virtual alanine mutation, the coordinates for the trajectory will be used to update m2 rather than m1. The SIE calculation will then be performed on m2 rather than m1. The results are printed on the screen and to a file. All energies are in kcal/mol. The change in area is in square Angstroms.
if {$mut2ala == 0} { SIE m1 bemvars amber TARGET LIGAND vdw coul rf_complex rf_target rf_ligand darea } else { SIE m2 bemvars amber TARGET LIGAND vdw coul rf_complex rf_target rf_ligand darea } set E_rf [expr $rf_complex - $rf_target - $rf_ligand] set deltaG [expr $alpha * ($vdw + $coul + $E_rf + $gamma * $darea) + $const] puts [format "%d %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f" \ $iframe $deltaG $E_rf $rf_complex $rf_target $rf_ligand $vdw $coul $darea] if {$dbflag == 0} { puts $fp_out "$iframe $deltaG $E_rf $rf_complex $rf_target $rf_ligand $vdw $coul $darea" flush $fp_out } else { db eval {insert into sie (batch,iframe,deltaG,E_rf,rfc,rft,rfl,vdw,elec,darea) values ($dbtag, $iframe, $deltaG, $E_rf, $rf_complex, $rf_target, $rf_ligand, $vdw, $coul, $darea)} } # Position file pointer to next frame if {[expr $iend - $iframe] >= $inc} { SkipFramesInTrj m1 $fp_trj [expr $inc - 1] } } if {$dbflag == 0} { close $fp_out } close $fp_trj
References
  1. Bhat, S.; Purisima, E. O. Molecular surface generation using a variable-radius solvent probe. Proteins: Structure, Function, and Bioinformatics 2006, 62, 244-261.
  2. Chan, S. L.; Purisima, E. O. Molecular Surface Generation Using Marching Tetrahedra. J. Comput. Chem. 1998, 19, 1268-1277.
  3. Nam, M.; Bhat, S.; Rankin, K. N.; Dennis, S.; Chowdhury, S. F.; Siddiqi, I.; Drabik, P.; Sulea, T.; Bayly, C.; Jakalian, A.; Purisima, E. O. Solvated interaction energy (SIE) for scoring protein-ligand binding affinities. 1. Exploring the parameter space. J. Chem. Inf. Model. 2007, 47, 122-133.
  4. Purisima, E. O. Fast Summation Boundary Element Method for Calculating Solvation Free Energies of Macromolecules. J. Comput. Chem. 1998, 19, 1494-1504.
  5. Purisima, E. O.; Nilar, S. H. A simple yet accurate boundary element method for continuum dielectric calculations. J. Comput. Chem. 1995, 16, 681-689.