- Contemporary messages sorted: [ by date ] [ by thread ] [ by subject ] [ by author ] [ by messages with attachments ]

From: Marek Maly <marek.maly.ujep.cz>

Date: Mon, 07 May 2012 03:26:46 +0200

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

Date: Mon, 07 May 2012 03:26:46 +0200

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/amberReceived on Sun May 06 2012 - 19:00:04 PDT

Custom Search