On Thu, 2014-01-02 at 18:48 +0800, Nitin Sharma wrote:
> Dear amber users,
>
> I used following protocol for MD simulation for several ligands against same protein
>
> run_1minimize
>
> restrained energy minimization to minimize water molecules and counterions only
>
> run_2minimize
>
> energy minimization to minimize the entire system
>
> run_3md
>
> heat up system with protein restraints 20ps
>
> run_4md
>
> equilibrate density with protein restraints 1000 ps
>
> run_5md
>
> production run 30ns
>
>
> I used the production run trajectory, receptor, ligand, complex prmtop files, complex-solvated pmtop files to calculate binding free energy using MM-GBSA and entropy was calculated by Quasi-Harmonic entropy calculation.
Any solute entropy approximation is problematic, including both the
quasi-harmonic and normal mode approximations. Because you included the
entropic contributions in this way, it does not surprise me that your
total binding free energy came out positive. There have been some
attempts (some more successful than others) at using end-state free
energy methods (like MM/PBSA) to compute absolute binding free energies,
but they will require approaches to computing solute entropy not
currently implemented in any of the MM/PBSA scripts bundled with Amber.
MM/PBSA analyses are typically poor predictors of absolute binding
affinities, and are better used as a relative metric for rank-ordering
similar drug-target affinities (similar ligands with the same receptor,
or the same ligand with various mutants of the same receptor, for
example). In your case, the relative binding free energy (7.8 vs. 6.6
kcal/mol) suggests that the second system has a stronger binding
affinity. This trend reverses when the quasi-harmonic entropy is not
included, however, so I would not trust MM/GBSA to correctly distinguish
the better binder in this example (since quasi-harmonic entropies
typically require very long simulations to converge).
> However, I am getting positive value for Delta G
>
[snip]
>
> Differences (Complex - Receptor - Ligand):
> Energy Component Average Std. Dev. Std. Err. of Mean
> -------------------------------------------------------------------------------
> VDWAALS -43.9240 3.6080 0.0659
> EEL -5.7609 3.3182 0.0606
> EGB 21.8073 3.2526 0.0594
> ESURF -5.2260 0.3916 0.0071
>
> DELTA G gas -49.6849 5.2006 0.0949
> DELTA G solv 16.5813 3.1147 0.0569
>
> DELTA TOTAL -33.1035 3.4451 0.0629
>
>
> -------------------------------------------------------------------------------
> -------------------------------------------------------------------------------
> Using Quasi-harmonic Entropy Approximation: DELTA G binding = 7.8414
> -------------------------------------------------------------------------------
> -------------------------------------------------------------------------------
>
>
> Using different ligand against protein and follwing same protocol I got WARNING: INCONSISTENCIES EXIST WITHIN INTERNAL POTENTIAL TERMS. THE VALIDITY OF THESE RESULTS ARE HIGHLY QUESTIONABLE.
> In addition to VDWAALS, EEL, EGB and ESURF I got values for BOND , ANGLE, DIHED, 1-4 VDW, 1-4 EEL too. The Delta G was positive in this case also.
>
> | Run on Thu Jan 2 07:11:57 2014
> |
> |Input file:
> |--------------------------------------------------------------
> |Input file for running PB and GB
> |&general
> | endframe=3000, verbose=1,netcdf=1,
> | entropy=1,
> |/
> |&gb
> | igb=2, qm_theory = 'PM3', saltcon=0.100
> |/--------------------------------------------------------------
> |MMPBSA.py Version=13.0
> |Solvated complex topology file: complex_solvated.prmtop
> |Complex topology file: complex.prmtop
> |Receptor topology file: receptor.prmtop
> |Ligand topology file: ligand.prmtop
> |Initial mdcrd(s): 5md.mdcrd
> |
> |Receptor mask: ":1-154"
> |Ligand mask: ":155"
> |Ligand residue name is "UNK"
> |
> |Calculations performed using 3000 complex frames.
> |
> |Generalized Born ESURF calculated using 'LCPO' surface areas
> |
> |All units are reported in kcal/mole.
> |All entropy results have units kcal/mol (Temperature is 298.15 K).
> -------------------------------------------------------------------------------
> -------------------------------------------------------------------------------
> ENTROPY RESULTS (QUASI-HARMONIC APPROXIMATION) CALCULATED WITH PTRAJ:
> Translational Rotational Vibrational Total
> Complex: 16.4531 16.6222 2422.9884 2456.0637
> Receptor: 16.4251 16.5858 2373.0474 2406.0583
> Ligand: 13.3714 11.4856 63.6693 88.5264
>
> DELTA S: -13.3434 -11.4493 -13.7283 -38.5210
> -------------------------------------------------------------------------------
> -------------------------------------------------------------------------------
>
> GENERALIZED BORN:
>
> WARNING: INCONSISTENCIES EXIST WITHIN INTERNAL POTENTIAL
> TERMS. THE VALIDITY OF THESE RESULTS ARE HIGHLY QUESTIONABLE
> Complex:
> Energy Component Average Std. Dev. Std. Err. of Mean
> -------------------------------------------------------------------------------
> BOND 502.3616 18.8777 0.3447
> ANGLE 1372.5913 30.2290 0.5519
> DIHED 2775.0726 23.1393 0.4225
> VDWAALS -1234.5779 22.3262 0.4076
> EEL -10728.2002 165.9006 3.0289
> 1-4 VDW 599.5937 10.5785 0.1931
> 1-4 EEL 4764.7011 37.4519 0.6838
> EGB -2120.6543 164.5357 3.0040
> ESURF 64.6350 2.1923 0.0400
>
> G gas -1948.4579 180.1894 3.2898
> G solv -2056.0194 163.0035 2.9760
>
> TOTAL -4004.4772 39.7967 0.7266
>
>
> Receptor:
> Energy Component Average Std. Dev. Std. Err. of Mean
> -------------------------------------------------------------------------------
> BOND 486.3270 18.6276 0.3401
> ANGLE 1338.2341 29.7430 0.5430
> DIHED 2734.3702 23.0076 0.4201
> VDWAALS -1182.3317 21.4647 0.3919
> EEL -10612.6402 163.2855 2.9812
> 1-4 VDW 573.6460 10.4485 0.1908
> 1-4 EEL 4752.3195 37.9485 0.6928
> EGB -2110.2623 164.6202 3.0055
> ESURF 65.0172 2.0815 0.0380
>
> G gas -1910.0751 178.6794 3.2622
> G solv -2045.2451 163.1962 2.9795
>
> TOTAL -3955.3202 38.4262 0.7016
>
>
> Ligand:
> Energy Component Average Std. Dev. Std. Err. of Mean
> -------------------------------------------------------------------------------
> BOND 16.0346 3.3559 0.0613
> ANGLE 34.3571 4.3665 0.0797
> DIHED 40.7029 3.5357 0.0646
> VDWAALS -9.0212 1.5387 0.0281
> EEL -101.6651 3.0428 0.0556
> 1-4 VDW 25.9476 1.7246 0.0315
> 1-4 EEL 12.3816 5.4734 0.0999
> EGB -40.6722 3.8112 0.0696
> ESURF 4.7384 0.0582 0.0011
>
> G gas 18.7376 6.8351 0.1248
> G solv -35.9338 3.7993 0.0694
>
> TOTAL -17.1961 6.1897 0.1130
>
>
> Differences (Complex - Receptor - Ligand):
> Energy Component Average Std. Dev. Std. Err. of Mean
> -------------------------------------------------------------------------------
> BOND -0.0000 0.0001 0.0000
> ANGLE -0.0000 0.0001 0.0000
> DIHED -0.0005 0.0079 0.0001
> VDWAALS -43.2251 4.3363 0.0792
> EEL -13.8948 8.8425 0.1614
> 1-4 VDW -0.0000 0.0000 0.0000
> 1-4 EEL -0.0000 0.0001 0.0000
By default, MMPBSA.py makes sure that bonded interactions --
specifically the BOND, ANGLE, DIHED, and 1-4 nonbonded interactions --
completely cancel out in the total binding free energy when you use a
single trajectory approach. This should be rigorously correct, since
the protein and ligand conformations are identical in the bound and
unbound states, and MMPBSA.py is valid only for noncovalent binding.
Making sure these contributions to the binding free energy are 0 is a
useful error check. Therefore, when MMPBSA.py verifies these
contributions are 0, it removes them from the results (unless you set
change that behavior with the 'verbose' variable as described in the
manual).
I suspect what you are hitting here is a minor issue in tleap when it
deals with improper torsions with 2 atoms of the same type. Overall,
the DIHED contribution to the binding free energy is negligible (but
still has at least one frame larger than the 0.001 kcal/mol cutoff that
MMPBSA.py uses to issue a warning). You should be safe to ignore this
warning in this case.
> I read in one forum that such error can be due to small production run. Is 30ns not good enough to calculate binding free energy
In general, the problem is using any end-state free energy method for
calculating a binding free energy. The question "is 30 ns long enough"
depends strongly on the system you are simulating. Sometimes the answer
is 'yes', other times it is 'no'. Typically, larger systems require
longer simulations to converge observable properties. Checking that
various properties (radius of gyration, native contact populations, RMSD
distributions, comparisons of average structures, etc.) have converged
on the time scale of your simulation would be helpful to answer this
question. Energies are actually rather poor indicators of convergence,
in my experience (since biomolecules often have numerous states that are
nearly degenerate).
HTH,
Jason
--
Jason M. Swails
BioMaPS,
Rutgers University
Postdoctoral Researcher
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu Jan 02 2014 - 05:30:03 PST
