Respected sir,
>>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.
Is there anything I can check to negate this error? Or I can ignore this error?
>> 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.
Which value are you referring to? Also what should be safe zone to ignore such error.
Thanks and regards,
Nitin Sharma
-----Original Message-----
From: Jason Swails [mailto:jason.swails.gmail.com]
Sent: Thursday, January 02, 2014 9:24 PM
To: amber.ambermd.org
Subject: Re: [AMBER] binding free energy analysis
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
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Wed Jan 08 2014 - 09:00:02 PST
