Molecular Modeling Team

sietraj

sietraj is a collection of scripts for calculating binding free energies from Amber-generated molecular dynamics trajectories of protein-ligand or protein-protein complexes. The scripts make use of a Tcl runtime library (BRIMM) that provides Tcl-accessible functions for analyzing molecular structures and carrying out energy and solvation calculations. The library is a set of C++ classes and functions with a Tcl wrapper generated by Swig.

sietraj is an alternative to the MM-PBSA software provided by the AMBER distribution. Solvated interaction energies (SIE) are calculated using parameters that have been fitted to reproduce binding free energies of a data set of 99 protein-ligand complexes (Naim et al., 2007). Note that the underlying physics behind the electrostatic components of SIE and MM-PBSA is the same. Hence, all the caveats of implicit solvation apply. The main differences lie in the choice of parameters, surface generation method (Bhat & Purisima, 2007; Chan & Purisima, 1998) and Poisson solver (Purisima & Nilar, 1995; Purisima, 1998).

Click here for download and installation instructions.

New feature
Virtual alanine mutations are now possible. (See below for the new command line option, -mut2ala.)

Sample run 1
Two sample prmtop and short trajectory files can be found in $BRIMM_ROOT/apps/sietraj/examples. (In typing the commands, replace $BRIMM_ROOT with the installation directory for the software.) If you add $BRIMM_ROOT/bin to your search path, you will not need to use the full pathname in the future. bcl2.trj is a trajectory file with 11 snapshots with all water molecules stripped away. This is a protein-ligand complex with 2611 protein atoms and 43 ligand atoms.

sietraj is invoked as shown below. The output is displayed on the screen as well as written to a file. The screen output is rounded to 3 decimal places. The output file, /tmp/bcl2.out, retains full precision. The output columns are: snapshot frame number, binding free energy, reaction field energy (binding), reaction field energy (complex), reaction field energy (target), reaction field energy (ligand), intermolecular vdW, intermolecular coulomb, change in surface area. Energies are in kcal/mol and areas are in 2.
prompt% cd $BRIMM_ROOT/apps/sietraj/examples prompt% $BRIMM_ROOT/bin/sietraj -pt bcl2.prmtop -trj bcl2.trj -sf 1 -ef 11 -inc 4 \ -tr 1-2611 -lr 2612-2654 -o /tmp/bcl2.out -sie Loaded Britcl-v2.3c-64.so Britcl_functions.tcl 338 2011-01-12 1 -5.727 9.775 -1274.224 -1278.546 -5.453 -29.364 -2.066 -420.831 5 -5.591 8.220 -1285.478 -1288.667 -5.031 -26.969 -1.747 -409.848 9 -5.931 9.110 -1274.969 -1279.056 -5.023 -30.162 -2.219 -446.502
The arguments given to sietraj are explained below. The order does not matter. 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 carrying out an SIE calculation: -pt prmtop_file (prmtop file to use) -trj trajectory_file (trajectory coordinate file) -sf start_frame (First snapshot to include) -ef end_frame (Last snapshot ≤ end_frame) -inc increment (Snapshot increment) -tr target_range (Range of atom numbers comprising the target, separated by a "-", e.g., 1-1947. If the numbering is not contiguous, enter all ranges separated with a comma with no spaces, e.g., 1-1947,2000-2034,2100-2234.) -lr ligand_range (Range of atom numbers comprising the ligand. Same format as target_range.) -o output_file (Output file containing SIE energy terms.) -sie (Carry out SIE calculation.) For carrying out a virtual alanine mutation: -mut2ala resnum (Mutate residue resnum to Ala on the fly.) For calculating averages from the SIE output: -ave filename (Process the output of the SIE calculation. filename is the output file from the SIE calculation.)
The output file can be loaded into a spreadsheet for processing or you can use the sietraj command as shown below. Of course, in an actual application you would be averaging over more than just 3 snapshots. It is useful to plot a time series of the calculated snapshot ΔG values to make sure that you are averaging over a stable part of the trajectory, i.e., with little drift in the mean value over time.
prompt% sietraj -ave /tmp/bcl2.out Energies in kcal/mol: Average StdErr Stdev Inter vdW -28.83 0.78 1.36 Inter Coulomb -2.01 0.11 0.20 Reaction Field 9.04 0.37 0.64 Cavity -5.49 0.11 0.20 Constant -2.89 0.00 0.00 ----------------------- Delta G -5.75 0.08 0.14 Sample size= 3 Coefficients used: alpha= 0.104758 gamma= 0.012894 const= -2.89 Delta_G = alpha * (vdw + Coul + RF + Cav) + constant Cav = gamma * Delta_SA
The solvated interaction energy (SIE) is calculated for each snapshot via rigid infinite separation of the target and ligand. SIE is the sum of the intermolecular van der Waals and coulomb interactions plus the change in reaction field energy and nonpolar solvation energy. In this approximation, the internal energies of the target and ligand cancel out and do not contribute. The SIE is then scaled by an empirically determined factor, α, obtained by fitting to a training set of 99 protein-ligand complexes. The scaling can be considered a crude treatment of entropy-enthalpy compensation (Nam et al., 2007; Chen et al., 2004).

Sample run 2
The second example is a protein-protein complex with 3844 atoms. Each protein has about 1900 atoms. 1VET.trj is a trajectory file with 11 snapshots with all water molecules removed. 1VET.prmtop is the corresponding prmtop file for the complex.
prompt% cd $BRIMM_ROOT/apps/sietraj/examples prompt% $BRIMM_ROOT/bin/sietraj -pt 1VET.prmtop -trj 1VET.trj -sf 1 -ef 11 -inc 4 \ -tr 1-1947 -lr 1948-3844 -o /tmp/1VET.out -sie Loaded Britcl-v2.3c-64.so Britcl_functions.tcl 338 2011-01-12 1 -21.931 97.425 -1119.768 -601.006 -616.187 -151.209 -99.597 -2201.314 5 -22.446 106.714 -1129.132 -619.217 -616.629 -158.472 -107.433 -2131.850 9 -23.452 103.110 -1166.878 -635.682 -634.305 -161.664 -110.797 -2088.921 prompt% $BRIMM_ROOT/bin/sietraj -ave /tmp/1VET.out Energies in kcal/mol: Average StdErr Stdev Inter vdW -157.12 2.53 4.37 Inter Coulomb -105.94 2.71 4.69 Reaction Field 102.42 2.21 3.82 Cavity -27.60 0.34 0.60 Constant -2.89 0.00 0.00 ----------------------- Delta G -22.61 0.36 0.63 Sample size= 3 Coefficients used: alpha= 0.104758 gamma= 0.012894 const= -2.89 Delta_G = alpha * (vdw + Coul + RF + Cav) + constant Cav = gamma * Delta_SA

Virtual alanine mutation
Virtual alanine mutations can be performed by adding "-mut2ala resnum" to the command line, where resnum is the residue number of the amino acid to be mutated. The mutation is carried out on the fly with partial charges and atom types assigned according to the FF99SB parameters. Note that disulfide bonded CYS cannot be mutated. Also, at the moment GLY cannot be mutated to ALA.
prompt% cd $BRIMM_ROOT/apps/sietraj/examples prompt% $BRIMM_ROOT/bin/sietraj -pt bcl2.prmtop -trj bcl2.trj -sf 1 -ef 11 -inc 4 \ -tr 1-2611 -lr 2612-2654 -o /tmp/bcl2.out -sie -mut2ala 94 Loaded Britcl-v2.3c-64.so Britcl_functions.tcl 338 2011-01-12 Res # 94 (LEU) will be mutated to Ala 1 -5.580 10.277 -1274.816 -1279.641 -5.453 -28.192 -2.075 -440.886 5 -5.497 8.231 -1286.386 -1289.585 -5.031 -25.911 -1.777 -421.056 9 -5.802 9.330 -1276.603 -1280.910 -5.023 -29.161 -2.210 -446.509

Notes
  • The prmtop file format is that of Amber 7 or later.
  • If you wish to run the calculation on a single-snapshot restrt/inpcrd file, convert the file to a 10-column trajectory file using the Amber's ptraj utility. You can also use the crd2trj.sh and crd2trj.awk scripts found in the apps/sietraj subdirectory.
  • The trajectory file must be generated with the nobox option in ptraj. Otherwise, the snapshot frames will not be read in properly.
  • sietraj calls a Tcl script to carry out the SIE calculation. For more details about that script and the internal settings used, click here.
How to cite
If you use sietraj in your work, please cite the following papers:
  • Cui, Q.; Sulea, T.; Schrag, J. D.; Munger, C.; Hung,M.-N.; Nam, M., Cygler, M.; Purisima, E. O. Molecular Dynamics and Solvated Interaction Energy Studies of Protein-Protein Interactions: the MP1-p14 Scaffolding Complex. J. Mol. Biol. 2008, 379, 787-802. [link]
  • 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. [link]
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. Chen, W.; Chang, C.E.; Gilson, M. K.. Calculation of Cyclodextrin Binding Affinities: Energy, Entropy, and Implications for Drug Design. Biophys.J. 2004, 87, 3035-3049.
  4. 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.
  5. Purisima, E. O. Fast Summation Boundary Element Method for Calculating Solvation Free Energies of Macromolecules. J. Comput. Chem. 1998, 19, 1494-1504.
  6. Purisima, E. O.; Nilar, S. H. A simple yet accurate boundary element method for continuum dielectric calculations. J. Comput. Chem. 1995, 16, 681-689.
Recent publications that use sietraj
  1. Lill, M., and Danielson, M. (2011) Computer-aided drug design platform using PyMOL, J. Comput.-Aided Mol. Des. 25, 13-19.
  2. Yang, B., Hamza, A., Chen, G., Wang, Y., and Zhan, C.-G. (2010) Computational Determination of Binding Structures and Free Energies of Phosphodiesterase-2 with Benzo[1,4]diazepin-2-one Derivatives, J. Phys. Chem. B 114, 16020-16028.
  3. Rodriguez-Granillo, A., Crespo, A., and Wittung-Stafshede, P. (2010) Interdomain Interactions Modulate Collective Dynamics of the Metal-Binding Domains in the Wilson Disease Protein, J. Phys. Chem. B 114, 1836-1848.
  4. Rodriguez-Granillo, A., Crespo, A., Estrin, D. A., and Wittung-Stafshede, P. (2010) Copper-Transfer Mechanism from the Human Chaperone Atox1 to a Metal-Binding Domain of Wilson Disease Protein, J. Phys. Chem. B 114, 3698-3706.
  5. Okamoto, M., Takayama, K., Shimizu, T., Muroya, A., and Furuya, T. (2010) Structure-activity relationship of novel DAPK inhibitors identified by structure-based virtual screening, Bioorg. Med. Chem. 18, 2728-2734.
  6. Wimmerov, M., Mishra, N., Pokorn, M., and Koca, J. (2009) Importance of oligomerisation on Pseudomonas aeruginosaLectin-II binding affinity. In silico and in vitro mutagenesis, J. Mol. Model. 15, 673-679.
  7. Wang, Y. T., Su, Z. Y., Hsieh, C. H., and Chen, C. L. (2009) Predictions of Binding for Dopamine D2 Receptor Antagonists by the SIE Method, J. Chem. Inf. Model. 49, 2369-2375.
  8. Rangarajan, E. S., Proteau, A., Cui, Q., Logan, S. M., Potetinova, Z., Whitfield, D., Purisima, E. O., Cygler, M., Matte, A., Sulea, T., and Schoenhofen, I. C. (2009) Structural and Functional Analysis of Campylobacter jejuni PseG: A UDP-Sugar Hydrolase from the Pseudaminic Acid Biosynthetic Pathway, J. Biol. Chem. 284, 20989-21000.
  9. Liu, Z., Zheng, X., Yang, X., Wang, E., and Wang, J. (2009) Affinity and specificity of levamlodipine-human serum albumin interactions: insights into its carrier function, Biophys. J. 96, 3917-3925.
  10. Rodriguez-Granillo, A., Sedlak, E., and Wittung-Stafshede, P. (2008) Stability and ATP Binding of the Nucleotide-binding Domain of the Wilson Disease Protein: Effect of the Common H1069Q Mutation, J. Mol. Biol. 383, 1097-1111.
  11. Cui, Q., Sulea, T., Schrag, J. D., Munger, C., Hung, M.-N., Nam, M., Cygler, M., and Purisima, E. O. (2008) Molecular Dynamics - Solvated Interaction Energy Studies of Protein-Protein Interactions: the MP1-p14 Scaffolding Complex, J. Mol. Biol. 379, 787-802.
How to contact us
Send feedback, suggestions or queries to ccb@bri.nrc.ca.