Electrostatic binding energy

For the impatient

In the following example, we fetch a complex structure from the web (1AVZ, or use a local file), assign the first and second chain as receptor (-r 0 1) and the third chain as ligand (-l 2) and then calculate the electrostatic free energy of binding between the two:

~> delphiBinding.py -i 1AVZ -r 0 1 -l 2

This will (after some minutes at least) generate a file dg_delphi.out with the delta G and some components and it will create some pickled python objects containing the complex with additional data as a Biskit.Dock.Complex object. Add the -v option for additional output along the way and call the script without any arguments to get a long help screen with all the options and parameters.

Now read on for some more details.

Prepare complex structures

You can skip this step and directly start the delphi calculation (delphiBinding.py) from a PDB file. However, it's a good idea to create and inspect the starting complex in a separate step first. So, let's assume we have a PDB file with a protein-protein complex. An example would be the Barnase/Barstar structure 1BGS. This structure contains several chains (crystal copies) and we have to decide which ones we are going to use as a receptor and ligand structure. Inspect the structure with your favorite viewer. You can also issue some quick biskit commands to get an overview. For example, load the structure into the biskit interactive python shell with bis.py:

~> bis.py 1BGS
bispy [pickled_object or PDB_file]

Opens an interactive python prompt with all Biskit/* modules and Biskit.tools
functions imported into the default namespace. Tab-completion is enabled.

A single argument can point either to a pickled python object which is then
unpickled and put into a variable x, or can point to a PDB file which is then
loaded into a PDBModel object x.

Loaded 1BGS and put it into variable x.
>>>
>>> x.report()
[PDBModel 1BGS  4987 atoms,  796 residues, 12 chains]
      * chain 0 (A): AQVINTFDGVADYLQTYHKLPDNYITKSEAQALGWVASKGNLADVAPGKSIGGDIFSN..
      * chain 1 (B): AQVINTFDGVADYLQTYHKLPDNYITKSEAQALGWVASKGNLADVAPGKSIGGDIFSN..
      * chain 2 (C): AQVINTFDGVADYLQTYHKLPDNYITKSEAQALGWVASKGNLADVAPGKSIGGDIFSN..
      * chain 3 (E): KKAVINGEQIRSISDLHQTLKKELALPEYYGENLDALWDALTGWVEYPLVLEWRQFEQ..
      * chain 4 (F): KKAVINGEQIRSISDLHQTLKKELALPEYYGENLDALWDALTGWVEYPLVLEWRQFEQ..
      * chain 5 (G): KKAVINGEQIRSISDLHQTLKKELALPEYYGENLDALWDALTGWVEYPLVLEWRQFEQ..
      * chain 6 (A): ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      * chain 7 (B): ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      * chain 8 (C): ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      * chain 9 (E): ~~~~~~~~~~~~~~~~~~~~~
      * chain 10(F): ~~~~~~~~~~~~~~~~~~~~~~~~
      * chain 11(G): ~~~~~~~~~~~~~~~~~~
 source: '1BGS'
 15 atom profiles:    ['name', 'residue_number', 'insertion_code', 'alternate..
  1 residue profiles: ['biomol']
  4 info records:     ['date', 'pdb_code', 'resolution', 'compound']

So what we can see is that there are 4 copies of two different protein chains. Chain 0 (A) contains the first copy of the receptor and chain 3 (E) the first copy of the ligand. The protein chains are followed by additional "chains" of solvent associated with each copy. So 0 and 3 are the chains we want. You could use x.takeChains( [0,3] ) to isolate those chains but for now we stick to ready-made scripts rather than programming. Close bis.py with CTRL-D.

Create a Biskit.Complex object from the PDB file and save it to the file input.complex:

~> pdb2complex.py -c 1BGS -r 0 -l 3 -o input.complex
Loading...
Removing water...
Extracting rec and lig...
Saving...
~>

Note that we have here fetched the PDB from the web. You can of course also supply a PDB file or (which is by far the fastest option) a pickled Biskit.PDBModel object as an input to -c. Call pdb2complex.py without any arguments to see an help screen with the available options. For now just note that we now have a file input.complex (as well as rec.model and lig.model with the separated structures).

Perform DelPhi calculations

The actual calculation is performed with delphiBinding.py (in scripts/Dock). Call delphiBinding.py without any arguments to see an help screen with the available options. For now, let's be brave and call this script with all options left at default values:

~> delphiBinding.py -c input.complex -v

The -v switches into verbose mode. A good idea, because your patience is now going to be tested. Depending on the size of the complex and the hardware you are sitting on the following 6 DelPhi runs can easily take 20 min or much longer. It's a good idea to first make a much quicker test run with a sparse grid. This can be done by setting the scale option to, for example, 0.8 (0.8 grid points per Angstroem):

~> delphiBinding.py -c input.complex -v -scale 0.8

Re-building hydrogen atoms...
executing: /usr/local/bin/reduce  -BUILD -Nterm2 -DB
/data/raik/py/biskit/Biskit/data/reduce/reduce_wwPDB_het_dict.txt /tmp/reduce_1kqcP5_in.pdb

... [ many more detailed lines and then, finally: ] ...

Final Result
============
 Rec :  Lig  dG(kT)    solv    coul    ions dG(kcal/mol)
1BGS : 1BGS   12.34  234.68 -227.32    4.97    7.32

energy values written to  /home/raik/dg_delphi.out

Here we go, the electrostatic free energy of binding is calculated to be 7.32 kcal/mol -- that means this component is supposed to be unfavorable. There is a large favorable electrostatic interaction (coulomb) between charges but it is more than cancelled out by the cost of desolvation. The Barnase/Barstar complex has a highly charged interface. All these charges can interact with water and ions and the burial of these charges comes at a significant cost. In reality of course, this will be compensated by solvent entropy gains (hydrophobic effect) and (to a lesser degree) vdW interactions but DelPhi knows nothing about these. In any case, -scale 0.8 is really a brutal approximation. Run the calculation once more with a proper scaling of -scale 2.3 (default) or similar. This will change the result substantially (but take much longer).

Inspect results

The last two lines of the output is also what you will find in the file dg_delphi.out. Another output of the script is a new version of the Biskit.Complex (delphi.complex) which now has the detailed Delphi results stored within. Python programmers can load (unpickle) this object and retrieve the dictionaries with all the data from the complex. For interactive work, use again bis.py:

~> bis.py delphi.complex

Loaded delphi.complex and put it into variable x.
>>> com = x
>>> com.info
{'date': '2011/11/01:23.28.36', 'dG_delphi': {'ions': 4.9700000000002547, 'dG_kt':
12.337200000001076, 'dG_kcal': 7.3159999999999998, 'solv': 234.68419999999992, 'dG_kJ': 30.584,
'coul': -227.3169999999991}, 'delphi_0salt': {'egrid': 6878.0029999999997, 'erxn': -
1048.0840000000001, 'scharge': 1.0465, 'ecoul': -22155.66, 'erxnt': -46044.0, 'eself': -
44995.910000000003}}

>>> print 'dG electrostatic in kJ/mol: ', com['dG_delphi']['dG_kJ']
dG electrostatic in kJ/mol:  30.584

The result dictionary is described in the documentation of the Biskit.Dock.DelphiBindingEnergy class. So here is what it contains:

{'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 'raw' delphi results are associated to receptor and ligand models within the Complex. Ligand and receptor contain two delphi result dictionaries within their info record -- one with the results from the zero-salt calculation and one with the result from the delphi run with salt:

>>> rec = com.rec()
>>> lig = com.lig()
>>> rec.info['delphi_0salt']
{'egrid': 3633.1379999999999, 'erxn': -569.98929999999996, 'scharge': -0.4128, 'ecoul': -12190.84,
'erxnt': -25220.169999999998, 'eself': -24650.18}
>>> rec.info['delphi_0.15salt']
{'egrid': 3631.0529999999999, 'erxn': -570.53099999999995, 'scharge': -0.4128, 'ecoul': -12190.84,
'erxnt': -25220.709999999999, 'eself': -24650.18}

That means the reaction field energy at 0.15M salt can be retrieved from the complex like this:

>>> com.rec()['delphi_0.15salt']['erxn']
-570.53099999999995

One last thing to note is that the structure of the delphi.complex object is not any longer identical to the input complex. The reason is the addition of hydrogens and the correction (flipping) of charged some moieties with the reduce program. You can generate a new PDB model with the complex used for the delphi calculation directly from the python prompt:

>> com.writeComplex('delphi.pdb')

A customized workflow

The delphiBinding.py script is a quite minimalist wrapper around the DelphiBindingEnergy class. Depending on the input structure you may need some extra pre-cautions to clean up your PDB, customize the capping of chain ends. In that case it is better to build your own pipeline from the components used by DelphiBindingEnergy. Have a look at the example workflow in

• Custom delphi workflow

I used this script to calculate electrostatic free energies for a dimer of GFP (Green Fluorescent Protein) structures and to explore whether there is any electrostatic steering at different distances of receptor and ligand.

The script first cleans up the input structure with PDBCleaner, then uses the Reduce wrapper and the AtomCharger explicitly to assign hydrogens and partial atomic charges. It only then creates a Biskit.Dock.Complex object from the two chains of the modified PDBModel and creates a gap between the receptor and ligand. Only then the calculation is delegated to DelphiBindingEnergy and the default calls to reduce and AtomCharger are switched off (because that has been done explicitly before).

There are some hard-coded numbers in this script. So it can not be used without modifications. Namely which chains to use as receptor and ligand (0, and 1) is hard-coded. Another issue was the capping of chain breaks. The GFP input structures have unresolved N-termini and the removal of the (non-standard) chromophore residues creates a short gap in each protein chain. Because of the tight gap, the capping residue actually often clashed into other parts of the structure which messed up the Delphi calculation. The workaround was to first add the capping residues to each N-terminus, let reduce do its job adding hydrogens and map charges as before. But then remove the ACE residues before the actual delphi calculation. That prevents any terminal charge without introducing clashes.

Download the script and run line by line interactively to get an idea of what's going on. Please also check out the (ample) documentation in the python code of the Biskit classes used to better understand what's going on. Of course this is just an example and you may have completely other needs.

Good luck!

Dependencies

Two external programs are required for this pipeline:

• DelPhi -- see setup Delphi
• reduce -- see setup Reduce

Before reporting any problems, please test your local Biskit and delphi installation with the built-in unittest case:

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