U.S. patent application number 14/138186 was filed with the patent office on 2015-06-25 for methods and systems for calculating free energy differences using a modified bond stretch potential.
This patent application is currently assigned to SCHRODINGER, INC.. The applicant listed for this patent is SCHRODINGER, INC.. Invention is credited to Robert Abel, Lingle Wang.
Application Number | 20150178442 14/138186 |
Document ID | / |
Family ID | 53400322 |
Filed Date | 2015-06-25 |
United States Patent
Application |
20150178442 |
Kind Code |
A1 |
Abel; Robert ; et
al. |
June 25, 2015 |
METHODS AND SYSTEMS FOR CALCULATING FREE ENERGY DIFFERENCES USING A
MODIFIED BOND STRETCH POTENTIAL
Abstract
A method and system for calculating the free energy difference
between a target state and a reference state. The method includes
determining one or more intermediate states using a coupling
parameter, performing molecular simulations to obtain ensembles of
micro-states for each of the system states, and calculating the
free energy difference by an analysis of the ensembles of
micro-states of the system states. The method can be particularly
suited for calculating physical or non-physical transformation of
molecular systems such as ring-opening, ring-closing, and other
transformations involving bond breaking and/or formation. A soft
bond potential dependent on a bond stretching component of the
coupling parameter and different from the conventional harmonic
potential is used in the molecular simulations of the system states
for the bond being broken or formed during the transformation.
Inventors: |
Abel; Robert; (Brooklyn,
NY) ; Wang; Lingle; (New York, NY) |
|
Applicant: |
Name |
City |
State |
Country |
Type |
SCHRODINGER, INC. |
New York |
NY |
US |
|
|
Assignee: |
SCHRODINGER, INC.
New York
NY
|
Family ID: |
53400322 |
Appl. No.: |
14/138186 |
Filed: |
December 23, 2013 |
Current U.S.
Class: |
703/2 ;
703/11 |
Current CPC
Class: |
G16C 20/10 20190201;
G16C 20/30 20190201; G16B 15/00 20190201 |
International
Class: |
G06F 19/12 20060101
G06F019/12 |
Claims
1. A computer-implemented method for computing the free energy
difference between a reference state and a target state, wherein
the reference state and target state each include a common set of
atoms P.sub.AB, and wherein the reference state further includes a
set of atoms P.sub.A, the target state further includes a set of
atoms P.sub.B, the set P.sub.A being present only in the reference
state and not in the target state, and the set P.sub.B being
present only in the target state and not the reference state, where
there exist at least two atoms A.sub.a and A.sub.b, A.sub.a and
A.sub.b being either: (1) not valence-bonded to each other in the
reference state and valence-bonded in the target state, or (2)
valence-bonded to each other in reference state and not
valence-bonded to each other in the target state, the method
comprising: (a) providing a topology, including the bonded
connections between the atoms and the relative spatial arrangements
of the atoms, for all the atoms in P.sub.A, P.sub.B, and P.sub.AB;
(b) determining one or more intermediate system states along a
transformation path between the reference state and the target
state, the transformation path defined by a coupling parameter
.lamda. that modulates the energies arising from inter-atom
interactions for each system state, the coupling parameter .lamda.
including a plurality of components each having a value belonging
to [0,1] and modulating a different type of interaction energy; (c)
performing, using at least one computer processor, molecular
simulations to obtain ensembles of micro-states for the reference
state, the target state, and the intermediate states, wherein
performing molecular simulations for each of the system states
includes calculating a bonded stretch interaction energy between
the atoms A.sub.a and A.sub.b, the bonded stretch interaction
energy being defined by a soft bond potential, wherein the soft
bond potential is a function of a bonded stretch component,
.lamda..sub.sbs, of the coupling parameter .lamda., and does not
include any singular regions for all values of .lamda..sub.sbs
within [0,1] and for all values of the distance r between A.sub.a
and A.sub.b, the soft bond potential further satisfies the
following conditions: when .lamda..sub.sbs is within (0,1), the
soft bond potential is flat when the distance between A.sub.a and
A.sub.b approaches infinity; when A.sub.a and A.sub.b are not
valence bonded in either the reference state or the target state,
the soft bond potential is flat and zero for all distances between
A.sub.a and A.sub.b; and when A.sub.a and A.sub.b are valence
bonded in either the target state or the reference state, the soft
bond potential reverts to a harmonic potential; and wherein the
soft bond potential is wherein the soft bond potential is a
function of (r-r.sub.0).sup.2, where r.sub.0 is the equilibrium
distance between A.sub.a and A.sub.b, and is expressed by: U sbs (
r , .lamda. sbs ) = 1 2 kf ( .lamda. sbs ) ( r - r 0 ) 2 1 1 + g (
.lamda. sbs ) .alpha. ( k , .lamda. sbs ) ( r - r 0 ) 2
##EQU00031## where k is a constant, and the functions f, g and
.alpha. are each continuous functions and satisfy the following
conditions: f(.lamda..sub.sbs=0)=0, f(.lamda..sub.sbs=1)=1,
g(.lamda..sub.sbs=0)=1, g(.lamda..sub.sbs=1)=0,
.alpha.(k,.lamda..sub.sbs<1)>0; and (d) calculating, using at
least one computer processor, the free energy difference between
the reference state and the target state, by way of an analysis of
the ensembles of micro-states obtained at the target state, the
reference state, and the intermediate states.
2. The method of claim 1, wherein the first derivative of the soft
bond potential with respect to .lamda..sub.sbs is continuous and
bounded for all values of .lamda..sub.sbs between 0 and 1 and for
all values of the distance r between A.sub.a and A.sub.b.
3. The method of claim 1, wherein both the first derivative and the
second derivative of the soft bond potential with respect to the
distance r between A.sub.a and A.sub.b are continuous and bounded
for all values of .lamda..sub.sbs, and approach zero when r
approaches infinity.
4. (canceled)
5. (canceled)
6. The method of claim 1, wherein
f(.lamda..sub.sbs)=.lamda..sub.sbs,
g(.lamda..sub.sbs)=1-.lamda..sub.sbs, and
.alpha.(k,.lamda..sub.sbs)=const.
7. The method of claim 1, wherein performing molecular simulations
for each of the system states comprises: if A.sub.a and A.sub.b are
valence-bonded to each other in the reference state and not
valence-bonded in the target state, using a schedule of
.lamda..sub.sbsA and a corresponding soft bond potential for
calculating the bonded stretch interaction energy between A.sub.a
and A.sub.b for each of the intermediate states, wherein
.lamda..sub.sbsA is 1 at the reference state, 0 at the target
state, and varied from 1 to 0 at each intermediate state along the
transformation from the reference state to the target state; if
A.sub.a and A.sub.b are not valence-bonded to each other in the
reference state and valence-bonded in the target state, using a
schedule of .lamda..sub.sbsB and a corresponding soft bond
potential for each of the intermediate states and a soft bond
potential corresponding to the .lamda..sub.sbsB for calculating the
bonded stretch interaction between A.sub.a and A.sub.b, wherein
.lamda..sub.sbsB is 0 at the reference state, 1 at the target
state, and varied from 0 to 1 at each intermediate state along the
transformation from the reference state to the target state.
8. The method of claim 7, wherein performing molecular simulations
for each of the system states further comprises: (a) computing a
bonded angle interaction, using applicable parameters for bonded
angle interactions of a force field, between (i) a bond formed by
A.sub.a and another atom A.sub.c, and (ii) the bond between A.sub.a
and A.sub.b that is being broken or formed by the transformation
from the reference state to the target state; (b) if A.sub.a and
A.sub.b are valence-bonded to each other in the reference state and
not valence-bonded to each other in the target state, multiplying
the computed bonded angle interaction obtained in (a) by a bonded
angle coupling parameter .lamda..sub.baA, wherein .lamda..sub.baA
is 1 at the reference state, 0 at the target state, and varied from
1 to 0 at each intermediate state along the transformation from the
reference state to the target state; and if A.sub.a and A.sub.b are
not valence-bonded to each other in the reference state and
valence-bonded to each other in the target state, multiplying the
computed bonded angle interaction obtained in (a) by a bonded angle
coupling parameter .lamda..sub.baB, wherein .lamda..sub.baB is 0 at
the reference state, 1 at the target state, and varied from 0 to 1
at each intermediate state along the transformation from the
reference state to the target state; and (c) including the bonded
angle interaction obtained in (b) into the total energy of a
simulation step of the corresponding system state,
9. The method of claim 8, wherein performing molecular simulations
for each of the system states further includes: (a) computing a
dihedral angle interaction, using applicable parameters for
dihedral interactions of a force field, of a group of four
connected atoms {A.sub.i, A.sub.j, A.sub.k, A.sub.l}, the group
including both A.sub.a and A.sub.b; (b) if A.sub.a and A.sub.b are
valence-bonded to each other in the reference state and not
valence-bonded to each other in the target state, multiplying the
computed dihedral interaction obtained in (a) by a dihedral angle
coupling parameter .lamda..sub.bdA, wherein .lamda..sub.bdA is 1 at
the reference state, 0 at the target state, and varied from 1 to 0
at each intermediate state along the transformation from the
reference state to the target state; and if A.sub.a and A.sub.b are
not valence-bonded to each other in the reference state and
valence-bonded to each other in the target state, multiplying the
computed dihedral interaction obtained in (a) by a dihedral angle
coupling parameter .lamda..sub.bdB, wherein .lamda..sub.bdB is 0 at
the reference state, 1 at the target state, and varied from 0 to 1
at each intermediate state along the transformation from the
reference state to the target state; and (c) including the dihedral
interaction obtained in (b) into the total energy of the simulation
step of the corresponding system state.
10. The method of claim 9, wherein: if A.sub.a and A.sub.b are
valence-bonded to each other in the reference state and not
valence-bonded to each other in the target state, the bonded angle
interaction and the bonded dihedral interaction coupling parameters
.lamda..sub.baA and .lamda..sub.bdA, are each selected to be 0 when
.lamda..sub.sbsA is smaller than a predefined threshold, and if
A.sub.a and A.sub.b are not valence-bonded to each other in the
reference state and valence-bonded to each other in the target
state, the bonded angle interaction and the bonded dihedral
interaction coupling parameters .lamda..sub.baB and .lamda..sub.bdB
are each selected to be 0 when .lamda..sub.sbsB is smaller than a
predefined threshold.
11. The method of claim 9, wherein performing molecular simulations
for all of the states further includes: (a) computing nonbonded
electrostatic interactions and van der Waals interactions, using
applicable parameters for electrostatic interactions and van der
Waals of a force field, between two atoms A.sub.i and A.sub.j and
the non-bonded exclusion status of the pair (A.sub.i, A.sub.j) is
affected by the transformation from the reference state to the
target state; (b) if A.sub.a and A.sub.b are valence-bonded to each
other in the reference state and not valence-bonded in the target
state, and the nonbonded interactions between A.sub.i and A.sub.j
are excluded in the reference state but not excluded in the target
state, multiplying the nonbonded electrostatic interactions and van
der Waals interactions between A.sub.i and A.sub.j obtained in (a)
by coupling parameters .lamda..sub.elecAex and .lamda..sub.vdwAex,
respectively, wherein both of .lamda..sub.elecAex and
.lamda..sub.vdwAex are 0 at the reference state, 1 at the target
state, and varied from 0 to 1 at each intermediate state along the
transformation from the reference state to the target state; if
A.sub.a and A.sub.b are not valence-bonded to each other in the
reference state and valence-bonded in the target state, and the
nonbonded interactions between A.sub.i and A.sub.j are not excluded
in the reference state but excluded in the target state,
multiplying the nonbonded electrostatic interactions and van der
Waals interactions between A.sub.i and A.sub.j obtained in (a) by
coupling parameters .lamda..sub.elecBex and .lamda..sub.vdwBex
respectively, wherein both of .lamda..sub.elecBex and
.lamda..sub.vdwBex are 1 at the reference state, 0 at the target
state, and varied from 1 to 0 at each intermediate state along the
transformation from the reference state to the target state; and
(c) including the calculated nonbonded electrostatic interactions
and van der Waals interactions obtained in (b) into the total
energy of the simulation step of the corresponding system
state.
12. The method of claim 11, wherein performing molecular
simulations for all of the states further includes: if A.sub.a and
A.sub.b are valence-bonded to each other in the reference state and
not valence-bonded in the target state, and the nonbonded
interactions between A.sub.i and A.sub.j are excluded in the
reference state but not excluded in the target state, varying each
of .lamda..sub.elecAex and .lamda..sub.vdwAex according to a
schedule for each of the intermediate states along the
transformation from the reference state to the target state such
that when .lamda..sub.vdwAex is smaller than 1 for an intermediate
state, .lamda..sub.elecAex is 0 for that intermediate state; and if
A.sub.a and A.sub.b are not valence-bonded to each other in the
reference state and valence-bonded in the target state, and the
nonbonded interactions between A.sub.i and A.sub.j are not excluded
in the reference state but excluded in the target state, varying
each of .lamda..sub.elecBex and .lamda..sub.vdwBex according to a
schedule for each of the intermediate states along the
transformation from the reference state to the target state such
that when .lamda..sub.vdwAex is smaller than 1 for an intermediate
state, .lamda..sub.elecAex is 0 for that intermediate state.
13. The method of claim 11, wherein performing molecular
simulations for all of the states further includes: (a) computing
nonbonded electrostatic 1-4 pair interactions and van der Waals 1-4
pair interactions, using applicable parameters for electrostatic
1-4 pair interactions and van der Waals 1-4 pair interactions of a
force field, between two atoms A.sub.i and A.sub.j which together
with another two intervening atoms forms a bonded dihedral angle
interaction in either the reference state or the target state; (b)
if A.sub.a and A.sub.b are valence-bonded to each other in the
reference state and not valence-bonded in the target state: if the
nonbonded 1-4 pair interactions between A.sub.i and A.sub.j are
included in the reference state but not included in the target
state, multiplying the nonbonded electrostatic 1-4 pair
interactions and van der Waals 1-4 pair interactions between
A.sub.i and A.sub.j obtained in (a) by coupling parameters
.lamda..sub.elecA14 and .lamda..sub.vdwA14, respectively, wherein
both of .lamda..sub.elecA14 and .lamda..sub.vdwA14 are 1 at the
reference state, 0 at the target state, and varied from 1 to 0 at
each intermediate state along the transformation from the reference
state to the target state, and if the nonbonded 1-4 pair
interactions between A.sub.i and A.sub.j are not included in the
reference state but included in the target state, multiplying the
nonbonded electrostatic 1-4 pair interactions and van der Waals 1-4
pair interactions between A.sub.i and A.sub.j obtained in (a) by
coupling parameters .lamda..sub.elecB14 and .lamda..sub.vdwB14,
respectively, wherein both of .lamda..sub.elecB14 and
.lamda..sub.vdwB14 are 0 at the reference state, 1 at the target
state, and varied from 1 to 0 at each intermediate state along the
transformation from the reference state to the target state, and if
A.sub.a and A.sub.b are not valence-bonded to each other in the
reference state and valence-bonded in the target state: if the
nonbonded 1-4 pair interactions between A.sub.i and A.sub.j are not
included in the reference state but included in the target state,
multiplying the nonbonded electrostatic 1-4 pair interactions and
van der Waals 1-4 pair interactions between A.sub.i and A.sub.j
obtained in (a) by coupling parameters .lamda..sub.elecB14 and
.lamda..sub.vdwB14, respectively, wherein both of
.lamda..sub.elecB14 and .lamda..sub.vdwB14 are 0 at the reference
state, 1 at the target state, and varied from 0 to 1 at each
intermediate state along the transformation from the reference
state to the target state, and if the nonbonded 1-4 pair
interactions between A.sub.i and A.sub.j are included in the
reference state but not included in the target state, multiplying
the nonbonded electrostatic 1-4 pair interactions and van der Waals
1-4 pair interactions between A.sub.i and A.sub.j obtained in (a)
by the coupling parameters .lamda..sub.elecA14 and
.lamda..sub.vdwA14, respectively, wherein both of
.lamda..sub.elecA14 and .lamda..sub.vdwA14 are 1 at the reference
state, 0 at the target state, and varied from 1 to 0 at each
intermediate state along the transformation from the reference
state to the target state; and (c) including the calculated
electrostatic 1-4 pair interactions and van der Waals 1-4 pair
interactions obtained in (b) into the total energy of the
simulation step of the corresponding system state.
14. The method of claim 13, wherein performing molecular
simulations for all of the states further includes: if A.sub.a and
A.sub.b are valence-bonded to each other in the reference state and
not valence-bonded in the target state, if the nonbonded 1-4 pair
interactions between A.sub.i and A.sub.j are included in the
reference state but not included in the target state, varying each
of .lamda..sub.elecA14 and .lamda..sub.vdwA14 according to a
schedule for each of the intermediate states along the
transformation from the reference state to the target state such
that when .lamda..sub.vdwA14 is smaller than 1 for an intermediate
state, .lamda..sub.elecA14 is 0 for that intermediate state; if the
nonbonded 1-4 pair interactions between A.sub.i and A.sub.j are not
included in the reference state but included in the target state,
varying each of .lamda..sub.elecB14 and .lamda..sub.vdwB14
according to a schedule for each of the intermediate states along
the transformation from the reference state to the target state
such that when .lamda..sub.vdwB14 is smaller than 1 for an
intermediate state, .lamda..sub.elecB14 is 0 for that intermediate
state; and if A.sub.a and A.sub.b are not valence-bonded to each
other in the reference state and valence-bonded in the target
state, if the nonbonded 1-4 pair interactions between A.sub.i and
A.sub.j are not included in the reference state but included in the
target state, varying each of .lamda..sub.elecB14 and
.lamda..sub.vdwB14 according to a schedule for each of the
intermediate states along the transformation from the reference
state to the target state such that when .lamda..sub.vdwB14 is
smaller than 1 for an intermediate state, .lamda..sub.elecB14 is 0
for that intermediate state; if the nonbonded 1-4 pair interactions
between A.sub.i and A.sub.j are included in the reference state but
not included in the target state, varying each of
.lamda..sub.elecA14 and .lamda..sub.vdwA14 according to a schedule
for each of the intermediate states along the transformation from
the reference state to the target state such that when
.lamda..sub.vdwA14 is smaller than 1 for an intermediate state,
.lamda..sub.elecA14 is 0 for that intermediate state.
15. The method of any of claims 11-14, wherein at least one of the
computing of the van der Waals interactions includes using a
soft-core LJ interaction potential.
16. The method of claim 1, wherein at least one of the reference
state and the target state includes a molecule having a ring
structure in which the atoms A.sub.a and A.sub.b are bonded to each
other and form a part of the ring structure.
17. The method of claim 1, wherein the molecular simulations
include at least one of molecular dynamic simulations and Monte
Carlo simulations.
18. The method of claim 1, wherein calculating the free energy
difference between the reference state and the target state
comprises performing an analysis of the ensembles of micro-states
obtained at the target state, the reference state, and the
intermediate states by way of a determination and analysis of the
work associated with the variation of coupling parameter
.lamda..
19. The method of claim 1, wherein calculating the free energy
difference between the reference state and the target state
comprises performing an analysis of the ensembles of micro-states
obtained at the target state, the reference state, and the
intermediate states by way of an analysis of the differences in a
thermodynamic property of an ensemble of the micro-states obtained
at the target state, the reference state, and the intermediate
states as coupling parameter .lamda. is instantaneously varied for
the selected ensemble of micro-states.
20. The method of claim 19, wherein the ensemble is selected from
an NVT ensemble, a NPT ensemble, a NVE ensemble, and a .mu.VT
ensemble.
21. The method of claim 19, wherein performing the analysis of the
differences in a thermodynamic property comprises applying an
estimator selected from BAR, MBAR, WHAM, Zwanzig average
estimators.
22. The method of claim 19, wherein performing the analysis of the
differences in a thermodynamic property comprises applying one of
an FEP-family estimators.
23. The method of claim 1, wherein calculating the free energy
difference between the reference state and the target state
comprises performing a thermodynamic integration analysis of the
derivative of a thermodynamic property of an ensemble of
micro-states obtained for the target state, the reference state,
and the intermediate states with respect of the coupling parameter
.lamda..
24. The method of claim 23, wherein the ensemble is selected from
an NVT ensemble, a NPT ensemble, a NVE ensemble, and a .mu.VT
ensemble.
25. The method of claim 1, wherein at least one of P.sub.A and
P.sub.B is null.
26. A nontransitory computer readable medium, storing instructions
which when executed by one or more processors, carry out the method
of computing the free energy difference between a reference state
and a target state, wherein the reference state and target state
each include a common set of atoms P.sub.AB, and wherein the
reference state further includes a set of atoms P.sub.A, the target
state further includes a set of atoms P.sub.B, the set P.sub.A
being present only in the reference state and not in the target
state, and the set P.sub.B being present only in the target state
and not the reference state, where there exist at least two atoms
A.sub.a and A.sub.b, A.sub.a and A.sub.b being either: (1) not
valence-bonded to each other in the reference state and
valence-bonded in the target state, or (2) valence-bonded to each
other in reference state and not valence-bonded to each other in
the target state, the method comprising: (a) providing a topology,
including the bonded connections between the atoms and the relative
spatial arrangements of the atoms, for all the atoms in P.sub.A,
P.sub.B, and P.sub.AB; (b) determining one or more intermediate
system states along a transformation path between the reference
state and the target state, the transformation path defined by a
coupling parameter .lamda. that modulates the energies arising from
inter-atom interactions for each system state, the coupling
parameter .lamda. including a plurality of components each having a
value belonging to [0,1] and modulating a different type of
interaction energy; (c) performing, using at least one computer
processor, molecular simulations to obtain ensembles of the
micro-states for the reference state, the target state, and the
intermediate states, wherein performing molecular simulations for
each of the system states includes calculating a bonded stretch
interaction energy between the atoms A.sub.a and A.sub.b, the
bonded stretch interaction energy being defined by a soft bond
potential, wherein the soft bond potential is a function of a
bonded stretch component, .lamda..sub.sbs, of the coupling
parameter .lamda., and does not include any singular regions for
all values of .lamda..sub.sbs within [0,1] and for all values of
the distance r between A.sub.a and A.sub.b, the soft bond potential
further satisfies the following conditions: when .lamda..sub.sbs is
within (0,1), the soft bond potential is flat when the distance
between A.sub.a and A.sub.b approaches infinity; when A.sub.a and
A.sub.b are not valence bonded in either the reference state or the
target state, the soft bond potential is flat and zero for all
distances between A.sub.a and A.sub.b; and when A.sub.a and A.sub.b
are valence bonded in either the target state or the reference
state, the soft bond potential reverts to a harmonic potential; and
wherein the soft bond potential is wherein the soft bond potential
is a function of (r-r.sub.0).sup.2, where r.sub.0 is the
equilibrium distance between A.sub.a and A.sub.b, and is expressed
by: U sbs ( r , .lamda. sbs ) = 1 2 kf ( .lamda. sbs ) ( r - r 0 )
2 1 1 + g ( .lamda. sbs ) .alpha. ( k , .lamda. sbs ) ( r - r 0 ) 2
##EQU00032## where k is a constant, and the functions f, g and
.alpha. are each continuous functions and satisfy the following
conditions: f(.lamda..sub.sbs=0)=0, f(.lamda..sub.sbs=1)=1,
g(.lamda..sub.sbs=0)=1, g(.lamda..sub.sbs=1)=0,
.alpha.(k,.lamda..sub.sbs<1)>0; and (d) calculating, using at
least one computer processor, the free energy difference between
the reference state and the target state, by way of an analysis of
the ensembles of micro-states obtained at the target state, the
reference state, and the intermediate states.
27. The method of claim 1, wherein valence-bonded is
covalence-bonded.
28. The nontransitory computer readable medium storing instructions
of claim 26, wherein valence-bonded is covalence-bonded.
Description
BACKGROUND
[0001] Free energy is a fundamental molecular property that plays
an essential role in characterizing chemical and biological
systems. An understanding of the free energy behavior of many
chemical and biochemical processes, such as protein-ligand binding,
can be of critical importance in endeavors such as rational drug
design (which involves the design of small molecules that bind to a
biomolecular target).
[0002] Computer modeling and simulations are often used in free
energy studies. In most instances, evaluation of accurate absolute
free energies from simulations is extremely difficult, if at all
possible. Hence, the free energy difference between two
well-delineated thermodynamic states, or relative free energy, are
often used as a study system to provide insight to particular
systems, such as a relative binding affinity of a ligand predicated
on the measured affinity of a different but similar ligand (e.g., a
congeneric ligand).
[0003] In the relative free energy calculations, the two
thermodynamic states can be referred to as a reference system state
and a target system state, which can represent respectively an
initial state of a molecular system, such as a first molecule, and
an ending state of the molecule after one or more transformations
have taken place (such as a conformational change, topological
change, or a replacement of one atom or chemical group with another
(i.e., a mutation)). It is noted that such transformations may not
always represent realistic physical transformations, but may
involve nonphysical or "alchemical" transformations. Different
frameworks have been developed for calculating free energy
differences, such as free energy perturbations (FEP), thermodynamic
integrations (TI), and umbrella sampling.
[0004] Within the FEP framework, the free energy difference
.DELTA.F.sub.a.fwdarw.b between the two system states a and b can
be expressed by:
.DELTA. F a -> b = - 1 .beta. ln exp { - .beta. [ H b ( x , p x
) - H a ( x , p x ) ] } a ( 1 ) ##EQU00001##
where .beta..sup.-1=k.sub.BT, where is k.sub.B is the Boltzmann
constant, T is the temperature. .sub.a(x,p.sub.x) and .sub.b
(x,p.sub.x) are the Hamiltonians characteristic of states a and b
respectively. . . . .sub.a denotes an ensemble average over
configurations representative of the initial, reference state,
a.
[0005] In practical applications of FEP, the transformation between
the two thermodynamic states is usually achieved by a series of
transformations between non-physical, intermediate states along a
well-delineated pathway that connects a to b. This pathway is often
characterized by a general extent parameter, often referred to as a
coupling parameter, .lamda., which varies from 0 to 1 from the
reference state to the target state, and relates the Hamiltonians
of the two states by:
(.lamda.)=(1-.lamda.).sub.a+.sub.b (2)
where (.lamda.) is the .lamda.-coupled or hybrid Hamiltonian of the
system between the two states (including the two states, when
.lamda. takes the end values of 0 and 1). Hence, the free energy
difference .DELTA.F.sub.a.fwdarw.b between a and b will be:
.DELTA. F a -> b = - 1 .beta. ln exp { - .beta. [ H ( .lamda. =
1 ) - H ( .lamda. = 0 ) ] } .lamda. = 0 = - 1 .beta. i = 0 N - 1 ln
exp { - .beta. [ H ( x , p x ; .lamda. i + 1 ) - H ( x , p x ;
.lamda. i ) ] } i ( 3 ) ##EQU00002##
where N stands for the number of "windows" between neighboring
states between the reference (initial) state and the target (final)
state, and .lamda..sub.i is the values of the coupling parameter in
the initial, intermediate, and final state.
[0006] The free energy difference between the reference system
state a and the target system state b can also be calculated using
thermodynamic integration method, where the free energy difference
is calculated using the following formula:
.DELTA. F a .fwdarw. b = .intg. .lamda. = 0 .lamda. = 1 .lamda.
.differential. ( .lamda. ) .differential. .lamda. .lamda. ( 4 )
##EQU00003##
where .lamda. is the coupling parameter which varies from 0 to 1
from the reference state to the target state, (.lamda.) is the
.lamda.-coupled or hybrid Hamiltonian of the system between the two
states (including the two states, when .lamda. takes the end values
of 0 and 1), and
.differential. ( .lamda. ) .differential. .lamda. ##EQU00004##
is the first derivative of the coupled Hamiltonian with respect to
the coupling parameter .lamda.. In practical applications of TI,
the transformation between the reference system state and the
target system state is achieved by a series transformations along a
well-delineated pathway that connects a to b, and the ensemble
average of
.differential. ( .lamda. ) .differential. .lamda. ##EQU00005##
is calculated for all the states sampled, including the reference
system state, the intermediate non-physical states, and the target
system state. The free energy difference between the reference
system state and the target system state is then approximated by
numerical integration of the above integral based on the value of
the
.differential. ( .lamda. ) .differential. .lamda. .lamda. i ,
##EQU00006##
where .lamda..sub.i is the values of the coupling parameter in the
initial, intermediate, and final states.
[0007] Under conventional methods, calculating the free energy to
open a ring of a molecule into a linear structure or close a linear
structure of a molecule to form a ring can be difficult. In these
cases, the reversible work of turning on and off a valence bond
connecting two ring atoms needs to be computed. Although one
possible approach might be to annihilate a whole ring and grow a
corresponding linear structure from scratch, it is computationally
very inefficient and sometimes impossible when the ring is very
large or the ring is fused with other rings.
SUMMARY OF THE INVENTION
[0008] In some embodiments of the invention, a computer-implemented
method for computing the free energy difference between a reference
state and a target state is provided. The reference state and
target state each include a common set of atoms P.sub.AB. The
reference state further includes a set of atoms P.sub.A, and the
target state further includes a set of atoms P.sub.B. The set
P.sub.A is present only in the reference state and not in the
target state, and the set P.sub.B being present only in the target
state and not the reference state. There exist at least two atoms
A.sub.a and A.sub.b, A.sub.a and A.sub.b are either: (1) not
valence-bonded to each other in the reference state and
valence-bonded in the target state, or (2) valence-bonded to each
other in reference state and not valence-bonded to each other in
the target state. The method comprises:
[0009] (a) providing a topology, including the bonded connections
between the atoms and the relative spatial arrangements of the
atoms, for all the atoms in P.sub.A, P.sub.B, and P.sub.AB;
[0010] (b) determining one or more intermediate system states along
a transformation path between the reference state and the target
state, the transformation path defined by a coupling parameter
.lamda. that modulates the energies arising from inter-atom
interactions for each system state, the coupling parameter .lamda.
including a plurality of components each having a value belonging
to [0,1] and modulating a different type of interaction energy;
[0011] (c) performing, using at least one computer processor,
molecular simulations, such as Monte Carlo simulations or molecular
dynamics simulations, to obtain ensembles of micro-states for the
reference state, the target state, and the intermediate states,
wherein performing molecular simulations for each of the system
states includes calculating a bonded stretch interaction energy
between the atoms A.sub.a and A.sub.b, the bonded stretch
interaction energy being defined by a soft bond potential; and
[0012] (d) calculating, using at least one computer processor, the
free energy difference between the reference state and the target
state, by way of an analysis of the ensembles of micro-states
obtained at the target state, the reference state, and the
intermediate states.
[0013] In some embodiments of the method, the soft bond potential
is a function of a bonded stretch component, .lamda..sub.sbs, of
the coupling parameter .lamda., and does not include any singular
regions for all values of .lamda..sub.sbs within [0,1] and for all
values of the distance r between A.sub.a and A.sub.b. The soft bond
potential further satisfy the following conditions: when
.lamda..sub.sbs is within (0,1), the soft bond potential is flat
when the distance between A.sub.a and A.sub.b approaches infinity;
when A.sub.a and A.sub.b are not valence bonded in either the
reference state or the target state, the soft bond potential is
flat and zero for all distances between A.sub.a and A.sub.b; and
when A.sub.a and A.sub.b are valence bonded in either the target
state or the reference state, the soft bond potential reverts to a
harmonic potential.
[0014] In one embodiment of the method, the first derivative of the
soft bond potential with respect to .lamda..sub.sbs is continuous
and bounded for all values of .lamda..sub.sbs between 0 and 1 and
for all values of the distance r between A.sub.a and A.sub.b. In
another embodiment, both the first derivative and the second
derivative of the soft bond potential with respect to the distance
r between A.sub.a and A.sub.b are continuous and bounded for all
values of .lamda..sub.sbs, and approach zero when r approaches
infinity.
[0015] In some embodiments of the method, the soft bond potential
is a function of (r-r.sub.0).sup.2, where r.sub.0 is the
equilibrium distance between A.sub.a and A.sub.b. In one
embodiment, the soft bond potential is expressed by:
U sbs ( r , .lamda. sbs ) = 1 2 kf ( .lamda. sbs ) ( r - r 0 ) 2 1
1 + g ( .lamda. sbs ) .alpha. ( k , .lamda. sbs ) ( r - r 0 ) 2
##EQU00007##
where k is a constant, and the functions f, g and .alpha. are each
continuous functions and satisfy the following conditions:
f(.lamda..sub.sbs=0)=0, f(.lamda..sub.sbs=1)=1,
g(.lamda..sub.sbs=0)=1, g(.lamda..sub.sbs=1)=0,
.alpha.(k,.lamda..sub.sbs<1)>0. In particular embodiments,
f(.lamda..sub.sbs)=.lamda..sub.sbs,
g(.lamda..sub.sbs)=1-.lamda..sub.sbs, and .alpha.(k,
.lamda..sub.sbs) is a nonnegative constant.
[0016] In some embodiments, performing molecular simulations for
each of the system states comprises: if A.sub.a and A.sub.b are
valence-bonded to each other in the reference state and not
valence-bonded in the target state, using a schedule of
.lamda..sub.sbsA and a corresponding soft bond potential for
calculating the bonded stretch interaction energy between A.sub.a
and A.sub.b for each of the intermediate states, wherein
.lamda..sub.sbsA is 1 at the reference state, 0 at the target
state, and varied from 1 to 0 at each intermediate state along the
transformation from the reference state to the target state; and if
A.sub.a and A.sub.b are not valence-bonded to each other in the
reference state and valence-bonded in the target state, using a
schedule of .lamda..sub.sbsB and a corresponding soft bond
potential for each of the intermediate states and a soft bond
potential corresponding to the .lamda..sub.sbsB for calculating the
bonded stretch interaction between A.sub.a and A.sub.b, wherein
.lamda..sub.sbsB is 0 at the reference state, 1 at the target
state, and varied from 0 to 1 at each intermediate state along the
transformation from the reference state to the target state.
[0017] In some embodiments of the method, performing molecular
simulations for each of the system states further comprises:
[0018] (a) computing a bonded angle interaction, using applicable
parameters for bonded angle interactions of a force field, between
(i) a bond formed by A.sub.a and another atom A.sub.c, and (ii) the
bond between A.sub.a and A.sub.b that is being broken or formed by
the transformation from the reference state to the target
state;
[0019] (b) if A.sub.a and A.sub.b are valence-bonded to each other
in the reference state and not valence-bonded to each other in the
target state, multiplying the computed bonded angle interaction
obtained in (a) by a bonded angle coupling parameter
.lamda..sub.baA, wherein .lamda..sub.baA is 1 at the reference
state, 0 at the target state, and varied from 1 to 0 at each
intermediate state along the transformation from the reference
state to the target state; and
[0020] if A.sub.a and A.sub.b are not valence-bonded to each other
in the reference state and valence-bonded to each other in the
target state, multiplying the computed bonded angle interaction
obtained in (a) by a bonded angle coupling parameter
.lamda..sub.baB, wherein .lamda..sub.baB is 0 at the reference
state, 1 at the target state, and varied from 0 to 1 at each
intermediate state along the transformation from the reference
state to the target state; and
[0021] (c) including the bonded angle interaction obtained in (b)
into the total energy of a simulation step of the corresponding
system state,
[0022] In some embodiments of the method, performing molecular
simulations for each of the system states further includes:
[0023] (a) computing a dihedral angle interaction, using applicable
parameters for dihedral interactions of a force field, of a group
of four connected atoms {A.sub.i, A.sub.j, A.sub.k, A.sub.l}, the
group including both A.sub.a and A.sub.b;
[0024] (b) if A.sub.a and A.sub.b are valence-bonded to each other
in the reference state and not valence-bonded to each other in the
target state, multiplying the computed dihedral interaction
obtained in (a) by a dihedral angle coupling parameter
.lamda..sub.bdA, wherein .lamda..sub.bdA is 1 at the reference
state, 0 at the target state, and varied from 1 to 0 at each
intermediate state along the transformation from the reference
state to the target state; and
[0025] if A.sub.a and A.sub.b are not valence-bonded to each other
in the reference state and valence-bonded to each other in the
target state, multiplying the computed dihedral interaction
obtained in (a) by a dihedral angle coupling parameter
.lamda..sub.bdB, wherein .lamda..sub.bdB is 0 at the reference
state, 1 at the target state, and varied from 0 to 1 at each
intermediate state along the transformation from the reference
state to the target state; and
[0026] (c) including the dihedral interaction obtained in (b) into
the total energy of the simulation step of the corresponding system
state.
[0027] In the above embodiments, if A.sub.a and A.sub.b are
valence-bonded to each other in the reference state and not
valence-bonded to each other in the target state, the bonded angle
interaction and the bonded dihedral interaction coupling parameters
.lamda..sub.baA and .lamda..sub.bdA, can be each selected to be 0
when .lamda..sub.sbsA is smaller than a predefined threshold, and
if A.sub.a and A.sub.b are not valence-bonded to each other in the
reference state and valence-bonded to each other in the target
state, the bonded angle interaction and the bonded dihedral
interaction coupling parameters .lamda..sub.baB and .lamda..sub.bdB
can be each selected to be 0 when .lamda..sub.sbsB is smaller than
a predefined threshold.
[0028] In some embodiments of the method, performing molecular
simulations for all of the states further includes:
[0029] (a) computing nonbonded electrostatic interactions and van
der Waals interactions, using applicable parameters for
electrostatic interactions and van der Waals interactions of a
force field, between two atoms A.sub.i and A.sub.j and the
non-bonded exclusion status of the pair (A.sub.i, A.sub.j) is
affected by the transformation from the reference state to the
target state;
[0030] (b) if A.sub.a and A.sub.b are valence-bonded to each other
in the reference state and not valence-bonded in the target state,
and the nonbonded interactions between A.sub.i and A.sub.j are
excluded in the reference state but not excluded in the target
state, multiplying the nonbonded electrostatic interactions and van
der Waals interactions between A.sub.i and A.sub.j obtained in (a)
by coupling parameters .lamda..sub.elecAex and .lamda..sub.vdwAex,
respectively, wherein both of .lamda..sub.elecAex and
.lamda..sub.vdwAex are 0 at the reference state, 1 at the target
state, and varied from 0 to 1 at each intermediate state along the
transformation from the reference state to the target state;
[0031] if A.sub.a and A.sub.b are not valence-bonded to each other
in the reference state and valence-bonded in the target state, and
the nonbonded interactions between A.sub.i and A.sub.j are not
excluded in the reference state but excluded in the target state,
multiplying the nonbonded electrostatic interactions and van der
Waals interactions between A.sub.i and A.sub.j obtained in (a) by
coupling parameters .lamda..sub.elecBex and .lamda..sub.vdwBex,
respectively, wherein both of .lamda..sub.elecBex and
.lamda..sub.vdwBex are 1 at the reference state, 0 at the target
state, and varied from 1 to 0 at each intermediate state along the
transformation from the reference state to the target state;
and
[0032] (c) including the calculated nonbonded electrostatic
interactions and van der Waals interactions obtained in (b) into
the total energy of the simulation step of the corresponding system
state.
[0033] In some embodiments of the method, wherein performing
molecular simulations for all of the states further includes:
[0034] if A.sub.a and A.sub.b are valence-bonded to each other in
the reference state and not valence-bonded in the target state, and
the nonbonded interactions between A.sub.i and A.sub.j are excluded
in the reference state but not excluded in the target state,
varying each of .lamda..sub.elecAex and .lamda..sub.vdwAex
according to a schedule for each of the intermediate states along
the transformation from the reference state to the target state
such that when .lamda..sub.vdwAex is smaller than 1 for an
intermediate state, .lamda..sub.elecAex is 0 for that intermediate
state; and
[0035] if A.sub.a and A.sub.b are not valence-bonded to each other
in the reference state and valence-bonded in the target state, and
the nonbonded interactions between A.sub.i and A.sub.j are not
excluded in the reference state but excluded in the target state,
varying each of .lamda..sub.elecBex and .lamda..sub.vdwBex
according to a schedule for each of the intermediate states along
the transformation from the reference state to the target state
such that when .lamda..sub.vdwAex is smaller than 1 for an
intermediate state, .lamda..sub.elecAex is 0 for that intermediate
state.
[0036] In some embodiments of the method, performing molecular
simulations for all of the states further includes:
[0037] (a) computing nonbonded electrostatic 1-4 pair interactions
and van der Waals 1-4 pair interactions, using applicable
parameters for electrostatic 1-4 pair interactions and van der
Waals 1-4 pair interactions of a force field, between two atoms
A.sub.i and A.sub.j which together with another two intervening
atoms forms a bonded dihedral angle interaction in either the
reference state or the target state;
[0038] (b) if A.sub.a and A.sub.b are valence-bonded to each other
in the reference state and not valence-bonded in the target state:
[0039] if the nonbonded 1-4 pair interactions between A.sub.i and
A.sub.j are included in the reference state but not included in the
target state, multiplying the nonbonded electrostatic 1-4 pair
interactions and van der Waals 1-4 pair interactions between
A.sub.i and A.sub.j obtained in (a) by coupling parameters
.lamda..sub.elecA14 and .lamda..sub.vdwA14, respectively, wherein
both of .lamda..sub.elecA14 and .lamda..sub.vdwA14 are 1 at the
reference state, 0 at the target state, and varied from 1 to 0 at
each intermediate state along the transformation from the reference
state to the target state, and [0040] if the nonbonded 1-4 pair
interactions between A.sub.i and A.sub.j are not included in the
reference state but included in the target state, multiplying the
nonbonded electrostatic 1-4 pair interactions and van der Waals 1-4
pair interactions between A.sub.i and A.sub.j obtained in (a) by
coupling parameters .lamda..sub.elecB14 and .lamda..sub.vdwB14,
respectively, wherein both of .lamda..sub.elecB14 and
.lamda..sub.vdwB14 are 0 at the reference state, 1 at the target
state, and varied from 1 to 0 at each intermediate state along the
transformation from the reference state to the target state,
and
[0041] if A.sub.a and A.sub.b are not valence-bonded to each other
in the reference state and valence-bonded in the target state:
[0042] if the nonbonded 1-4 pair interactions between A.sub.i and
A.sub.j are not included in the reference state but included in the
target state, multiplying the nonbonded electrostatic 1-4 pair
interactions and van der Waals 1-4 pair interactions between
A.sub.i and A.sub.j obtained in (a) by coupling parameters
.lamda..sub.elecB14 and .lamda..sub.vdwB14, respectively, wherein
both of .lamda..sub.elecB14 and .lamda..sub.vdwB14 are 0 at the
reference state, 1 at the target state, and varied from 0 to 1 at
each intermediate state along the transformation from the reference
state to the target state, and [0043] if the nonbonded 1-4 pair
interactions between A.sub.i and A.sub.j are included in the
reference state but not included in the target state, multiplying
the nonbonded electrostatic 1-4 pair interactions and van der Waals
1-4 pair interactions between A.sub.i and A.sub.j obtained in (a)
by the coupling parameters .lamda..sub.elecA14 and
.lamda..sub.vdwA14, respectively, wherein both of
.lamda..sub.elecA14 and .lamda..sub.vdwA14 are 1 at the reference
state, 0 at the target state, and varied from 1 to 0 at each
intermediate state along the transformation from the reference
state to the target state; and
[0044] (c) including the calculated electrostatic 1-4 pair
interactions and van der Waals 1-4 pair interactions obtained in
(b) into the total energy of the simulation step of the
corresponding system state.
[0045] In some embodiments of the method, performing molecular
simulations for all of the states further includes:
[0046] (a) if A.sub.a and A.sub.b are valence-bonded to each other
in the reference state and not valence-bonded in the target state,
[0047] if the nonbonded 1-4 pair interactions between A.sub.i and
A.sub.j are included in the reference state but not included in the
target state, varying each of .lamda..sub.elecA14 and
.lamda..sub.vdwA14 according to a schedule for each of the
intermediate states along the transformation from the reference
state to the target state such that when .lamda..sub.vdwA14 is
smaller than 1 for an intermediate state, .lamda..sub.elecA14 is 0
for that intermediate state; [0048] if the nonbonded 1-4 pair
interactions between A.sub.i and A.sub.j are not included in the
reference state but included in the target state, varying each of
.lamda..sub.elecB14 and .lamda..sub.vdwB14 according to a schedule
for each of the intermediate states along the transformation from
the reference state to the target state such that when
.lamda..sub.vdwB14 is smaller than 1 for an intermediate state,
.lamda..sub.elecB14 is 0 for that intermediate state
[0049] (b) if A.sub.a and A.sub.b are not valence-bonded to each
other in the reference state and valence-bonded in the target
state, [0050] if the nonbonded 1-4 pair interactions between
A.sub.i and A.sub.j are not included in the reference state but
included in the target state, varying each of .lamda..sub.elecB14
and .lamda..sub.vdwB14 according to a schedule for each of the
intermediate states along the transformation from the reference
state to the target state such that when .lamda..sub.vdwB14 is
smaller than 1 for an intermediate state, .lamda..sub.elecB14 is 0
for that intermediate state, [0051] if the nonbonded 1-4 pair
interactions between A.sub.i and A.sub.j are included in the
reference state but not included in the target state, varying each
of .lamda..sub.elecA14 and .lamda..sub.vdwA14 according to a
schedule for each of the intermediate states along the
transformation from the reference state to the target state such
that when .lamda..sub.vdwA14 is smaller than 1 for an intermediate
state, .lamda..sub.elecA14 is 0 for that intermediate state;
[0052] In some embodiments of the method, computing van der Waals
interactions can include using a soft-core LJ interaction
potential.
[0053] In some embodiments of the method, either the reference
state or the target state includes a molecule having a ring
structure in which the atoms A.sub.a and A.sub.b are bonded to each
other and form a part of the ring structure.
[0054] In some embodiments of the method, calculating the free
energy difference between the reference state and the target state
comprises performing an analysis of the ensemble of micro-states
obtained at the target state, the reference state, and the
intermediate states by way of a determination and analysis of the
work associated with the variation of coupling parameter
.lamda..
[0055] In some embodiments of the method, calculating the free
energy difference between the reference state and the target state
comprises performing an analysis of the ensemble of micro-states
obtained at the target state, the reference state, and the
intermediate states by way of an analysis of the differences in a
thermodynamic property of a suitable ensemble of the micro-states
obtained at the target state, the reference state, and the
intermediate states as the coupling parameter .lamda. is
instantaneously varied for the selected ensemble of micro-states.
Performing the analysis of the differences in the thermodynamic
property comprises applying an estimator selected from BAR, MBAR,
WHAM, Zwanzig average estimators, or one of an FEP-family
estimators. The ensemble can be, for example, an NVT ensemble, a
NPT ensemble, a NVE ensemble, and a .mu.VT ensemble.
[0056] In some embodiments of the method, calculating the free
energy difference between the reference state and the target state
comprises performing a thermodynamic integration analysis of the
derivative of a thermodynamic property of a suitable ensemble of
micro-states obtained for the target state, the reference state,
and the intermediate states with respect of the coupling parameter
.lamda.. The ensemble can be, for example, an NVT ensemble, a NPT
ensemble, a NVE ensemble, and a .mu.VT ensemble.
[0057] The invention also provides an apparatus including one or
more processors, a memory operably coupled to the one or more
processors comprising instructions executable by the processors,
the one or more processors being operable when executing the
instructions to perform the various embodiments of the method as
described herein. The invention further provides non-transitory
computer readable media storing the instructions which when
executed by one or more processors, carry out the various
embodiments of the method as described herein.
[0058] Other features and advantages of the invention are apparent
from the following description, and from the claims.
BRIEF DESCRIPTION OF THE DRAWINGS
[0059] The present invention will be better understood by reference
to the accompanying drawings, wherein:
[0060] FIGS. 1a-1c are diagrams showing illustrative ring-closing
transformations between a reference system state to a target system
state according to some embodiments of the present invention;
[0061] FIG. 2 depicts plots of a soft bond potential as a function
of inter-atom distance between the two atoms involved in a bond
formation or breaking transformation at different values of the
coupling parameter .lamda. according to some embodiments of the
present invention;
[0062] FIG. 3 depicts plots of a soft bond potential as a function
of inter-atom distance between the two atoms involved in a bond
formation or breaking transformation at different values of the
coupling parameter .lamda., as well as the derivative of the
potential with respect to the coupling parameter .lamda., and the
first and second derivative of the potential with respect to the
inter-particle distance according to some embodiments of the
present invention;
[0063] FIG. 4 is a diagram showing an illustrative ring-opening
transformation occurring in a cyclic structure, the free energy
change of which is amenable to the application of embodiments of
the present invention;
[0064] FIG. 5 is a diagram showing an illustrative bond breaking
between two connected ring structures, the free energy change of
which is amenable to the application of embodiments of the present
invention;
[0065] FIG. 6 is a depiction of a protein residue mutation
transformation where a non-proline residue in the reference system
is transformed into a proline residue in the target system, the
free energy change of which is amenable to the application of
embodiments of the present invention; and
[0066] FIG. 7 is a depiction of a protein residue deletion
transformation where a residue (proline) in the reference system is
deleted in the target system, resulting in a direct connection
between its two neighboring residues, the free energy change of
which is amenable to the application of embodiments of the present
invention.
DETAILED DESCRIPTION
[0067] The present application discloses computer-implemented
methods and systems for computing free energy difference between a
reference system state and a target system state. In particular, to
address the issues in determining free energy difference arising
from bond breaking in a ring structure that transforms a ring
structure to a linear structure and bond formation that transforms
a linear structure into a ring structure (each of which is further
discussed below), the methods and systems disclosed in the present
application utilize a functional form for bond stretching that
allows a rigorous connection to the harmonic bond functional form
to be maintained at any points in the alchemical transformation for
calculating the free energy difference. Accordingly, the methods
and systems of the present application can advantageously improve
numerical stability and accuracy of the free energy calculations.
However, it is to be recognized that the general principles of the
free energy calculations using such modified bond stretching
potentials disclosed herein can be applied generally in any bond
formation and breaking situations, and not limited to ring closing
or ring opening.
[0068] As with traditional free energy difference calculations, the
atoms in the system can be categorized into different groups for
evaluating the system energy in different system states. The
reference state and target state both include a common set of atoms
P.sub.AB. The reference state further includes a set of atoms
P.sub.A, and the target state further includes a set of atoms
P.sub.B. The set of atoms P.sub.A are present only in the reference
state and not in the target state, and the set of atoms P.sub.B are
present only in the target state and not the reference state. In a
ring formation scenario, P.sub.A can be the atoms connected to the
two terminal atoms to form a bond. During the course of the
transformation, the atoms in P.sub.A and in P.sub.B interact with
other atoms within their own set as well as with those in P.sub.AB,
but the atoms in P.sub.A do not interact with any atoms in P.sub.B,
or vice versa. For example, for a molecule having a structure shown
in FIG. 1a, where a bond is about to be formed between two atoms
A.sub.1 and A.sub.2, P.sub.AB={A.sub.1, A.sub.2, A.sub.3, A.sub.4,
A.sub.5, A.sub.6, A.sub.7, A.sub.8, A.sub.9} (the atoms are shown
as having valence of 4 for purpose of illustration only; it is
understood that other atoms or groups can be represented in this
structure), P.sub.A={A.sub.10, A.sub.11}, and P.sub.B=O. By way of
another example, for a molecule having a backbone structure
A.sub.1-A.sub.3-A.sub.2 shown in FIG. 1b, where a molecular
fragment including A.sub.12, A.sub.13, A.sub.14 is to be inserted
between A.sub.1 and A.sub.3 to form two bonds A.sub.1-A.sub.12 and
A.sub.3-A.sub.12, P.sub.AB={A.sub.1, A.sub.2, A.sub.3, A.sub.4,
A.sub.5, A.sub.6, A.sub.7, A.sub.8, A.sub.9}; P.sub.A={A.sub.10,
A.sub.11}, and P.sub.B={A.sub.12, A.sub.13, A.sub.14}. In the
latter example, the insertion of the fragment to form the two bonds
can be considered as taking two steps, the first being a linear
growth of the chain A.sub.3-A.sub.2-A.sub.1 on A.sub.1 (or A.sub.3)
by the fragment containing A.sub.12, and the second step being the
closure of the ring between A.sub.12-A.sub.3 (or A.sub.12-A.sub.1).
In another example as illustrated in FIG. 1c, where a molecule
having a closed 4-membered ring structure is alchemically
transformed into a 5-membered ring structure, P.sub.AB={A.sub.1,
A.sub.2, A.sub.3, A.sub.4, A.sub.5, A.sub.6, A.sub.7, A.sub.8,
A.sub.9, A.sub.10, A.sub.11, A.sub.12}, P.sub.A.sup.=O,
P.sub.B={A.sub.13, A.sub.14, A.sub.15}.
[0069] As described herein, the coupling Hamiltonian (.lamda.) for
alchemical transformation involving ring opening and closing can
generally include the following terms, .sub.bs(.lamda..sub.bs),
.sub.ba(.lamda..sub.ba), .sub.bd(.lamda..sub.bd) and
.sub.nb(.lamda..sub.nb), corresponding to the bonded stretch terms,
the bonded angle terms, the bonded dihedral angle terms, and the
nonbonded exclusion and 1-4 pair interaction terms respectively. To
simplify the discussion, in the following the case where a valence
bond between two ring atoms is being formed (i.e., the
transformation from the reference state to the target state
involves ring closing) is described in detail.
[0070] As an initial step, the topology of the system is provided,
including the bonded connections between the atoms in the system
and the relative spatial arrangements of the atoms forming each of
P.sub.AB, P.sub.A, and P.sub.B.
[0071] One or more, e.g., a plurality of intermediate states
between the reference state and the target state can be determined
along a path defined by different values of the coupling parameter
.lamda., where the increments of .lamda. in value move the system
from the reference state to the target state. While .lamda. can be
a scalar variable that varies from 0 to 1, in some embodiments of
the present invention, such as those further discussed below,
.lamda. can be a vector containing different components for
different types of interactions within the system. Computer
molecular simulations, such as, but not limited to, molecular
dynamics or Monte Carlo simulations, can be performed to obtain
ensembles of the micro-states for the reference state, the target
state, and each of the intermediate states. The .lamda. values of
the intermediate states can be chosen by known techniques such that
between each neighboring .lamda. windows on the "reaction pathway"
from the reference state to the target state there is substantial
overlap between the micro-states in the successive .lamda. windows
that are sampled by the molecular simulations.
[0072] In performing molecular simulations for all these states,
the bonded stretch interaction energy between the two atoms A.sub.a
and A.sub.b that are to form a bond (e.g., A.sub.1 and A.sub.3 in
FIG. 1a) can be defined by a soft bond potential which is modulated
by 2 (or the bond stretch component thereof). When .lamda.=0
(A.sub.a and A.sub.b are completely nonbonded in the reference
state), the soft bond potential is a flat potential for all
distances r between A.sub.a and A.sub.b. When 0<.lamda.<1,
(the bond between A.sub.a and A.sub.b is being "partially formed"
in the alchemical transformation), the soft bond potential levels
off to a flat potential when r.fwdarw..infin., i.e., the partial
derivative of the potential with respect to the distance r between
A.sub.a and A.sub.b is zero when r.fwdarw..infin.. When .lamda.=1
(A.sub.a and A.sub.b are fully valence bonded in the target state),
the soft bond potential reverts to a harmonic potential. Further,
the potential energy function for the bond stretch term does not
have any singular regions for all values of the bonded stretch
component, .lamda..sub.sbs, of the coupling parameter .lamda.
within [0,1] and for all values of the distance r between Aa and
Ab. The details of developing the soft bond potential and some
properties of the soft bond potential are provided below.
[0073] In popular molecular mechanics force fields, such as OPLS,
CHARMM, and AMBER, the bonded stretch interactions between two
atoms are modeled by a harmonic bond of the following form:
U bs ( .lamda. , r ) = 1 2 k ( r - r 0 ) 2 ( 5 ) ##EQU00008##
where k is the "force constant" or "Hookean constant" which defines
the strength (or rigidity) of the bond of the force field used, r
is the instantaneous distance between the two atoms and r.sub.0 is
the equilibrium distance between the two atoms. In conventional
method of linear scaling of the coupling parameter between the
Hamiltonians of the reference state and the target state, the
bonded stretch term has the following form in the coupling
Hamiltonian:
U bs ( .lamda. , r ) = 1 2 .lamda. k ( r - r 0 ) 2 ( 6 )
##EQU00009##
Therefore,
[0074] .differential. F ( .lamda. ) .differential. .lamda. =
.differential. U bs ( .lamda. , r ) .differential. .lamda. .lamda.
= 1 2 k ( r - r 0 ) 2 .lamda. ( 7 ) ##EQU00010##
The integrand in the above equation approaches infinity when r is
very large. In the limit when .lamda., approaches 0, there is no
bonded stretch interaction between the two atoms, and the distance
between the two atoms can be very large, leading to singularity and
numerical instability problem in the calculation of the above
integral. In practice, the distance between the two atoms is
limited by the size of the simulation box (a three-dimensional
volume or unit cell in which the simulation is conducted and
boundary conditions, such as periodic boundary conditions, can be
applied), the singularity problem can be avoided. However, since
the integrand in the above equation is unbounded, it can still
cause numerical instability and inaccuracy problems in the free
energy calculations.
[0075] Additionally, using the above conventional coupling
Hamiltonian functional form, the potential is undefined for
.lamda.=0 when r approaches infinity. The limiting value depends on
how .lamda. approaches 0 and how r approaches infinity, i.e.:
lim r .fwdarw. .infin. lim .lamda. .fwdarw. 0 U bs ( .lamda. , r )
= lim r .fwdarw. .infin. lim .lamda. .fwdarw. 0 1 2 .lamda. k ( r -
r 0 ) 2 diverges ( 8 ) ##EQU00011##
[0076] In order to obtain a pathway that allows stable and
efficient simulations from which reliable free energies involving
the annihilation and/or formation of a bond between two atoms can
be determined, the present inventors have discovered the following
coupling potential (referred to herein as the soft bond potential)
to connect the two physical systems where the harmonic interactions
between the two atoms are fully turned on and off when the coupling
parameter .lamda., changes between 0 and 1:
U sbs ( r , .lamda. ) = 1 2 kf ( .lamda. ) ( r - r 0 ) 2 1 1 + g (
.lamda. ) .alpha. ( k , .lamda. ) ( r - r 0 ) 2 ( 9 )
##EQU00012##
where the functions f, g and .alpha. are each continuous functions
and simultaneously satisfy the following conditions:
f(.lamda.=0)=0; f(.lamda.=1)=1; g(.lamda.=0)=1; g(.lamda.=1)=0;
.alpha.(k,.lamda.<1)>0. It is noted that for all the
discussions herein regarding the soft bond potential for the bonded
stretch interactions, .lamda. as used in the equations (from
Equation 9 and onwards) refers to the bonded stretch component,
.lamda..sub.sbs, of the coupling parameter .lamda..
[0077] In particular example of the soft bond potential described
by Eq. 9, f(.lamda.)=.lamda., g(.lamda.)=1-.lamda., and
.alpha.(k,.lamda.)=.alpha.=const (a constant number), i.e., the
soft bond potential takes the following form:
U sbs ( r , .lamda. ) = 1 2 .lamda. k ( r - r 0 ) 2 1 1 + .alpha. (
1 - .lamda. ) ( r - r 0 ) 2 ( 10 ) ##EQU00013##
[0078] It can be seen that the soft bond interaction of Eq. 10
correctly recovers the two physical end states when .lamda.=0 and
.lamda.=1:
U sbs ( r , .lamda. = 1 ) = 1 2 k ( r - r 0 ) 2 U sbs ( r , .lamda.
= 0 ) = 0 U sbs ( r .fwdarw. .infin. , .lamda. ) = k .lamda. 2
.alpha. ( 1 - .lamda. ) when .lamda. .noteq. 1 ( 11 )
##EQU00014##
[0079] The introduction of .alpha.(1-.lamda.)(r-r.sub.0).sup.2 in
the denominator of Eq. 10 changes the harmonic interaction into a
soft bond interaction (the interaction is bounded when r approaches
infinity) when .lamda. is smaller than 1, and at .lamda.=1 (the
bond is formed) the function has the exact harmonic potential
form.
[0080] The above functional form removes the singularity and
numerical instability problems associated with the conventional
harmonic potentials. In the following, some properties of the soft
bond interaction functional form as exemplified by Eq. 10 are
discussed.
Property 1:
[0081] The soft bond interaction functional form does not have any
singular regions for all values of .lamda. between 0 and 1 and for
all values of r. From the above description, i.e.:
U sbs ( r , .lamda. = 1 ) = 1 2 k ( r - r 0 ) 2 U sbs ( r , .lamda.
= 0 ) = 0 U sbs ( r , .lamda. ) = 1 2 .lamda. k ( r - r 0 ) 2 1 1 +
.alpha. ( 1 - .lamda. ) ( r - r 0 ) 2 for .lamda. .di-elect cons. (
0 , 1 ) ( 12 ) ##EQU00015##
[0082] The soft bond potential for a model system with the force
constant k=20 kcalmol.sup.-1.ANG..sup.-2 and the soft bond
parameter .alpha.=1 at different values of .lamda. are given in
FIG. 2. It can be seen clearly that the potential is continuous for
all values of .lamda. and r. It changes slowly from a harmonic
potential at .lamda.=1 to a soft bond potential at intermediate
values of .lamda., and goes to 0 at .lamda.=0.
Property 2:
[0083] .differential. F sbs ( .lamda. ) .differential. .lamda. =
.differential. U sbs ( .lamda. ) .differential. .lamda. .lamda. = 1
2 k ( r - r 0 ) 2 1 + .alpha. ( r - r 0 ) 2 ( 1 + ( 1 - .lamda. )
.alpha. ( r - r 0 ) 2 ) 2 .lamda. ( 13 ) ##EQU00016##
[0084] As discussed in the above section, it is
.differential. F ( .lamda. ) .differential. .lamda.
##EQU00017##
which determines the numerical stability and accuracy of the free
energy calculations. In the above formulation, when
.lamda..epsilon.[0,1), the thermodynamic property to be averaged in
the rightmost bracket in Eq. 13 is continuous and bounded for all
values of r. When .lamda.=1, the soft bond potential recovers the
harmonic potential, so only phase space regions where r is close to
r.sub.0 are sampled and taken into the ensemble average. In the
phase space regions where r is close to r.sub.0, the integrand in
Eq. 13 is also bounded. Therefore, the quantity in the bracket of
Eq. 13 is bounded for all values of .lamda. between 0 and 1.
Since
.differential. F ( .lamda. ) .differential. .lamda.
##EQU00018##
does not have any singular region for all values of .lamda. between
0 and 1, accurate and reliable free energy results can be obtained
using the above soft bond interaction functional form.
[0085] As an example, the derivative of the soft bond potential of
Eq. 10 with respect to the coupling parameter .lamda. at
.lamda.=0.5 for a model system with force constant k=20
kcalmol.sup.-1.ANG..sup.-2 and the soft bond parameter .alpha.=1 is
plotted in FIG. 3. It can be seen that
.differential. F ( .lamda. ) .differential. .lamda.
##EQU00019##
does not have any singular region for all values of .lamda. and r,
allowing reliable and accurate free energy calculations.
Property 3:
[0086] .differential. U sbs ( .lamda. , r ) .differential. r = k
.lamda. r - r 0 ( 1 + ( 1 - .lamda. ) .alpha. ( r - r 0 ) 2 ) 2
.differential. 2 U sbs ( .lamda. , r ) .differential. r 2 = k
.lamda. 1 - 3 ( 1 - .lamda. ) .alpha. ( r - r 0 ) 2 ( 1 + ( 1 -
.lamda. ) .alpha. ( r - r 0 ) 2 ) 3 ( 14 ) ##EQU00020##
[0087] Both the first derivative and the second derivative of the
soft bond potential of Eq. 10 with respect to the inter-particle
(or inter-atom) distance r are continuous and bounded for all
values of .lamda..epsilon.[0,1] and they are short ranged and
approaching 0 when r.fwdarw..infin.. When .lamda.=1, the soft bond
potential reverts to the harmonic potential, and only phase space
regions where r is close to r.sub.0 are sampled in the molecular
simulations. Thus the first and the second derivatives of the
potential with respect to the inter-particle distance r are also
continuous and bounded in the physically accessible phase space
regions when .lamda.=1. Therefore, the forces and acceleration on
the atoms are always continuous and bounded for all values of
.lamda. between 0 and 1, allowing stable molecular dynamics
simulations to be performed for all values of .lamda.. As an
example, the first and second derivative of the soft bond potential
with respect to the inter-particle distance r at .lamda.=0.5 for a
model system with force constant k=20 kcalmol.sup.-1.ANG..sup.-1
and the soft bond parameter .alpha.=1 are plotted in FIG. 3. It can
be seen that those derivatives are continuous over the whole space,
allowing stable molecular dynamics simulations.
[0088] It is clear from the above description that the soft bond
potential removes the singularity and numerical instability
problems associated with the traditional methods, and it allows
stable molecular dynamics simulations and more convergent Monte
Carlo to be performed for all system states, including the
reference system state, the intermediate system states, and the
target system state. Using the soft bond potential described
herein, the free energy difference between the reference system
state and the target system state involving breaking and form
valance bond can be accurately and reliable calculated. The free
energy calculations utilizing such a soft bond potential are
particularly advantageous for alchemical transformations involving
ring opening, ring closing, or ring rearrangement, where the
computational efficiency and convergence of the free energy
calculations can be significantly improved over the conventional
methods.
[0089] As described in above section, in the conventional coupling
Hamiltonian, the limit does not exist for the initial state where
.lamda.=0 when r approaches infinity, while in the soft bond
potential the limit exists for all values of r. The conventional
coupling Hamiltonian is a special case of the soft bond potential
where .alpha.=0. While the conventional coupling Hamiltonian and
the soft bond potential reach the same end state when .lamda.=1,
the initial states when .lamda.=0 are different depending on the
value of the soft bond parameter .alpha.. In the following, it is
shown that the free energy of the initial state does not depend on
the soft bond parameter .alpha..
[0090] Consider a Hamiltonian with the following potential energy
term (i.e., Eq. 10 where .lamda.=0):
U ( r , .alpha. , .lamda. ) = 1 2 k .lamda. ( r - r 0 ) 2 1 1 +
.alpha. ( r - r 0 ) 2 ( 15 ) ##EQU00021##
[0091] When .alpha.=0, it becomes the conventional harmonic
potential, and when .alpha..noteq.0 it becomes the soft bond
potential. We need to prove that
F(.alpha.=0,.lamda.=0)=F(.alpha.=.alpha.'>0,.lamda.=0) (16)
Note that
F ( .alpha. = .alpha. ' > 0 , .lamda. = 0 ) - F ( .alpha. = 0 ,
.lamda. = 0 ) = .intg. 0 .alpha. ' .alpha. lim .lamda. -> 0
.differential. F ( .alpha. , .lamda. ) .differential. .alpha. ( 17
) ##EQU00022##
where
.differential. F ( .alpha. , .lamda. ) .differential. .alpha. =
.differential. U ( r , .alpha. , .lamda. ) .differential. .alpha.
.alpha. = k .lamda. ( r - r 0 ) 4 2 ( 1 + .alpha. ( r - r 0 ) 2 ) 2
.alpha. = .lamda. k ( r - r 0 ) 4 2 ( 1 + .alpha. ( r - r 0 ) 2 ) 2
= .lamda. I ( .alpha. ) ( 18 ) ##EQU00023##
[0092] The ensemble average I(.alpha.) is a finite number since
I ( .alpha. ) = k ( r - r 0 ) 4 2 ( 1 + .alpha. ( r - r 0 ) 2 ) 2
.alpha. = .intg. 0 r ma x r k ( r - r 0 ) 4 2 ( 1 + .alpha. ( r - r
0 ) 2 ) 2 exp ( - .beta. U ( r , .alpha. , .lamda. ) ) .intg. 0 r
ma x exp ( - .beta. U ( r , .alpha. , .lamda. ) ) .ltoreq. .intg. 0
r ma x rk ( r - r 0 ) 4 1 exp ( - .beta. U ( r ma x , .alpha. ,
.lamda. ) ) V .ltoreq. kr ma x 4 exp ( - .beta. U ( r ma x ,
.alpha. , .lamda. ) ) V ( 19 ) ##EQU00024##
[0093] Therefore,
lim .lamda. -> 0 .differential. F ( .alpha. , .lamda. )
.differential. .alpha. = lim .lamda. -> 0 .lamda. I ( .alpha. )
= 0 F ( .alpha. = .alpha. ' > 0 , .lamda. = 0 ) - F ( .alpha. =
0 , .lamda. = 0 ) = 0 ( 20 ) ##EQU00025##
[0094] In view of the above, the initial state for the conventional
coupling Hamiltonian and soft bond potential when .lamda.=0 have
the same free energy.
[0095] In free energy calculations, the convergence of the free
energy can be affected by the overlap of phase space regions
between neighboring or successive .lamda. windows. Empirically, a
suitable path from the initial state to the final state can be
achieved when the change of the free energy as a function of
.lamda. is continuous and smooth for all values of .lamda.,
i.e.:
.differential. F sbs ( .lamda. ) .differential. .lamda. =
.differential. U sbs ( .lamda. ) .differential. .lamda. .lamda. = 1
2 k ( r - r 0 ) 2 1 + .alpha. ( r - r 0 ) 2 ( 1 + ( 1 - .lamda. )
.alpha. ( r - r 0 ) 2 ) 2 .lamda. .apprxeq. const ( 21 )
##EQU00026##
[0096] Thus, in preferred embodiments, the value of the soft bond
parameter a can be obtained when Eq. 21 is satisfied.
[0097] As mentioned above, in ring opening and closing free energy
calculations, in addition to the bonded stretch interactions, there
can be other interaction energy terms which are different in the
initial and final states, including the bonded angle terms, the
bonded dihedral angle terms, and the nonbonded interactions. These
different types of interactions can be treated differently during
the transformation process. As described herein, in some
embodiments, the coupling parameter .lamda. can take a vector form
which includes components for different types of interactions. For
example, the coupling vector .lamda. can include the following
components for interactions affected by the formation (or the
breaking) of the valence bond between the reference state and the
target state: a component .lamda..sub.sbs for bonded stretch energy
term; a component .lamda..sub.ba for the bonded angle energy term;
a component .lamda..sub.bd for bonded dihedral angle energy term;
components .lamda..sub.vdw and .lamda..sub.elec for the van de
Waals and electrostatic energy terms of the nonbonded exclusion and
1-4 pair interaction, respectively. The coupling vector .lamda. can
further include components .lamda..sub.other for other interactions
not affected by the formation and/or breaking of the valence bond
between the reference state and the target state. Each component of
the coupling vector .lamda. belongs to [0, 1]. During the course of
the alchemical transformation from the reference state to the
target state, the interactions unique in the reference state can be
turned off according to a first set of schedules for different
.lamda. components, and the interactions unique in the target state
can be turned on according to a second set of schedules for
different .lamda. components, as will be further described
below.
[0098] In popular molecular mechanics force fields, the bonded
angle and bonded dihedral angle interactions usually have the
following potential energy form:
U ba ( .theta. ) = k .theta. ( .theta. - .theta. 0 ) 2 ( 22 ) U bd
( .phi. ) = k .phi. n f ( cos ( n .phi. ) ) ( 23 ) ##EQU00027##
where .theta. is the bond angle, .theta..sub.0 is the equilibrium
bond angle, k.sub..theta. is the angle force constant (both
.theta..sub.0 and k.sub..theta. depend on the atoms forming the
bond angle); .phi. is the dihedral angle, k.sub..phi. is the
dihedral angle force constant (which depends on the atoms forming
the dihedrals). With the opening or closing of a ring, the bonded
angle and dihedral angle terms that are affected by the breaking or
forming of the bond can be modulated by components .lamda..sub.ba
and .lamda..sub.bd of the coupling parameter .lamda., respectively.
Although the bonded angle and dihedral interaction terms are
bounded for all .lamda. values, the absolute values of these terms
can be very large if the molecule is in a very twisted geometry. To
improve the accuracy of the free energy calculation, in some
embodiments, the bonded stretch interaction can be first turned on
to a significant degree (e.g., .lamda..sub.sbs=0.5) before turning
on the bonded angle and bonded dihedral angle interaction (during
bond formation). In this way, the bonded stretch interaction will
steer the molecule clear from a severely twisted geometry,
improving the inaccuracy problem caused by the bonded angle and
bonded dihedral angle interactions in the free energy
calculations.
[0099] The nonbonded electrostatic and van der Waals interactions
between two atoms are usually modeled by the following potential
energy form in popular molecular mechanics force fields:
U elec ( r ) = C q 1 q 2 r ( 24 ) U vdw ( r ) = 4 [ ( .sigma. r )
12 - ( .sigma. r ) 6 ] ( 25 ) ##EQU00028##
where in Eqs. 24 and 25, r is the inter-atom distance, q.sub.1 and
q.sub.2 are the charges of the two atoms, C is a constant,
.epsilon. is the depth of the potential well of U.sub.vdw(r), and
.sigma. is the finite distance at which the inter-atom potential
U.sub.vdw(r) is zero.
[0100] Many force fields, including OPLS, CHARM, AMBER, exclude or
modify the nonbonded interactions between atoms separated by one,
two, or three bonds. In particular, when two atoms are separated by
three bonds (e.g., the two atoms A.sub.1 and A.sub.8 shown in FIG.
1a), the conventional non-bonded electrostatic and van der Waals
interactions as described in Eq. 24 between the two atoms are not
calculated (i.e., excluded); instead, they are substituted by
electrostatic 1-4 pair and van der Waals 1-4 pair interactions.
[0101] With the opening or closing of a ring, the nonbonded
electrostatic and van der Waals interactions for atoms that span
the bond that is broken or forming can be modulated for by
components .lamda..sub.elec and .lamda..sub.vdw, respectively,
which will be further discussed below. As used herein, in some
embodiments of the invention, the van der Waals and/or the
electrostatic interaction can be described by a soft-core Lennard
Jones (LJ) and/or soft-core Coulomb potential that keeps pairwise
interaction energies finite for all configurations. By way of
example, a soft-core LJ potential may take the following functional
form:
U LJ ( r , .lamda. ) = 4 .lamda. ( 1 ( .alpha. LJ ( 1 - .lamda. ) n
+ ( r / .sigma. ) 6 ) 2 - 1 .alpha. LJ ( 1 - .lamda. ) n + ( r /
.sigma. ) 6 ) ( 26 ) ##EQU00029##
where it .alpha..sub.LJ is the soft core parameter, .epsilon. and
.sigma. are the standard Lennard Jones interaction parameters, n is
the order parameter which usually takes values between 1 and 6. See
Beutler et al., "Avoiding singularities and numerical instabilities
in free energy calculations based on molecular simulations," Chem.
Phys. Lett., 222 (1994) 529-539, the disclosure of which is
incorporated by reference herein. The soft-core LJ interaction
recovers the standard Lennard Jones interactions when .lamda.=1 and
it becomes 0 when .lamda.=0.
[0102] In some embodiments, for ring opening and closing free
energy calculations, a .lamda. schedule for all the system states
can be used to treat the bonded stretch, bonded angle, bonded
dihedral angle, nonbonded exclusion and 1-4 pair interactions, that
are affected by the formation and/or annihilation of a bond that
results in the opening or the closing of a ring structure. The
coupling parameter .lamda. for those interactions affected by the
bond formation or breaking may include seven components for bond
breaking transformation and seven components for bond formation
transformation from the reference state to the target state.
[0103] The seven components of the coupling parameter .lamda.
applicable for the bond breaking transformation include the
following terms: .lamda..sub.sbsA, which modulates the bonded
stretch interactions, respectively, that are present in the initial
state (the bond is present and yet to be broken) but not in the
final state (the bond is broken); .lamda..sub.baA and
.lamda..sub.bdA, which modulate the bonded angle and bonded
dihedral angle interactions that are present in the initial state
but not in the final state; .lamda..sub.elecA.sub.ex and
.lamda..sub.vdwA.sub.ex, which modulate the nonbonded electrostatic
interactions and van der Waals interactions present in the final
state but excluded (nonbonded exclusions) in the initial state due
to the breaking of a bond present in the initial state during the
transformation from the initial state to the final state, and
.lamda..sub.elecA.sub.14 and .lamda..sub.vdwA.sub.14 which modulate
the electrostatic 1-4 pair and van der Waals 1-4 pair interactions
present in the initial state but not included in the final state
due to the breaking or formation of the bond during the
transformation from the initial state to the final state.
[0104] Similarly, the seven components for bond formation
transformation include the following terms: .lamda..sub.sbsB, which
modulates the bonded stretch interactions that are present in the
final state (the bond is formed) but not in the initial state (the
bond is yet to be formed), .lamda..sub.baB and .lamda..sub.bdB,
which modulate the bonded angle and bonded dihedral angle
interactions that are present in the final state but not in the
initial state, respectively, .lamda..sub.elecB.sub.ex and
.lamda..sub.vdwB.sub.ex which modulate the nonbonded electrostatic
and van der Waals interactions present in the initial state but
excluded in the final state due to the formation of a bond in the
final state, and A.sub.elecB.sub.14 and A.sub.vdwB.sub.14, which
modulate the electrostatic 1-4 pair and van der Waals 1-4 pair
interactions present in the final state but excluded in the initial
state due to the formation or breaking of the bond during the
transformation from the initial state to the final state.
[0105] As used herein, the term "modulate" when used in connection
with a component of the coupling parameter .lamda. means that the
interaction energy for that particular component is calculated
using parameters of a conventional force field and the
corresponding component coupling parameter. To be specific, for
bonded stretch interaction and the LJ interaction, (or the
electrostatic interaction) where the soft-core potentials are used,
the interaction energies are calculated according to equations 9
and 26 respectively, while for other types of the interactions, the
interaction energy is calculated using parameters of a conventional
force field multiplied by the corresponding component coupling
parameter .lamda..
[0106] One example of .lamda. schedules discussed above is shown in
Scheme 1 below, in which the superscript (0, 1, . . . , m, m+1, . .
. , n) indicate the indexes of the reference system state, the
intermediate states, and the target system state, and
.lamda..sub.comp.sup.i is the respective .lamda. component value
for state with index "i" and component "comp".
.lamda. sbsA = [ .lamda. sbsA 0 = 1 , .lamda. sbsA 1 .lamda. sbsA m
= 0.5 , .lamda. sbsA m + 1 .lamda. sbsA n = 0 ] .lamda. baA = [
.lamda. baA 0 = 1 , .lamda. baA 1 .lamda. baA m = 0 , .lamda. baA m
+ 1 = 0 .lamda. baA n = 0 ] .lamda. bdA = [ .lamda. bdA 0 = 1 ,
.lamda. baA 1 .lamda. baA m = 0 , .lamda. baA m + 1 = 0 .lamda. baA
n = 0 ] .lamda. elecA ex = [ .lamda. elecA ex 0 = 0 , .lamda. elecA
ex 1 = 0 .lamda. elecA ex m = 0 , .lamda. elecA ex m + 1 .lamda.
elecA ex n = 1 ] .lamda. vdwA ex = [ .lamda. vdwA ex 0 = 0 ,
.lamda. vdwA ex 1 .lamda. vdwA ex m = 1 , .lamda. vdwA ex m + 1 = 1
.lamda. vdwA ex n = 1 ] .lamda. elecA 14 = [ .lamda. elecA 14 0 = 1
, .lamda. elecA 14 1 .lamda. elecA 14 m = 0 , .lamda. elecA 14 m +
1 = 0 .lamda. elecA 14 n = 0 ] .lamda. vdwA 14 = [ .lamda. vdwA 14
0 = 1 , .lamda. vdwA 14 1 .lamda. vdwA 14 m = 1 , .lamda. vdwA 14 m
+ 1 .lamda. vdwA 14 n = 0 ] .lamda. sbsB = [ .lamda. sbsB 0 = 0 ,
.lamda. sbsB 1 .lamda. sbsB m = 0.5 , .lamda. sbsB m + 1 .lamda.
sbsB n = 1 ] .lamda. baB = [ .lamda. baB 0 = 0 , .lamda. baB 1 = 0
.lamda. baB m = 0 , .lamda. baB m + 1 .lamda. baB n = 1 ] .lamda.
bdB = [ .lamda. bdB 0 = 0 , .lamda. bdB 1 = 0 .lamda. bdB m = 0 ,
.lamda. bdB m + 1 .lamda. bdB n = 1 ] .lamda. elecB ex = [ .lamda.
elecB ex 0 = 1 , .lamda. elecB ex 1 .lamda. elecB ex m = 0 ,
.lamda. elecB ex m + 1 .lamda. elecB ex n = 0 ] .lamda. vdwB ex = [
.lamda. vdwB ex 0 = 1 , .lamda. vdwB ex 1 = 1 .lamda. vdwB ex m = 1
, .lamda. vdwB ex m + 1 .lamda. vdwB ex n = 0 ] .lamda. elecB 14 =
[ .lamda. elecB 14 0 = 0 , .lamda. elecB 14 1 = 0 .lamda. elecB 14
m = 0 , .lamda. elecB 14 m + 1 .lamda. elecB 14 n = 1 ] .lamda.
vdwB 14 = [ .lamda. vdwB 14 0 = 0 , .lamda. vdwB ex 1 .lamda. vdwB
ex m = 1 , .lamda. vdwB ex m + 1 = 1 .lamda. vdwB ex n = 1 ] Scheme
1 ##EQU00030##
[0107] .lamda..sub.sbsA as shown in Scheme 1 can be varied from 1
to 0 over the bond breaking (e.g., a ring opening) transformation.
In some embodiments, the variation of .lamda..sub.sbsA can be
linear and/or monotonic over the transformation. In other
embodiments, the variation of .lamda..sub.sbsA can be non-linear
and/or non-monotonic over the transformation. Although it is shown
that .lamda..sub.sbsA=0.5 at the intermediate system state indexed
by m, this schedule for .lamda..sub.sbsA is merely illustrative and
non-limiting (e.g., other values smaller or greater than 0.5 can
also be used for the intermediate system state indexed by m).
[0108] .lamda..sub.baA as shown in Scheme 1 can be varied from 1 to
0 over the bond breaking (e.g., a ring opening) transformation. In
some embodiments, the variation of .lamda..sub.baA can be linear
and/or monotonic over the transformation or a portion of the
transformation. In other embodiments, the variation of
.lamda..sub.baA can be non-linear and/or non-monotonic over the
transformation or a portion of the transformation.
[0109] .lamda..sub.bdA as shown in Scheme 1 can be varied from 1 to
0 over the bond breaking (e.g., a ring opening) transformation. In
some embodiments, the variation of .lamda..sub.bdA can be linear
and/or monotonic over the transformation or a portion of the
transformation. In other embodiments, the variation of
.lamda..sub.bdA can be non-linear and/or non-monotonic over the
transformation or a portion of the transformation.
[0110] Further, in certain embodiments, for improved sampling
efficiency in the molecular simulations, the bonded angle and
bonded dihedral interactions between two atoms can be turned off
more quickly to 0 before the bonded stretch interactions are turned
off during the bond breaking transformation (conversely, the bonded
angle and bonded dihedral interactions can be turned on only after
the bonded stretch interactions are turned on to a predetermined
degree). Although it is shown that .lamda..sub.baA/.lamda..sub.bdA
can be decreased to 0 at the intermediate system state indexed by m
(meaning that .lamda..sub.baA/.lamda..sub.bdA can be varied from 1
at the initial state to 0 at this intermediate state by a more
rapid decrease than that of .lamda..sub.sbsA), the schedule for
.lamda..sub.baA/.lamda..sub.bdA are merely illustrative and
non-limiting (e.g. .lamda..sub.baA/.lamda..sub.bdA can be decreased
to 0 more rapidly or slowly from the initial state). Also, the
.lamda..sub.baA and .lamda..sub.bdA can be varied separately
according to their own respective schedules and do not need to be
synchronized.
[0111] .lamda..sub.elecA.sub.ex as shown in Scheme 1 can be varied
from 0 to 1 over the bond breaking (e.g., a ring opening)
transformation. In some embodiments, the variation of
.lamda..sub.elecA.sub.ex can be linear and/or monotonic over the
transformation or a portion of the transformation. In other
embodiments, the variation of .lamda..sub.elecA.sub.ex can be
non-linear and/or non-monotonic over the transformation or a
portion of the transformation. Although it is shown that
.lamda..sub.elecA.sub.ex=0 at the intermediate system state indexed
by m (i.e., .lamda..sub.elecA.sub.ex can be 0 from the initial
state through the intermediate state indexed by m, and increased
from this intermediate state to the final state), this schedule is
merely illustrative and non-limiting (e.g., can be kept at 0 until
an intermediate state that precedes or subsequent to the
intermediate state indexed by m).
[0112] .lamda..sub.vdwA.sub.ex as shown in Scheme 1 can be varied
from 0 to 1 over the bond breaking (e.g., a ring opening)
transformation. In some embodiments, the variation of
.lamda..sub.vdwA.sub.ex can be linear and/or monotonic over the
transformation or a portion of the transformation. In other
embodiments, the variation of .lamda..sub.vdwA.sub.ex can be
non-linear and/or non-monotonic over the transformation or a
portion of the transformation. Although it is shown that
.lamda..sub.vdwA.sub.ex=1 at the intermediate system state indexed
by m (.lamda..sub.vdwA.sub.ex can be increased from the initial
state at the value of 0 to the intermediate state indexed by m at
the value of 1, and stay at 1 thereon until the final state), this
schedule is merely illustrative and non-limiting (e.g.,
.lamda..sub.vdwA.sub.ex can be increased to 1 at an intermediate
state that precedes or follows the intermediate state indexed by
m). Further, in certain embodiments, for improved sampling
efficiency in the molecular simulations, the van der Waals
interaction between two atoms can be fully turned on before the
electrostatic interactions are turned on during the transformation.
This is illustrated in Schedule 1, where the schedules of
A.sub.elecA.sub.ex and .lamda..sub.vdwA.sub.ex are coordinated such
that when .lamda..sub.vdwA.sub.ex is smaller than 1 for an
intermediate state, .lamda..sub.elecA.sub.ex is 0 for that
intermediate state.
[0113] .lamda..sub.elecA.sub.14 as shown in Scheme 1 can be varied
from 1 to 0 over the bond breaking (e.g., a ring opening)
transformation. In some embodiments, the variation of
.lamda..sub.elecA.sub.14 can be linear and/or monotonic over the
transformation or a portion of the transformation. In other
embodiments, the variation of .lamda..sub.elecA.sub.14 can be
non-linear and/or non-monotonic over the transformation or a
portion of the transformation. Although it is shown that
A.sub.elecA.sub.14=0 at the intermediate system state indexed by m
(i.e., .lamda..sub.elecA.sub.14 can be varied from 1 at the initial
state to 0 at this intermediate state by a more rapid decrease than
that of .lamda..sub.sbsA), this schedule is merely illustrative and
non-limiting (e.g., .lamda..sub.elecA.sub.14 can be decreased to 0
more rapidly or slowly from the initial state).
[0114] .lamda..sub.vdwA.sub.14 as shown in Scheme 1 can be varied
from 1 to 0 over the bond breaking (e.g., a ring opening)
transformation. In some embodiments, the variation of
.lamda..sub.vdwA.sub.14 can be linear and/or monotonic over the
transformation or a portion of the transformation. In other
embodiments, the variation of .lamda..sub.vdwA.sub.14 can be
non-linear and/or non-monotonic over the transformation or a
portion of the transformation. Although it is shown that
.lamda..sub.vdwA.sub.14=1 at the intermediate system state indexed
by m (.lamda..sub.vdwA.sub.14 can be kept at 1 from the initial
state through the intermediate state indexed by m, and then
decreased to 0 from that intermediate state to the final state),
this schedule is merely illustrative and non-limiting (e.g.,
.lamda..sub.vdwA.sub.14 can be kept to 1 from the initial state
through an intermediate state that precedes or follows the
intermediate state indexed by m, and then decreased to 0 from that
intermediate state to the final state). Further, in certain
embodiments, for improved sampling efficiency in the molecular
simulations, the van der Waals interactions between two atoms are
turned off only after the electrostatic interactions are fully
turned off during the transformation. This is illustrated in
Schedule 1, where the schedules of .lamda..sub.elecA.sub.14 and
.lamda..sub.vdwA.sub.14 are coordinated such that when
.lamda..sub.vdwA.sub.14 is smaller than 1 for an intermediate
state, .lamda..sub.elecA.sub.14 is 0 for that intermediate
state.
[0115] The schedules of the .lamda. components for the bond
formation transformation, .lamda..sub.sbsB, .lamda..sub.baB and
.lamda..sub.bdB, .lamda..sub.elecB.sub.ex and
.lamda..sub.vdwB.sub.ex, .lamda..sub.elecB.sub.14 and
A.sub.vdwB.sub.14 can be readily understood by a person of ordinary
skill in the art in view of the above description of the .lamda.
schedule for the bond breaking transformation, as the two
transformation processes are essentially reverse to each other. It
is noted, however, that although the values of a .lamda. component
in the bond formation and the corresponding .lamda. component in
the bond breaking transformation (e.g., .lamda..sub.sbsA and
.lamda..sub.sbsB) appear to sum up to 1, this is merely
illustrative. For any given system, a bond cannot be both forming
and breaking simultaneously in the same alchemical transformation
considered in this application. Thus, the .lamda. component
schedules in a bond breaking transformation do not need to
coordinate with those in bond formation transformation. As an
example, .lamda..sub.sbsB can follow a schedule of increasing from
0 to 1 in a non-linear manner while .lamda..sub.sbsA follows a
schedule of decreasing from 1 to 0 in a linear manner. Further, as
shown in FIG. 1c, this application contemplates the breaking of one
bond (e.g., the bond between A.sub.1 and A.sub.4) and the
simultaneous formation of another bond (or bonds) (e.g., the bond
between A.sub.1 and A.sub.13 and the bond between A.sub.4 and
A.sub.13). In this case, the .lamda. component schedules for the
bond breaking transformation and the .lamda. component schedules
for the bond formation transformation can be selected separately
and not dependent upon each other.
[0116] In the following, energy calculations in a molecular
simulation of the transformation shown in FIG. 1a are illustrated
in connection with certain embodiments of the above-described
.lamda. schedule. [0117] Bonded stretch interactions that are
affected by the bond formation between A.sub.1 and A.sub.3 in FIG.
1a: (A.sub.1,A.sub.3) [0118] Bonded angle interactions that are
affected by the bond formation between A.sub.1 and A.sub.3 in FIG.
1a: [0119] (i) Bonded angle interactions only appearing in the
initial state and not in the final state:
(A.sub.2,A.sub.1,A.sub.10), (A.sub.4,A.sub.1,A.sub.10),
(A.sub.5,A.sub.1,A.sub.10), (A.sub.2,A.sub.3,A.sub.11),
(A.sub.8,A.sub.3,A.sub.11), (A.sub.9,A.sub.3,A.sub.11); [0120] (ii)
Bonded angle interactions only appearing in the final state and not
in the initial state: (A.sub.2,A.sub.1,A.sub.3),
(A.sub.4,A.sub.1,A.sub.3), (A.sub.5,A.sub.1,A.sub.3),
(A.sub.1,A.sub.3,A.sub.2), (A.sub.1,A.sub.3,A.sub.8),
(A.sub.1,A.sub.3,A.sub.9). [0121] Bonded dihedral interactions that
are affected by the bond formation between A.sub.1 and A.sub.3 in
FIG. 1a: [0122] (i) Bonded dihedral angle interactions involving
atoms only appearing in the initial state and not in the final
state: (A.sub.10,A.sub.1,A.sub.2,A.sub.3),
(A.sub.10,A.sub.1,A.sub.2,A.sub.6),
(A.sub.10,A.sub.1,A.sub.2,A.sub.7),
(A.sub.1,A.sub.2,A.sub.3,A.sub.11),
(A.sub.6,A.sub.2,A.sub.3,A.sub.11),
(A.sub.7,A.sub.2,A.sub.3,A.sub.11); [0123] (ii) Bonded dihedral
angle interactions involving atoms only appearing in the final
state and not in the initial state:
(A.sub.3,A.sub.1,A.sub.2,A.sub.6),
(A.sub.3,A.sub.1,A.sub.2,A.sub.7),
(A.sub.2,A.sub.1,A.sub.3,A.sub.8),
(A.sub.2,A.sub.1,A.sub.3,A.sub.9),
(A.sub.4,A.sub.1,A.sub.3,A.sub.2),
(A.sub.4,A.sub.1,A.sub.3,A.sub.8),
(A.sub.4,A.sub.1,A.sub.3,A.sub.9),
(A.sub.5,A.sub.1,A.sub.3,A.sub.2),
(A.sub.5,A.sub.1,A.sub.3,A.sub.8),
(A.sub.5,A.sub.1,A.sub.3,A.sub.9),
(A.sub.6,A.sub.2,A.sub.3,A.sub.1),
(A.sub.7,A.sub.2,A.sub.3,A.sub.1). [0124] Atom pairs excluded for
the calculation of nonbonded interactions in the initial state but
not excluded in the final state in FIG. 1a: (A.sub.1,A.sub.10),
(A.sub.1,A.sub.11), (A.sub.2,A.sub.10), (A.sub.2,A.sub.11),
(A.sub.3,A.sub.10), (A.sub.3,A.sub.11), (A.sub.4,A.sub.10),
(A.sub.5,A.sub.10), (A.sub.6,A.sub.10), (A.sub.6,A.sub.11),
(A.sub.7,A.sub.10), (A.sub.7,A.sub.11), (A.sub.8,A.sub.11),
(A.sub.9,A.sub.11). [0125] Atom pairs excluded for the calculation
of nonbonded interactions in the final state but not excluded in
the initial state in FIG. 1a: (A.sub.4,A.sub.8), (A.sub.4,A.sub.9),
(A.sub.5,A.sub.8), (A.sub.5,A.sub.9). [0126] 1-4 atom pairs
included in the initial state but not included in the final state
in FIG. 1a: (A.sub.1,A.sub.8), (A.sub.1,A.sub.9),
(A.sub.3,A.sub.4), (A.sub.3,A.sub.5). [0127] 1-4 atom pairs
included in the final state but not included in the initial state
in FIG. 1a: (A.sub.4,A.sub.8), (A.sub.4,A.sub.9),
(A.sub.5,A.sub.8), (A.sub.5,A.sub.9).
[0128] The other bonded and nonbonded interactions that are not
affected by the formation and annihilation of bonds for the ring
opening or closing transformation can be modulated by a regular
.lamda. schedule as in conventional free energy perturbations not
involving ring opening and closing (e.g., incrementing
.lamda..sub.other from the initial state to the final state, where
the interactions between P.sub.AB and P.sub.A is scaled by
1-.lamda..sub.other, and the interactions between P.sub.AB and
P.sub.B is scaled by .lamda..sub.other). Further, all the
interactions involving atoms that only appear in the initial state
but missing from the final state (i.e., those atoms that become
"dummy" atoms in the final state) will also be treated by a normal
.lamda. schedule. Conversely, for the reverse (bond formation)
transformation, all the interactions involving atoms that appear
only in the final state and are missing from the initial state
(e.g., A.sub.11 and A.sub.12 in FIG. 1a) will also be treated by
the normal .lamda. schedule. The remaining interactions, i.e.,
those affected by the formation and annihilation of bonds for the
ring opening or closing transformation that do not involve atoms
that only appear in the initial state or only appear in the final
state can be modulated by a "special" .lamda. schedule similar to
what is described in Schedule 1, that is, for the ring closing
transformation depicted in FIG. 1a: [0129] Bonded stretch
interaction included in the initial state but not in the final
state is modulated by .lamda..sub.sbsA using a soft bond
interaction potential described herein (in this example, no
interaction belongs to this category); [0130] Bonded stretch
interaction not included in the initial state but included in the
final state is modulated by .lamda..sub.sbsB using a soft bond
interaction potential described herein (in this example, the bonded
stretch interaction between atoms A.sub.1 and A.sub.3 belong to
this category); [0131] Bonded angle interactions included in the
initial state but not in the final state are modulated by
.lamda..sub.baA (in this example, no interaction belongs to this
category); [0132] Bonded angle interactions included in the final
state but not in the initial state are modulated by .lamda..sub.baB
(in this example, the bonded angle interactions between atoms
(A.sub.2,A.sub.1,A.sub.3), (A.sub.4,A.sub.1,A.sub.3),
(A.sub.5,A.sub.1,A.sub.3), (A.sub.1,A.sub.3,A.sub.2),
(A.sub.1,A.sub.3,A.sub.8), (A.sub.1,A.sub.3,A.sub.9) belong to this
category); [0133] Bonded dihedral angle interactions included in
the initial state but not in the final state are modulated by
.lamda..sub.bdA (in this example, no interaction belongs to this
category); [0134] Bonded dihedral angle interactions included in
the final state but not in the initial state are modulated by
.lamda..sub.bdB (in this example, the bonded dihedral angle
interactions between atoms (A.sub.3,A.sub.1,A.sub.2,A.sub.6),
(A.sub.3,A.sub.1,A.sub.2,A.sub.7),
(A.sub.2,A.sub.1,A.sub.3,A.sub.8),
(A.sub.2,A.sub.1,A.sub.3,A.sub.9),
(A.sub.4,A.sub.1,A.sub.3,A.sub.2),
(A.sub.4,A.sub.1,A.sub.3,A.sub.8),
(A.sub.4,A.sub.1,A.sub.3,A.sub.9),
(A.sub.5,A.sub.1,A.sub.3,A.sub.2),
(A.sub.5,A.sub.1,A.sub.3,A.sub.8),
(A.sub.5,A.sub.1,A.sub.3,A.sub.9),
(A.sub.6,A.sub.2,A.sub.3,A.sub.1),
(A.sub.7,A.sub.2,A.sub.3,A.sub.1) belong to this category); [0135]
Interactions excluded in the initial state but not excluded in the
final state are modulated by .lamda..sub.elecAex and
.lamda..sub.vdwAex (in this example, no interaction belongs to this
category); [0136] Interactions excluded in the final state but not
excluded in the initial state are modulated by .lamda..sub.elecBex
and .lamda..sub.vdwBex (in this example, the interactions excluded
in the final state for pairs (A.sub.4,A.sub.8), (A.sub.4,A.sub.9),
(A.sub.5,A.sub.8), (A.sub.5,A.sub.9) belong to this category);
[0137] 1-4 pair interactions included in the initial state but not
in the final state are modulated by .lamda..sub.elecA14 and
.lamda..sub.vdwA14 (in this example, the 1-4 pair interactions for
pairs (A.sub.1,A.sub.8), (A.sub.1,A.sub.9), (A.sub.3,A.sub.4),
(A.sub.3,A.sub.5) belong to this category); and [0138] 1-4 pair
interactions included in the final state but not in the initial
state are modulated by .lamda..sub.elecB14 and .lamda..sub.vdwB14
(in this example, the 1-4 pair interactions for pairs
(A.sub.4,A.sub.8), (A.sub.4,A.sub.9), (A.sub.5,A.sub.8),
(A.sub.5,A.sub.9) belong to this category).
[0139] With the energy terms defined by a suitable .lamda. schedule
(such as the one depicted in Scheme 1) for all the system states in
the transformation from the initial state to the final state,
molecular simulations can be run to sample the ensembles of
micro-states obtained at the reference state, the target state, and
the intermediate states according to the .lamda. schedule. For each
.lamda. window, the free energy difference can be calculated
between all the neighboring lambda windows
.DELTA.F.sub..lamda.i.fwdarw..lamda.i+1 and/or between any pair of
lambda windows .DELTA.F.sub.i.fwdarw.j, including between the
reference state and the target state. The total free energy
difference between the reference state and the target state can be
obtained by adding the free energy differences between each two
successive state along the transformation path defined by the
.lamda. schedule or directly obtained by analyzing the data from
all the sampled states.
[0140] The free energy difference between neighboring .lamda.
windows or generally between any pair of states including initial
state and the final state can be calculated by a variety of ways.
For example, e.g., by the use of internal energy difference (FEP
NVT ensemble), the enthalpy difference (FEP NPT ensemble), or other
related thermodynamic property difference (FEP other ensembles,
such as the NVE ensemble), of a suitable ensemble of the
micro-states obtained at the target state, the reference state, and
the intermediate states as coupling parameter is instantaneously
varied for the selected ensemble of micro-states. The analysis can
be further performed, for example, by way of Bannet Acceptance
Ratio (BAR), Multistate Bannet Acceptance Ratio (MBAR), Weighted
Histogram Analysis Method (WHAM), Zwanzig averaging, or other
similar FEP-family estimators. Alternatively, the free energy
difference between neighboring windows or generally between any
pairs of states including the initial state and the final state can
be calculated by way of an analysis the derivative of the energy
with respect of the coupling vector .lamda. (TI NVT ensemble), the
derivative of the enthalpy with respect to the coupling vector
.lamda. (TI NPT ensemble), or the derivative of other related
thermodynamic properties with respect to the coupling vector
.lamda. (TI other ensembles, such as the TI NVE ensemble), for each
microscopic state obtained. In other embodiments, the free energy
difference of each .lamda. window can be calculated by way of an
analysis of the potential of mean force (PMF) associated with
sampling of the coupling vector .lamda. as a dynamical variable
that can dynamically transition between the reference state, the
target state and intermediate states for example and without loss
of generality via the .lamda.-dynamics, the principle of which is
generally described in Knight et al., .lamda.-dynamics free energy
simulation methods, J. Comput. Chem., 2009, 30: 1692-1700, the
disclosure of which is incorporated by reference herein.
.lamda.-dynamics based sampling methods include, but are not
limited to, .lamda.-Monte Carlo, .lamda.-metadynamics,
.lamda.-OSRW, and other .lamda. PMF sampling family methods.
[0141] The methods for free energy difference calculations
described herein can be applied to a number of highly useful
applications, which include, for example: [0142] Relative
protein-ligand binding affinity and/or relative solvation free
energy calculations between congeneric ligands with ring opening or
closing (see FIG. 4); [0143] Relative protein-ligand binding
affinity and/or relative solvation free energy calculations between
congeneric ligands that differ by a macrocyclization (see FIG. 5);
[0144] The calculation of the effect of a non-proline to proline or
proline to non-proline residue mutation to protein thermodynamic
stability, protein-ligand binding affinity, or protein-protein
binding affinity (see FIG. 6); and [0145] The calculation of the
effect of a residue insertion or residue deletion to protein
thermodynamic stability, protein-ligand binding affinity, or
protein-protein binding affinity (see FIG. 7, which schematically
shows a transformation of the illustrative protein segment
structure depicted on the left by losing a residue numbered 129,
resulting in the structure depicted on the right).
[0146] Embodiments of the method for the free energy calculations
of the disclosed subject matter can be implemented in a computer
program, which can take the form of a software component of a
suitable hardware platform, for example, a standalone computer, a
networked computer, a network server computer, a handheld device,
or the like. Different aspects of the disclosed methods may be
implemented in different software modules and executed by one
processor or different processors, sequentially or in parallel,
depending on how the software is designed. The apparatus on which
the program can be executed can include one or more processors, one
or more memory devices (such as ROM, RAM, flash memory, hard drive,
optical drive, etc.), input/output devices, network interfaces, and
other peripheral devices. A computer readable non-transitory media
storing the program is also provided.
[0147] The disclosed subject matter is not to be limited in scope
by the specific embodiments described herein. Indeed, various
modifications of the disclosed subject matter in addition to those
described herein will become apparent to those skilled in the art
from the foregoing description and the accompanying figures. Such
modifications therefore fall within the scope of the appended
claims.
* * * * *