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 SIE and MM-PBSA is the same. Hence, all the caveats of implicit solvation and neglect of vibrational entropy 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.
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_ROT/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.1e.so
Britcl_functions.tcl,v 1.9 2007/11/28
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 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 (Naïm 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.1e.so
Britcl_functions.tcl,v 1.9 2007/11/28
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
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:
- Naïm, 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]
- Cui, Q.; Sulea, T.; Schrag, J. D.; Munger, C.; Hung,M.-N.; Naïm, 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]
References
- Bhat, S.; Purisima, E. O. Molecular surface generation using a variable-radius solvent probe. Proteins: Structure, Function, and Bioinformatics 2006, 62, 244-261.
- Chan, S. L.; Purisima, E. O. Molecular Surface Generation Using Marching Tetrahedra. J. Comput. Chem. 1998, 19, 1268-1277.
- 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.
- Naïm, 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.
- Purisima, E. O. Fast Summation Boundary Element Method for Calculating Solvation Free Energies of Macromolecules. J. Comput. Chem. 1998, 19, 1494-1504.
- Purisima, E. O.; Nilar, S. H. A simple yet accurate boundary element method for continuum dielectric calculations. J. Comput. Chem. 1995, 16, 681-689.
How to contact us
Send feedback, suggestions or queries to ccb@bri.nrc.ca.
|
