Hello all,
after the first successful steps with TI in Amber (tut. A9 and solvation
energy of small molecules) I
would like to solve the last problem "TI of constrained ligand" to be able
to calculate
the whole free energy of the receptor-ligand binding in solvent using
Amber.
Here I would like to follow approach described in this article:
--------------------------------------------------------------------------------------
Boresch, S.; Tettinger, F.; Leitgeb, M.; Karplus, M. Absolute binding free
energies:
A quantitative approach for their calculation. J. Phys. Chem. B 2003, 107,
9535-9551.
--------------------------------------------------------------------------------------
Here the idea of the transfer of the ligand (L) from the binding site on
receptor (R) in solvent (s) into vacuum (v) is expressed
like this:
dA* dAr
(RL)s -----> (R)s***(L)v----->(R)s + (L)v (1)
where we have three states:
state 1
(RL)s = receptor-ligand complex in solvent
state 2
(R)s***(L)v = artificial intermediate state where ligand is decoupled from
it's original environment (R and s) but
is restrained in it's native bound state (position) respect
to receptor.
state 3
(R)s + (L)v = final state when ligand is decoupled from R and s and is
free to move (free R in solvent + free L in vacuum).
and
dA* is change in free energy of the system while going from the state 1 to
state 2 (calculated using TI MD runs)
dAr is change in free energy of the system while going from the state 2 to
state 3 (possible to calculate analytically - formula (14) or (32) in
article)
the total change in free energy of the system connected with decoupling of
the ligand from the binding site of the receptor in solvent is then dA =
dA* + dAr.
My actual question is about practical realization of the first step ( 1
---> 2 ) or more precisely about the calculation of the dA* in Amber.
I suppose that the restraints (distance/angle) are zero in state 1 and
are gradually turned on during the first step (1--->2) with maximum/final
values in state 2
(these values (force constants) are then used in analyt. calc. of dAr) as
the whole process (1--->2--->3) should connect two real states. I hope
this assumption is OK.
So my question is how to do this first step with Amber assuming "One step"
approach (using simultaneously softcore potentials for vdw and el.
interactions).
I suppose that the "V0 system" is solvated complex and the "V1 system" is
just solvated receptor
(of course with all common atoms (R,s) in the same positions/order in both
systems) with the ligand scmask in V0. Am I right ?
Restraints for the ligand are defined (in proper RST file) just in case of
the "V0 system" as the ligand is not present in
the "V1 system". Am I right ?
I am not sure here how to set "dvdl_norest" parameter although in the
context of such calculations is mentioned non-zero value (Amber manual,
mailing list)
to ensure that reported dv/dl values will not contain restraint dv/dl
contribution but on the other hand in case of calculation of dA* from (1)
also
restraint dv/dl contribution should be taken in account in my opinion as
dA* represents the TOTAL free energy change while going from the state 1
to state 2 which
means dvdl_norest=0.
So dvdl_norest=0 or dvdl_norest=1 in this particular case ?
I already did my first trial using above mentioned setup on well known
testing case "T4 Lysozyme+benzen" for lambdas
0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9. using
restraint force constants of 20 kcal/(mol Å2) or kcal/(mol rad2) in case
of angle rest.
I did not define any weights so I included just minimal &wt record: &wt
type='END' / and defined all six required restraints in external RST file.
Here is one of such record from my RST file:
&rst
ixpk= 0, nxpk= 0, iat= 1394, 2604, r1= 0.1, r2= 4.1, r3= 4.1, r4= 8.1,
rk2=20.0, rk3=20.0,
&end
All six restraints were successfully read in and in different lambdas
multiplied with different number:
#lambda 0.1
Scaling up weight of 20.000 by 1.111
New weight 22.222
#lambda 0.2
Scaling up weight of 20.000 by 1.250
New weight 25.000
#lambda 0.3
Scaling up weight of 20.000 by 1.429
New weight 28.571
#lambda 0.4
Scaling up weight of 20.000 by 1.667
New weight 33.333
#lambda 0.5
Scaling up weight of 20.000 by 2.000
New weight 40.000
#lambda 0.6
Scaling up weight of 20.000 by 2.500
New weight 50.000
#lambda 0.7
Scaling up weight of 20.000 by 3.333
New weight 66.667
#lambda 0.8
Scaling up weight of 20.000 by 5.000
New weight 100.000
#lambda 0.9
Scaling up weight of 20.000 by 10.000
New weight 200.000
How the above reported scaling constants (e.g. 1.111, 1.250, 1.429, 1.667
...) are calculated/connected with the actual value of lambda ?
This info would be helpful to be able to set some proper initial value of
the force restraint constants (in RST file) to obtain their desired target
values (for the biggest lambda). Of course there is also possibility to do
very short pre-run (small nstlim) for the concrete system/lambda set and
read the scaling constants (respectively the biggest one) from the output
files of this short pre-run.
Thanks a lot in advance for any comment !
Best wishes,
Marek
--
Tato zpráva byla vytvořena převratným poštovním klientem Opery:
http://www.opera.com/mail/
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Sun May 06 2012 - 19:00:04 PDT
