Re: question of MMPBSA calculation combined NAMD and APBS software

From: Aron Broom (broomsday_at_gmail.com)
Date: Thu Oct 11 2012 - 22:18:13 CDT

First, are you trying to get the gas-phase Kd? If so, you don't need the
MMPBSA method, but I'll assume you want the solution Kd, to which some
points can be made.

1) It looks like you are just taking a single snapshot? If this is the
case it will almost certainly not give a reasonable answer. When you're
using the MMPBSA method you are calculating the binding energy as the
difference of VERY large numbers, so you need a well-converged average
value, otherwise your error will be massive in comparison to the number you
are hoping to calculate. You probably want something in the ns range for a
protein ligand system I would think. Or is your script just meant to
calculate from a list of snapshots you are giving NAMD?

2) Related to the first point, if you just calculate the potential energy
with NAMD, and the solvation energy with APBS, you are essentially getting
the enthalpy for the solute from NAMD, and the enthalpy for the solvent
from APBS. This ignores entirely the entropic terms, which are just as
important.

To get the solvent entropy, you need a method for accounting for the
solvent accessible surface area, which is proportional to the solvent
entropy. If you are doing real MMPBSA you will add this term in later
based on your snapshots. If you want a simpler but less accurate method,
you can run your simulation with implicit solvent and "SASA on" and then
those energy terms will automatically be added in, in this case you won't
use APBS and your calculation becomes MMGBSA rather than MMPBSA. Otherwise
you should have explicit waters in your simulation.

To get the solute entropy, you need to run a normal modes or quasi-harmonic
calculation on your series of snapshots. You can find tools for this
within the free AmberTools package (I think the application that handles
this is conveniently named MMPBSA).

If you are trying to get a single accurate Kd, then you absolutely MUST
calculate the entropy for both solvent and solute in addition to your
potential energy terms. If, however, you just want relative Kds between
several different ligands, you MIGHT be able to get away with only
calculating the solvent entropy term (the SASA related one), and ignore the
normal-modes or quasi-harmonic calculation, which I believe is somewhat
computationally intensive. This is just a MIGHT though, it depends on your
ligands with the assumption that their changes in entropy from free to
bound are roughly the same. If they are fairly different in size or
flexibility, this will likely not be true.

Best of luck! Free energy calculations are challenging!

~Aron

On Thu, Oct 11, 2012 at 9:52 PM, baogen duan <dbaogen_at_gmail.com> wrote:

> Dear NAMD users and developers,
>
> I want to calculate the binding energy between protein and small
> molecule using NAMD and APBS software. But the result is weird. The
> absolute value of binding energy is very large which is about thousands of
> kJ/mol. In general, the value is about tens of kJ/mol. whether did I make a
> mistake in the course of calculation. The calculation procedure is as
> follow:
>
> (1) Emm part using NAMD, whether is Emm part in gas phase calculated
> without PBC defined in configure file. And the confiure file is :
>
> *## ADJUSTABLE PARAMETERS*
> *set inputname
> atp
> structure ./$inputname.psf
> coordinates
> ./$inputname.pdb
> set temperature
> 310
> firsttimestep 0 *
> * *
> *## SIMULATION PARAMETERS
> *
> *# Input
> paraTypeCharmm on
> parameters ./par_all27_prot_na.prm
> temperature $temperature*
>
> *# Force-Field Parameters
> exclude scaled1-4
> 1-4scaling 1.0
> cutoff 14.0 *
> *switching on
> switchdist 12.0
> pairlistdist 16.0*
>
> *# Integrator Parameters
> timestep 2.0 ;# 2fs/step
> rigidBonds all ;# needed for 2fs steps
> nonbondedFreq 1
> fullElectFrequency 2
> stepspercycle 1*
>
> *# Constant Temperature Control
> langevin on ;# do langevin dynamics
> langevinDamping 1 ;# damping coefficient (gamma) of 1/ps
> langevinTemp $temperature
> langevinHydrogen off ;# don't couple langevin bath to hydrogens*
> **
> *## EXECUTION SCRIPT*
> *run 0*
> **
> Is the above configure file of Emm calculation in gas phase rignt?
>
> (2) solvation free energy using APBS, the atom charge, radius, epsilon
> are come from CHARMM 22 force field, the follow is a part of parameter file
> as a input for APBS calculation
> *# Format:
> # <RESIDUE> <ATOM> <CHARGE> <RADIUS> <EPSILON>
> # where <RESIDUE> is the residue name, <ATOM> is the atom name, <CHARGE>
> # is the atomic charge in e, <RADIUS> is the van der Waals radius in
> # A, and <EPSILON> is the van der Waals well depth in kJ/mol.*
> *ALA N -0.470000 1.8500 0.8360
> ALA HN 0.310000 0.2245 0.1923
> ALA CA 0.070000 2.2750 0.0836
> ALA HA 0.090000 1.3200 0.0920
> ALA CB -0.270000 2.0600 0.3344
> ALA HB1 0.090000 1.3200 0.0920
> ALA HB2 0.090000 1.3200 0.0920
> ALA HB3 0.090000 1.3200 0.0920
> ALA C 0.510000 2.0000 0.4598
> ALA O -0.510000 1.7000 0.5016*
>
> Would you like to give me some help. Thanks in advance.
>
>
> Best wishes,
> Yours sincerely,
>
> Duan Baogen
>
>

-- 
Aron Broom M.Sc
PhD Student
Department of Chemistry
University of Waterloo

This archive was generated by hypermail 2.1.6 : Tue Dec 31 2013 - 23:22:40 CST