Personal tools
You are here: Home / Documentation / Interfaces with external programs / Delphi -- binding energy

Delphi -- binding energy

How to calculate the electrostatic component of the free energy of binding (protein-protein or Dna-protein) with the Biskit Delphi wrappers.

The Biskit.Dock.DelphiBindingEnergy class automates the calculation of the electrostatic component of the free energy of binding using several rounds of Delphi calculations. The class is based on the Biskit.Delphi wrapper. For end-users, the class is bundled into the script delphiBinding.py (scripts/Dock) where it can be used from the command line.

DelphiBindingEnergy accepts a binary complex (Biskit.Dock.Complex) as input, performs several house keeping tasks (optional capping of free terminals, h-bond optimization and protonation with the reduce program, assignment of Amber partial atomic charges) and then calls Delphi six times:

  • 1 - 3, one run each for complex, receptor, and ligand with zero ionic strength
  • 4 - 6, one run each for complex, receptor, and ligand with custom salt concentration (default: physiological 0.15 M).

The grid position and dimensions are determined once for the complex and then kept constant between all runs so that receptor, ligand and complex calculations can indeed be compared to each other.

The free energy of binding is calculated as described in the Delphi documentation -- see: delphi2004/examples/binding.html -- according to the energy partioning scheme:

dG = dG(coulomb) + ddG(solvation) + ddG(ions)

Please have a look at the source code of DelphiBindingEnergy.bindingEnergy() for a detailed breakup of these terms.

Use:

>>> com = Biskit.Complex( 'myrec.pdb', 'mylig.pdb' )
>>> dg = DelphiBindingEnergy( com, verbose=True )
>>> r = dg.run()

Variable r will then receive a dictionary with the free energy of binding and its three components:

{'dG_kt': free E. in units of kT,
 'dG_kcal': free E. in units of kcal / mol,
 'dG_kJ': free E. in units of kJ / mol,
 'coul': coulomb contribution in kT ,
 'solv': solvation contribution in kT,
 'ions': salt or ionic contribution in kT }

The modified complex used for the calculation (with hydrogens added and many other small changes) can be recovered from the DelphiBindingEnergy instance:

>>> modified_com = dg.delphicom

The run() method assigns the result dictionary to the info record of both the original and the modified complex. That means:

>>> r1 = com.info['dG_delphi']
>>> r2 = modified_com.info['dG_delphi']
>>> r1 == r2
    True

The result of the original DelPhi calculations (with and without salt for receptor, ligand and complex) are assigned to the info dictionaries of the complex' receptor and ligand model as well as to the complex itself:

com.info['delphi_0.15salt'] ... delphi run with 0.15M salt on complex
com.info['delphi_0salt']    ... delphi run without salt on complex

com.lig().info['delphi_0.15salt'] ...delphi run with 0.15M salt on ligand
com.lig().info['delphi_0salt']    ...delphi run without salt on ligand
com.rec().info['delphi_0.15salt'] ...delphi run with 0.15M salt on receptor
com.rec().info['delphi_0salt']    ...delphi run without salt on receptor

From there the individual values for solvation, ionic and couloumb contributions can be recovered. See Biskit.Delphi for a description of this result dictionary.

Customization:

The most important Delphi parameters (dielectric, salt concentration, grid scales, ion and probe radius) can be adjusted by passing parameters to the constructor (see documentation of __init__). The default parameters should be reasonable. By default, we create a grid that covers every linear dimension to at least 60% (perfil=60) and has a density of 2.3 points per Angstroem (scale=2.3). Such high densities come at much larger computational cost. It is recommended to test different values and average results.

Note:
Any parameters that are not recognized by DelphiBindingEnergy() will be passed on to the Biskit.Delphi instance of each Delphi run and, from there, passed on to Biskit.Executor.

The internal pipeline of DelphiBindingEnergy consists of:

  • adding hydrogens and capping of protein chain breaks and terminals with Biskit.Reduce( autocap=True ). The capping is performed by L{Biskit.PDBCleaner} called from within Reduce.
  • mapping Amber partial atomic charges into the structure with Biskit.AtomCharger()
  • setting up the various delphi runs with Biskit.Delphi.

A point worth looking at is the automatic capping of protein chain breaks and premature terminals with ACE and NME residues. This should generally be a good idea but the automatic discovery of premature C-termini or N-termini is guess work at best. See Biskit.PDBCleaner for a more detailed discussion. You can override the default behaviour by setting autocap=False (no capping at all, this is now the default) and you can then provide a complex structure that is already pre-treated by Biskit.Reduce. For example:

>>> m = PDBModel('mycomplex.pdb')
>>> m = Reduce( m, capN=[0], capC=[2] )
>>> com = Biskit.Complex( m.takeChains([0]), m.takeChains([1,2]) )
>>> dg = DelphiBindingEnergy( com, protonate=False )

In this case, the original structure would receive a ACE N-terminal capping of the first chain and the second chain would receive a C-terminal NME capping residue. The first chain is considered receptor and the second and third chain are considered ligand.

See also:
  • Biskit.PDBCleaner
  • Biskit.Reduce
  • Biskit.AtomCharger
  • Biskit.Delphi

Dependencies:

Two external programs are required for this wrapper:

Test:

Run the built-in unit test to verify your installation:

cd biskit/Biskit/Dock
python delphiBindingEnergy.py

Should (after one or two minutes and lots of other output) yield:

Final result: dG = 12.79 kcal/mol
ok

----------------------------------------------------------------------
Ran 1 test in 144.785s

OK
Help:

The complete list of parameters is documented in the source code (at the moment, the online API documentation is lagging behind, you can generate your own with epydoc though). All methods are documented python-style and can be interrogated at the python prompt:

>>> from Biskit.Dock import *
>>> dG = DelphiBindingEnergy( Complex() )  # create Delphi instance with an empty structure
>>> print dG.__init__.__doc__

Which will yield something like that:

@param com: complex to analyze
@type  com: Biskit.Complex
@param protonate: (re-)build hydrogen atoms with reduce program (True)
                  see L{Biskit.Reduce}
@type  protonate: bool
@param addcharge: Assign partial charges from Amber topologies (True)
@type  addcharge: bool
@param indi: interior dilectric (4.0)
@param exdi: exterior dielectric (80.0)
@param salt: salt conc. in M (0.15)
@param ionrad: ion radius (2)
@param prbrad: probe radius (1.4)
@param bndcon: boundary condition (4, delphi default is 2)
@param scale:  grid spacing (2.3)
@param perfil: grid fill factor in % (for automatic grid, 60)

@param template: delphi command file template [None=use default]
@type  template: str
@param f_radii: alternative delphi atom radii file [None=use default]
@type  f_radii: str
@param topologies: alternative list of residue charge/topology files
                   [default: amber/residues/all*]
@type  topologies: [ str ]
@param f_charges: alternative delphi charge file
                  [default: create custom]
@type  f_charges: str

@param kw: additional key=value parameters for Delphi or Executor:
@type  kw: key=value pairs
::
  debug    - 0|1, keep all temporary files (default: 0)
  verbose  - 0|1, print progress messages to log (log != STDOUT)
  node     - str, host for calculation (None->local) NOT TESTED
                  (default: None)
  nice     - int, nice level (default: 0)
  log      - Biskit.LogFile, program log (None->STOUT) (default: None)