Calculation of Absolute Protein-Ligand Binding Affinity Using Path and Endpoint Approaches

15 Pages·2006·0.19 MB·English

864 BiophysicalJournal Volume90 February2006 864–877
Calculation of Absolute Protein-Ligand Binding Affinity Using Path
and Endpoint Approaches
Michael S. Lee*y and Mark A. Olson*
*DepartmentofCellBiologyandBiochemistry,U.S.ArmyMedicalResearchInstituteofInfectiousDiseases,Frederick,Maryland21702;
andyComputationalandInformationSciencesDirectorate,U.S.ArmyResearchLaboratory,AberdeenProvingGround,Maryland21005
ABSTRACT Acomparativeanalysisisprovidedofrigorousandapproximatemethodsforcalculatingabsolutebindingaffinities
oftwoprotein-ligandcomplexes:theFKBPproteinboundwithsmallmolecules4-hydroxy-2-butanoneandFK506.Ourrigorous
approachisanumbrellasamplingtechniquewhereapotentialofmeanforceisdeterminedbypullingtheligandoutoftheprotein
activesiteoverseveralsimulationwindows.Theresultsofthisapproachagreewellwithexperimentallyobservedbindingaffinities.
Alsoassessedisacommonlyusedapproximateendpointapproach,whichseparatelyestimatesenthalpy,solvationfreeenergy,
andentropy.Weshowthatthisendpointapproachhasnumerousvariations,allofwhicharepronetocriticalshortcomings.For
example,conventionalharmonicandquasiharmonicentropyestimationproceduresproducedisparateresultsfortherelatively
simpleprotein-ligandsystemsstudiedinthiswork.
INTRODUCTION
The accurate calculation of absolute binding affinities of ofliquids.Theexponentialimprovementsincomputerhard-
protein-ligand complexes is an important goal for the study ware as well as enhanced molecular dynamics algorithms
ofbiomolecularrecognition(1)andcomputationaldrugdesign makethePMFapproacharealityforprotein-ligandsystems.
(2). However, currently available computational methods Nevertheless, the computational requirements are still quite
often require some knowledge of experimental binding demanding. Izrailev et al. (10) have been using pulling
affinitiestocalibrateparametersforaparticularproteintarget methods for over a decade to study the nature of molecular
(3). Free-energy techniques, known as double decoupling recognition in protein-protein complexes. Fukunishi et al.
methods (4–6), have been developed to calculate the abso- (11) devised an approach to estimate the free energy of
lute binding affinities of complexes without a priori exper- bindinginprotein-ligandcomplexesutilizingaself-avoiding
imental information. The methods involve calculating the randomwalkprocedure.Also,inthelastyear,WooandRoux
free-energycyclefordecouplingtheproteinandligand,and (12)successfullyappliedaPMFapproachtothecalculation
then reintroducing the ligand to the bulk solvent. This rig- of the equilibrium binding constant of the phosphotyrosine
oroustechniquehasonlybeenusedforverysmallligands(4) peptidepYEEItotheSrchomology2domainofhumanLck.
or with simplistic implicit solvent models (6). One of the A commonly used approximate method for the calcula-
difficultiesinvolved inthisapproach isthat theligandmust tion of absolute binding affinities is the so-called molecular
be decoupled slowly enough from the binding pocket such mechanics-Poisson-Boltzmann-surface area (MM/PB-SA)
thatthemechanicalworkassociatedwiththeprocesscanbe method (13,14). In this approach, an explicit solvent simu-
performedreversibly.Newtechniqueshavebeendeveloped lationoftheboundstateiscarriedout.Thenthesimulationis
that can obtain free energies from repeated nonequilibrium postprocessed to determine the enthalpic differences be-
simulations(7,8)andmayhelpmakedoubledecouplingappli- tween the bound and unbound solute states. The solvation
cationsmoreefficient.Usingadifferentstrategy,Changetal. free energy of binding is obtained from a Poisson-based
enumerated the configuration integrals of the bound and solvation model (15). Separately, the binding entropy is
unbound state of simple host-guest complexes to calculate estimated by harmonic analysis using a simple r-dielectric
thefree energy of association (9). function to approximate solvent screening of charge-charge
Any alchemical pathway between bound and unbound interactions.
states can, in principle, be used to obtain free energies of Thereareafewvariationstothismethodintheliterature
complexformation.Oneofthemostobviouspathwaysisto asseen in Table 1. One idea is to usethe fastergeneralized
simplypullouttheligandfromtheactivesiteoftheprotein Born(GB)solventmodeltoscorestructuresversusthePois-
by a potential of mean force (PMF) approach. The PMF sonmethod,alsoknownasMM/GB-SA(16).Also,thesimu-
approach has existed since the early days of molecular lations can be run with explicit solvent or GB. Another
mechanicsandiswellgrounded inthestatistical mechanics method requires running both bound and unbound simu-
lations.Afurtherchoiceisthatthevibrationalentropycanbe
calculated using either the normal mode approach or the
quasiharmonicapproximation.Moreover,thequasiharmonic
SubmittedJuly28,2005,andacceptedforpublicationOctober20,2005.
entropy calculation can be performed in two ways: one
Address reprint requests to Michael S. Lee, E-mail: michael.lee@amedd.
army.mil.
(cid:1)2006bytheBiophysicalSociety
0006-3495/06/02/864/14 $2.00 doi:10.1529/biophysj.105.071589
Report Documentation Page Form Approved
OMB No. 0704-0188
Public reporting burden for the collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and
maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information,
including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington
VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to a penalty for failing to comply with a collection of information if it

does not display a currently valid OMB control number.
1. REPORT DATE 2. REPORT TYPE 3. DATES COVERED
1 FEB 2006 N/A -
4. TITLE AND SUBTITLE 5a. CONTRACT NUMBER
Calculation of absolute protein-ligand binding affinity using path and
5b. GRANT NUMBER
endpoint approaches. Biophysical Journal 90:864 - 877
5c. PROGRAM ELEMENT NUMBER
6. AUTHOR(S) 5d. PROJECT NUMBER
Lee, MS Olson, MA
5e. TASK NUMBER
5f. WORK UNIT NUMBER
7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION
United States Army Medical Institute of Infectious Diseases, Fort Detrick, REPORT NUMBER
RPP-04-355
MD
9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR’S ACRONYM(S)
11. SPONSOR/MONITOR’S REPORT
NUMBER(S)
12. DISTRIBUTION/AVAILABILITY STATEMENT
Approved for public release, distribution unlimited
13. SUPPLEMENTARY NOTES
14. ABSTRACT
A comparative analysis is provided of rigorous and approximate methods for calculating absolute binding
affinities of two protein-ligand complexes: the FKBP protein bound with small molecules BUQ and FK506.
Our rigorous approach is an umbrella sampling technique where a potential of mean force is determined
by pulling the ligand out of the protein active site over several simulation windows. The results of this
approach agree well with experimentally observed binding affinities. Also assessed is a commonly used
approximate endpoint approach, which separately estimates enthalpy, solvation free energy, and entropy.
We show that this endpoint approach has numerous variations, all of which are prone to critical
shortcomings. For example, conventional harmonic and quasi-harmonic entropy estimation procedures
produce disparate results for the relatively simple protein-ligand systems studied in this work.
15. SUBJECT TERMS
computer modeling, Poisson solvation models, explicit, implicit, methods
16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF 18. NUMBER 19a. NAME OF
ABSTRACT OF PAGES RESPONSIBLE PERSON

a. REPORT b. ABSTRACT c. THIS PAGE SAR 14
unclassified unclassified unclassified
Standard Form 298 (Rev. 8-98)
Prescribed by ANSI Std Z39-18
AbsoluteProtein-LigandBindingAffinity 865
TABLE1 VariationsoftechniquesfortheMM/PB(GB)-SA TABLE2 Glossaryofsomeoftheabbreviationsandterms
method(54) usedinthiswork
Decision Options Abbreviation/term Definition
1.Explicitsolvent BUQ 1-hydroxybutanone
Solvationmodelfor
2.GB-SA COM Center-of-mass
ensemblegeneration
3.r-dielectric COG Center-of-geometry
and/orharmonicentropy
cpx Complex
1.One-state,complexonly
FKBP FK506-bindingprotein
2.Two-state,complex/freeprotein
Numberofsimulations GBMV2 GeneralizedBornmolecularvolumemethod(34)
andligand
GB-SA GeneralizedBornplusasurface-area-based
1.Explicitsolventchargingfree nonpolarterm
energy(36)* Hybrid Hybridexplicit/implicitsolventmethodthathasafixed
Electrostaticsolvation
2.PoissonBoltzmann boundarybutauniformwaterlayersurrounding
freeenergy
3.GeneralizedBorn thesolute(35)
lig Ligand
1.Singlesurface-areatermof
NM Normalmodevibrationalentropyanalysis
cavitationandnonpolar
Nonpolarsolvation prot Protein
2.Surface-areatermofcavitation,
freeenergy QH Quasiharmonicvibrationalentropyanalysis
Born-liketermfornonpolar(53)*
QH1 Standardquasiharmonicvibrationalentropyanalysis
1.Harmonicapproximation1ideal QH2 Quasiharmonicvibrationalentropyanalysiswhere
gasR/Tcorrection ligandtranslationsandrotationsarecalculated
2.Quasiharmonic1idealgasR/T independently
Entropy
correction RBT Rigid-bodytranslationoftheligand
3.Quasiharmonicvibration1rigid RBR Rigid-bodyrotationoftheligand
bodyrotationandtranslation RD4 Distance-dependentdielectricwithaprefactorof4
SASA-1 Solvent-accessiblesurface-areamethodforestimating
*Optionwasnotevaluatedinthiswork.
nonpolarsolvationfreeenergy(34)
methodisanalogoustothenormalmodeentropycalculation,
which uses an ideal gas correction and the other involves ½P(cid:1)½L(cid:1)
separation of the ligand translational and rotational degrees Kd ¼ ½PL(cid:1): (2)
of freedom from theother coordinates (17).
Inthisstudy,wetestedseveralvariationsoftheMM/GB-
SA method on two protein-ligand systems and compared Experimentally, Kd is the molar concentration of ligand
them to the more rigorous PMF method and experimental necessary to make a 50:50 mixture of bound and unbound
results. The outline of the article is the following. We first protein.Thefreeenergyofbindinginthestandardstate(1M)
describe thegeneral theoryandcomputational methods that isdefinedintermsofthedissociationconstant,
wereusedinthisstudy.Wethenpresentourresultsfollowed
(cid:1) (cid:2)
by adiscussionand summary of conclusions. DG ¼(cid:2)RTln Kd ; (3)
bind ½1M(cid:1)
where R is the universal gas constant and T is the absolute
THEORY
temperature. A general expression for calculating DG
bind
Thetheoreticalderivationofexpressionsthatcanbeusedfor from a computer simulationis(5)
the calculation of absolute binding affinities of a protein-
ligand complex from computer simulations has been (cid:1)C0 Z Z (cid:2)
DG ¼(cid:2)RTln PL1S S 1P0ÆDV æ; (4)
expounded elsewhere (5,18). We will outline some of the bind 8p2Z Z U/B
P1S L1S
keypointspertaining totheissues that areaddressedinthis
work. For reference, Table 2 lists some of the terms and where C0 is 1 M (;1 mol/1660 A˚3), Z corresponds to the
acronyms used in this work. The chemical reaction of configuration integral over the subscripted coordinates,
a simple two-state binary protein-ligand complex formation and S corresponds to solvent degrees of freedom. The
inan aqueous environment (aq)is second term in Eq. 4, which corresponds to the work
associated with the volume difference between the bound
P 1L 4ðPLÞ ; (1)
aq aq aq andunboundstates,isnegligibleinwaterandcanbesafely
where the left-hand side corresponds to the unbound state omitted(19).NokineticenergycomponentispresentinEq.
(U), and the right-hand side corresponds to the bound state 4 because these terms cancel out in the classical thermody-
(B) of the protein, P, and the ligand, L. The dissociation namic limit (5). For implicit solvent simulations, Eq. 4
constantfor thisreaction is reduces to (5)
BiophysicalJournal90(3)864–877
866 LeeandOlson
(cid:1)C0 z (cid:2) Potential of mean force approach
DG ¼(cid:2)RTln PL ; (5)
bind 8p2z z Most rigorous free-energy methods (e.g., free-energy
P L
perturbation and thermodynamic integration) determine
wherezrepresentsconfigurationintegralsinthecasethatthe
the free-energy difference between two states of a system
solvent degrees of freedom are embedded in an external
through the use of an arbitrarily defined path. For absolute
potential actingonthe complex and unbound states.
binding-affinity calculations, the two states are the bound
It is convenient to separate out the six coordinates asso-
complex and the unbound protein and ligand. In this work,
ciatedwiththerelativepositionandorientationoftheligand
thepathwaybetweenthesetwoendpointswaschosenasthe
withrespecttotheprotein,namely,ðr;jÞ.Thedetermination
pullingoftheligandoutofthecomplexintothebulksolvent
ofDG viasimulationcanthenbeformulatedintermsofa
bind along a straight line (12). A curvilinear pathway would be
bindingzone,B,thatcoversarangeofthesesixcoordinates,
necessary for certain types of protein-ligand complexes,
(cid:3)C0 Z (cid:4) especially when the protein atoms are held rigid (20). We
DGbind ¼(cid:2)kBTln 8p2 e(cid:2)Wðr;jÞ=kBTdrdj ; (6) defineaunit-vectorpullingdirection,v,andacoordinate,l,
B
whichmeasurestheprogressalongthepathwayfrombound
where Wðr;jÞ is the work necessary to bring the center- tounbound.Assumingthecenter-of-massandorientationof
of-mass (COM) of the ligand from infinity to a point ðr;jÞ the protein is fixed, the center-of-geometry (COG) of the
relativetotheCOMandorientationoftheprotein,andkBis ligand can be pulled out of the binding groove with a
the Boltzmann constant. Equation 6 can be reduced further translation vector, lv,
by integrating over the orientational degrees of freedom of
theligandin thereference frame of theprotein, rCOG ¼r0 1lv; (9)
lig lig
(cid:3) Z (cid:4) wherer0 istheinitialboundposition.Inthiswork,weused
DG ¼(cid:2)k Tln C0 e(cid:2)WðrÞ=kBTdr : (7) lig
bind B umbrella sampling to visit overlapping regions along the
B
path. We then calculate the PMF along this reaction
Theworkfunction,WðrÞ,canbeobtainedfromtheprob-
coordinate by combining the probability distributions of
ability, rðrÞ, of finding a ligand in a given position relative
each sampling window. The protein is restrained to be in
tothe protein,
a standard frame of reference and the COG of the ligand is
WðrÞ¼(cid:2)k Tln½rðrÞ(cid:1): (8) restrainedforwindowjatspatialpoint,rj,alongthepathway
B
in Eq. 8, using a biasing potential, VðrÞ, ofthe form
j
A difficulty arises in how one should numerically define
1
the bound state. In principle, the bound state should mimic VðrÞ¼ kðrCOG(cid:2)rÞ2: (10)
j 2 j lig j
the experimental zone whereby a ‘‘signal’’ is detected cor-
responding to binding. However, this zone is theoretically Also,anunrestrainedcomplexsimulationisperformedto
unknown. This issue is most critical for weakly bound determine the bound state and to increase sampling of the
complexes (19) where the integrands of the configuration complex. The distributions in three-dimensional space of
integralsinEqs.6and7aresignificantoverarelativelylarge the separate simulation windows are merged by using the
range of space. One practical suggestion (19), called an weighted histogram method (WHAM) (21,22) to form a
‘‘exclusionzone’’,istodefinethebindingzoneasthestates singleunweighteddistribution,rðrÞ,thatspansthetwostates
of a ligand for which no second ligand can have a more ofinterest.Becauseweareonlysamplingoverasmallfrac-
favorableinteractionwiththeprotein(19).Forthiswork,we tion of the total translational space, we choose to integrate
willapproximatethebindingzoneastherangeofcoordinates rðrÞ over the two directions perpendicular to the pulling
for the ligand COM associated with a computer simulation coordinate, thus yielding rðlÞ and subsequently WðlÞ via
oftheunrestrainedcomplex.Thisdefinition,albeitimperfect, analogytoEq.8.Wealsomaketheassumptionthatthework
is the basis behind other absolute binding free-energy function is effectively quadratic in the binding groove. The
methods in the literature such as the double decoupling bindingfreeenergythencanbecalculatedfromWðlÞ,aswell
methods(4,6). asanequilibriumsimulationofthecomplexvia(17)
(cid:3) Z (cid:1)(cid:5) 1 1 1 (cid:6)(cid:7) (cid:2) (cid:4)
DG (cid:3)(cid:2)k Tln C0 exp (cid:2)w (cid:2) k x2(cid:2) k x2(cid:2) k x2 k T dx dx dx
bind B min 2 1 1 2 2 2 2 3 3 B 1 2 3
B
(cid:3) qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi(cid:4)
(cid:3)w (cid:2)k T+ln ð2pÞ3=2C0 ÆDx2æ ÆDx2æ ÆDx2æ
min B 1 c 2 c 3 c
c
(cid:3)w 1DGRBTðligÞ; (11)
min cpx
BiophysicalJournal90(3)864–877
AbsoluteProtein-LigandBindingAffinity 867
where k (cid:3)k T=ÆDx2æ based on the classical harmonic components for binding are determined by subtracting the
i B i
oscillatorapproximation(17),w isthediscreteminimum averagecomplexenergyfromtheaverageproteinandligand
min
valueofWðlÞwhentheasymptoticvalueofWðlÞatlargel energiesobtainedbyeffectivelypullingthespeciesinfinitely
is shifted to zero, and the summed term, DGRBTðligÞ, is far apart,
cpx
defined as the free energy of rigid-body translation (RBT)
DEXð1-stateÞ¼ÆEXæ (cid:2)½ÆEXðprotÞæ 1ÆEXðligÞæ (cid:1);
of the ligand in the complex compared to confinement in cpx cpx cpx
a 1660-A˚3 volume. The rigid-body translation term and the (15)
ligand center-of-geometry fluctuation eigenvectors, Dx, for
i where EXðprotÞ is evaluated over the complex trajectory
eachcluster,c,areobtainedusingaprescriptiondescribedin
with the ligand deleted, and EXðligÞ is evaluated over the
thenextsection(17).Weusedacombinationoftwocriteria
complex trajectory with the protein deleted. This method
to determine the location of the unbound endpoint of the
leads to significantly lower statistical uncertainty but com-
PMF. First, one can look for the PMF to flatten out over
pletelyneglectstheenthalpiesofproteinandligandreorgani-
several angstroms and assume this to be the asymptote.
zation proceeding from the bound to the unbound state,
Second, the nonbonded cutoffs provide an approximate
which, inmany cases, may notbe negligible.
estimate of the distance such that the ligand and protein no
Thegeneralapproachtonormalmode(NM)andstandard
longer directly interact with eachother.
quasiharmonic entropy (denoted here as QH1) methods
involves computation of the vibrational (vib) and ideal gas
Endpoint approaches
entropies inthe bound andunboundstates,
In an endpoint approach representing MM/GB-SA models,
DSNM=QH1¼SvibðcpxÞ(cid:2)SvibðprotÞ(cid:2)SvibðligÞ1DStrans1DSrot ;
the basic formula for the binding free energy involves the bind cpx prot lig free free
sum of the enthalpic, entropic, and free-energy differences (16)
between the boundand unbound states,
where S denotes absolute entropy, the simulation type
bind
DG ¼DUgas1DGsolv(cid:2)TDS; (12) is in subscript, the evaluated degrees of freedom are in
bind
parentheses, and the ideal gas entropy components, DStrans
where DS is the entropy difference between the two states and DSrot , are described below. Quasiharmonic entrofrpeey
either in the gas or aqueous phase, and DGsolv reflects the free
methods are based on either classical or quantum harmonic
solvation free-energy difference. The term DUgas is the
oscillator theory (18,24). Both quantum and classical forms
difference in gasphase enthalpies defined as
produceroughlythesameresultwhendifferencesinentropy
DUgas (cid:3)DEgas ¼DEinternal1DECoulomb1DEvdW; (13) are taken (results not shown). We employ the quantum
modelinthiswork.Allquasiharmonic methods involvethe
wherethetermsDEinternal,DECoulomb,andDEvdWcorrespond formationof a covariancefluctuation matrix, C(24),
to differences in internal energy, electrostatic energy, and
van derWaals(vdW) dispersion energy, respectively. C¼M1=2(cid:3) 1 N+confðx (cid:2)x Þ5ðx (cid:2)x Þ(cid:4)M1=2; (17)
In the two-state version of the MM/GB-SA model, N i 0 i 0
conf i¼1
simulations of the complex, and free protein and ligand are
wherex istheconformationofsnapshoti,x istheaverage
performed.ChangesinenergytermsofcomponentX(e.g.,X i 0
structure over the entire set of conformations, N is the
¼ Coulomb) are obtained by independently averaging over conf
numberofconformations,5denotestheouterproduct,and
bound and unbound simulations,
M is a diagonal matrix of the atomic masses (24). Each
DEXð2-stateÞ¼ÆEXæ (cid:2)ðÆEXæ 1ÆEXæ Þ; (14) conformationusedinEq.17isrotatedandtranslatedtobest
cpx prot lig
fittheaveragestructure.Finally,theCmatrixisdiagonalized
where the subscript ‘‘cpx’’ of the thermal average, Æ(cid:4)(cid:4)(cid:4)æ; toobtaineigenvalues,e.Theseeigenvaluesareconvertedto
i
denotesthecomplex,andsimilarlyforthefreeprotein(prot)
frequencies (24),
and free ligand (lig). The two-state model, although
rffiffiffiffiffiffiffiffi
admirableinitssimplicity,isactuallyverydifficulttopursue 1 k T
n ¼ B : (18)
inpractice, because of atleast twocritical issues. First,it is i 2p e
i
difficult to pinpoint the average energy of an ensemble of
NormalmodeentropymethodsdifferinthattheCmatrix
proteinorcomplexstructures,becausethesystemmayjump
is calculated as the mass-normalized second derivative
frommacrostatetomacrostateovertimeperiodslargerthan
matrix (24),
what is currently feasible to simulate by today’s computing
resources (23). Second, the statistical uncertainty is quite C¼M(cid:2)1=2ð=2EÞM(cid:2)1=2: (19)
largefortakingthedifferencebetweenaverageenthalpiesof
boundandunboundensembles(16,17). The C matrix is diagonalized to yield eigenvalues. The
Often researchers revert to an one-endpoint simulation frequencies are then calculated as the square root of the
methodwhereonlythecomplexissimulatedandtheenergy eigenvalues. The frequencies from either harmonic or
BiophysicalJournal90(3)864–877
868 LeeandOlson
quasiharmonic methods are substituted into the vibrational clustered sets of snapshots, indexed as c, of the COM of
entropy formula for a noninteracting collection of quantum the ligand. Clustering was performed in three-dimensional
harmonic oscillators(24), spaceusingastandardhierarchicalagglomerationalgorithm
(26). Clustering was used to ensure that quasiharmonic
TSvib ¼(cid:2)Gvib1Hvib
basins were found such that Eq. 24 is appropriate. We
¼+(cid:2)k Tlnð1(cid:2)e(cid:2)hni=kBTÞ1 hni ; (20) arbitrarily selected the hierarchical level so that there were
B 1(cid:2)e(cid:2)hni=kBT five clusters. Results varied little with larger numbers of
i
clusters.Therotationaltermiscalculatedthroughhistogram
where h is Planck’s constant. In harmonic entropy and
binning of the spherical-polar angles, ðf;c;uÞ, which are
standard quasiharmonic binding entropy calculations, the
necessary to orient each ligand conformation into the
ideal gas entropy terms for translation (trans) and rotation
standard reference frame(27),
(rot)areusedtoaccountforthesixdegreesoffreedomofthe
free ligand and free protein at 1 M (C0) concentration in DSRBRðligÞ¼k N+int rðf;c;u Þlnrðf;c;u Þ(cid:2)rfreelnrfree;
solution (25): cpx B i j m i j m ijm ijm
i;j;m¼1
DStransðbindÞ¼StransðcpxÞ(cid:2)StransðprotÞ(cid:2)StransðligÞ; (25)
free free free free
DSrot ðbindÞ¼Srot ðcpxÞ(cid:2)Srot ðprotÞ(cid:2)Srot ðligÞ; where
free free free free
where rfree¼ 1 ½cosðfÞ(cid:2)cosðf Þ(cid:1); (26)
5 "1(cid:1)2pm k T(cid:2)3=2# ijm 2Ni2nt j j11
StransðxÞ¼ k 1k ln x B ; (21)
free 2 B B C0 h2 such that each of the angles are divided into N equally
int
spaced intervals, and rfree is the probability of finding the
Srot ðxÞ¼3k 1k ln ðpIxxIyxIzxÞ1=2(cid:3)8p2kBT(cid:4)3=2!; (22) ligandinacertainbinasisjmumingitcanrotatefreely.Thereare
free 2 B B s h2 noprecisecriteriaforhowtochoosethenumberofintervals
x
to achieve the best entropy estimate. Based on empirical
and Ix, Ix, and Ixare the moments of inertia of species x
x y z observations,we chose Nint ¼ 15,but also looked atvalues
(x¼cpx,prot,orlig),s isthesymmetryfactorofx,which
x of10 and 20to calculate errors associated with binning.
issetto1forthiswork.Analternativemethodforcalculat-
Accurately calculating entropy is a nontrivial task. For
ing quasiharmonic entropy, dubbed here as QH2, involves
example,therecanbeasignificantsourceoferrorduetothe
a complete decoupling of the RBT and rigid-body rotations
above-statedproblemofinadequatesamplingofmacrostates.
(RBR) of theligandfrom thecomplex(17),
In addition, the standard methods used to evaluate entropy
DSQH2 ¼DSðprotÞ1DSðligÞ1DSRBTðligÞ1DSRBRðligÞ make sweeping assumptions of ideal harmonic behavior,
bind cpx cpx
which may be inadequate for flexible proteins (28). For
¼S ðprotÞ(cid:2)S ðprotÞ1S ðligÞ(cid:2)S ðligÞ
cpx prot cpx lig instance,thequasiharmonicentropymethodassumesthatthe
1DSRBTðligÞ1DSRBRðligÞ: (23) protein is vibrating in a 3N (cid:2) 6 dimensional harmonic
cpx cpx atom
well. A simple illustration of nonideal harmonic behavior
The benefit of this approach is that one can directly
would be a freely rotating side chain on the surface of a
estimate DSRBTðligÞ and DSRBRðligÞ. The drawback of this
cpx cpx protein. Also, harmonic and quasiharmonic methods assume
procedure is that coupled motions between the RBT and
thateachmodeisuncoupledtoallothermodessuchthatthe
RBRoftheligandandtheproteinareneglectedandthusthe
entropy of each mode is additive (28). This only seems to
binding entropy may be underestimated. Also, it should be
be a reasonable assumption when sampling a single local
noted that we use rigid-body entropy terms rather than free
minimum.
energy terms as in Swanson et al. (17) to avoid double-
counting the rigid-body enthalpy. The rigid-body trans-
lational entropy term, DSRBTðligÞ, is(17)
cpx METHODS
DSRBTðligÞ¼½HRBTðligÞ(cid:2)DGRBTðligÞ(cid:1)=T
cpx cpx cpx The two protein-ligand complexes that were evaluated in this work were
3 FK506-binding protein/4-hydroxy-2-butanone (FKBP-BUQ) (Protein Data
¼2kB1kBln½C0zRcpBxTðligÞ(cid:1) Bank(PDB)identifier,1D7J)(29),andFKBP-FK506(PDBidentifier,1FKF)
3 (cid:3) qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi(cid:4) (30). FKBP contains 107 residues, BUQ consists of six heavy atoms, and
¼ k 1k ln C +ð2pÞ3=2 ÆDx2æ ÆDx2æ ÆDx2æ ; FK506 consists of 57 heavy atoms. For the FKBP-BUQ simulations, the
2 B B 0 1 c 2 c 3 c proteinatomsweresimulatedwiththePARAM22forcefield(31),andthe
c
ligandparameterswerederivedfromtheanalogousfunctionalgroupsalready
(24)
availableinthePARAM22andPARAM27parameterfiles.FortheFKBP-
FK506 simulations, the CHARMm MSI force field was used and the
where zRBTðligÞ is the same integral as in Eq. 11 with
cpx parameterswereobtaineddirectlyfromtheligand-proteindatabase(LPDB)of
w ¼0,andtheeigenvaluesÆDx2æ ,ÆDx2æ ,andÆDx2æ are
min 1 c 2 c 3 c Rocheetal.(32).TheFKBP-BUQtestcasewasselectedmainlybecauseit
calculated by diagonalizing the covariance matrices of wasthesubjectofapreviousworkwheretheMM/PB-SAmethodwasused
BiophysicalJournal90(3)864–877
AbsoluteProtein-LigandBindingAffinity 869
(17), and therefore, various energy terms could be compared. The FKBP- computationaleffort,lackofprotein-ligandinteractions,andultimateflatness
FK506systemwaschosenbecauseFK506issignificantlylargerthanBUQ ofthePMFcurveatthisrange,welimitedourproductiontimeto500psper
andisanappreciablymorefavorablebindertoFKBP(DDGexp;(cid:2)8kcal/mol). window. WHAM analysis on the resulting three-dimensional probability
bind
Alsotheforce-fieldparameterswerereadilyavailablefromtheLPDB(32). distributionswasperformedwithanin-houseprogramusingahistogrambin
Allsingle-pointsimulationswererunforatotalof2nswherethefirst0.5ns sizeof0.5A˚ toobtaintheunweightedcompositeprobabilitydistribution,
of trajectorywasconsideredequilibrationand the last 1.5 ns was usedfor rðrÞ,asrequiredinEq.8.Aone-dimensionalprobabilitydistributionalong
production.Snapshotsweresavedevery0.1psforatotalof15,000structures thepullingcoordinate,rðlÞ,wasgeneratedbysummingrðrÞovertheplane
persimulation.Thehighfrequencyofsnapshotstoragewasselectedtoen- perpendicular to the pulling direction. Finally, the one-dimensional work
sureconvergenceofthequasiharmonicvibrationalentropymethods.Energy function,WðlÞwasobtainedfromrðlÞinanalogytoEq.8.
averagesweremadeonlyover1-psintervalsforatotalof1500snapshots. In the all of the hybrid solvent simulations, before the dynamics run,
Likewise,therigid-bodytranslation/rotationalentropiesoftheligandinthe explicit solvent molecules were added. For the PMF runs, this was per-
complexwereestimatedusingthese1500snapshots.TheSHAKEalgorithm formedaftertheligandhadbeenshiftedtothecenterofitsrestraintpotential.
wasusedforallhydrogens,andthetimestepofthesimulationwassetto2fs. Water molecules in the x-ray structures were preserved. The additional
ANose-Hooverthermostat(33)withathermalinertiaparameterof10was explicitsolventmoleculeswereputintoplaceusingaprescriptionoutlined
usedtokeepthesimulationtemperatureconstantat298K. elsewhere(35).
All implicit solvent simulations were run with the GB electrostatic Intheendpointstudies,thehybridsimulationswererescoredwiththe
solvation potential GBMV2 (34). Nonpolar (np) solvation was modeled implicitsolventmodeltoestimatebindingfreeenergies.Snapshotsat1-ps
as proportional to the solvent-accessible surface area (SASA), Gnp¼ intervalsstartingfrom500ps(1500structures)wereprocessedbydeleting
gA 1b,whereA isthesolvent-accessiblesurfaceareausingawater theexplicitwatermoleculesandevaluating,withinfinitenonbondedcutoffs,
SASA SASA
probe of 1.4 A˚ and the SASA-1 algorithm (34), g is the surface tension themolecularmechanicsenergyplusGBMV2andSASA-1.Forcomparison
parameterandwassetto5.42cal/(molA˚2),andbwassetto0.92kcal/mol. in the FKBP-BUQ endpoint study, molecular surface-based Poisson sol-
Notethatforsimulatingtwomoleculessuchasacomplex,theconstantb vationenergieswerealsoobtainedat0.5-A˚ gridresolutionusingthePBEQ
moduleinCHARMM(40).
termneedstobetreatedcarefully.Whentheproteinandligandarespatially
Normal mode analysis of the two complexes was performed on 16
separatedfromeachotherbutpartofthesamesimulationasinthePMF
complexstructuresextractedatequallyspacedintervalscoveringtheentire
simulations,asecondbtermmustbeaddedtoaccountforthetwodistinct
productionrun.Twoimplicitsolventmodelswereemployed:distance-based
cavitationpoints.
dielectric (r-dielectric) with a prefactor of 4 and GBMV2/SASA-1. Each
Explicit solvent simulations were performed with the hybrid explicit/
implicitsolventmethod(35)usinganexplicitsolventlayerwidthof10A˚.The structurewasminimizedusingtheadopted-basisNewton-Raphsonmethod
for6000stepsoruntilagradienttoleranceof10(cid:2)6wasreached.Secondde-
hybridsolventmodel(35)involvesencapsulatingabiologicalsolutebyalayer
rivativeswerecalculatedusinganalyticalformulasforther-dielectricrunsand
ofwatermolecules.Outsidethislayer,aGB-basedimplicitsolventreaction
finitedifferencewithastepsizeof10(cid:2)6A˚ fortheGBMV2/SASA-1runs.
fieldisusedtomodelthebulkwatercontinuum.Thehybridsolventmethod
has been shown to be significantly less computationally expensive than
conventional Ewald approaches, primarily because less explicit solvent
moleculesarerequired.However,someminordeviationstotheconventional RESULTS
methodsaretobeexpected(35)becauseofsurfaceboundaryartifacts.Details
onthetheoryandimplementationofthehybridmethodcanbefoundelsewhere PMF calculations
(35,36).Thepairwisemultigridmethod(35,37)wasusedtoapproximatethe
The one-dimensional profiles of the PMF results are pre-
long-rangeelectrostaticandGBterms.Theshort-rangenonbondedcutoffwas
setto12A˚.TheLennard-Jonesinteractionsweretruncatedatthisshortrange sentedinFigs.1and2.InFig.1,themostfavorableposition
withnolong-rangecorrection.Thecenter-of-massandrotationalmomentsof of the center of the ligand for the BUQ/GB simulations is
inertiaoftheproteinwereharmonicallyrestrained(35)withforceconstantsof roughly 1.5 A˚ distance from the most likely spot for the
500kcal/(molA˚2)and500kcal/(molrad2),respectively,topreventtheprotein
hybrid calculation. Also the BUQ/GB simulations result in
frommovingoutofitsfixedshellofwatermolecules.
generating a 1 kcal/mol barrier while pulling the ligand out
The umbrella-sampling procedure involved restraining the protein as
aboveandholdingtheCOGoftheligandatvariouspointsalongalinear ofthepocket.BecausethisbarrierisnotpresentintheBUQ/
pathway. In the standard molecular orientation frame of reference of the hybridresults,itislikelyaresultofthephysicsoftheGB-SA
FKBP protein, the pathway unit vector (in A˚), v, for BUQ was visually implicitsolventmodel.InFig.2,thelocationsofthebinding
chosen to be [(cid:2)0.5243, 0.1068, 0.8448] using the visual molecular minima are roughly the same, but the decays of the inter-
dynamics(VMD)moleculargraphicssoftware(38).FortheFKBP-FK506
actions are different. The FK506/hybrid PMF levels out at
complex,theunitvectorwaschosentobethezaxis.Theinitialboundligand
position,r0 ,wassimplytheCOGoftheligandinx-raycrystalstructure.In 10 A˚. The FK506/GB PMF decays comparatively slower
theFKBPlcigalculations,15simulationswereperformedindependentlywith untilitflattensoutat;12.5A˚.Ingeneral,onewouldexpect
the1kcal/molA˚2COGrestraintpotentialsatpointsdefinedbyEq.8,where thataPMFwouldflattenoutwhenthecorrespondingprotein-
l2½0;1;...14(cid:1). In the FK506 simulation windows were in the range, ligand interaction energy decays to zero because of non-
l2½0;1;...18(cid:1).Duetoinsufficientlysampledregions,whichpertainedto
bonded cutoffs. The fluctuations in Figs. 1 and 2 at larger
barriers in the PMF along the pathway in the FK506 cases, certain
values of the pulling coordinate are likely due to lack of
simulationwindowswereaddedintheGBandhybridsimulations.Forthe
FK506/GBsimulations,simulationwindowswereaddedatl2½1:5;2;2:5(cid:1) statistical convergence. Based on comparisons with a PMF
with 4 kcal/mol A˚2 COG restraints. For the FK506/hybrid simulations, derived from half of the production simulation length, the
windowsatl2½2:5;3;3:5(cid:1)wereaddedwith4kcal/molA˚2COGrestraints errors atany given binare between 0.1 and 0.5kcal/mol.
andl2½5;6(cid:1)with2kcal/molA˚2COGrestraints.Amoreelegantsolutionto Table 3 summarizes the free-energy results for all of the
this problem might involve an automated adaptive umbrella sampling
PMFcalculations.TheBUQresultsforbothsolventmodels
method(39).Also,fortheFK506/hybridsimulations,weadded8A˚ more
waterlayer(foratotalof18A˚)totheligandintherange,l2½10;11...;18(cid:1), are very close to the experimentally determined values.
toavoidsmallbutcumulativesurfaceartifacts(35).Becauseoftheextra The FK506 results are systematically too large by roughly
BiophysicalJournal90(3)864–877
870 LeeandOlson
was chosen led to mild steric clashes of the ligand with the
protein active site. It would certainly be desirable to have
complete automation in a PMF simulation. One could
envision using an adaptive umbrella-sampling procedure
(39) whereby the histograms of the original equally spaced
windowsimulationsareautomaticallyanalyzedtodetermine
wheretherearedeficits insampling.Then,new simulations
wouldbeissuedtofillinthegaps(11).Nevertheless,asfor
computationaltimings,theGB-SAsimulationsforbothsys-
temsrequired8CPU-dayspersimulation window, whereas
the hybrid simulations used ;20 CPU-days per window.
These timings were obtained from simulations run on a
cluster of 32 AMD Athlon MP 22001 processors (AMD,
Sunnyvale,CA).
FIGURE 1 One-dimensionalpotentialofmeanforceforpullingBUQout Endpoint calculations
of the FKBP binding pocket using GB-SA and hybrid explicit/implicit
solvent models. See Theory and Methods sections for simulation details. IncomparisontothePMF calculations,energycomponents
Solid line indicates a hybrid solvent and dashed line signifies GB-SA areextracted fromendpointsimulationsandcombinedwith
implicitsolvent.Arbitraryverticalshiftsofeachcurveweremanuallyad- entropyestimatestoformapproximatebindingfreeenergies.
justedsothattheperceivedasymptoticvaluesofthePMFarelocatednear0.
The enthalpy of complex formation for the two systems,
which for simplicity includes the solvation free energy, is
1–3kcal/mol.TheDGRcpBxTðligÞtermsareverysimilarinboth broken down into components (41) in Tables 4 and 5. As
solvent models. The statistical error associated with the expected, the two-state energy differences have a much
accumulation of errors in the PMF is estimated to be ;0.5 higher degree ofuncertainty than theone-state models. The
kcal/mol, based on the above analysis. Also, in general, mean6SEsintheone-statemodelarequitesmallataround
errors will be better or worse depending on the height of a few tenths of a kcal/mol, and thus, the one-state model is
the barriers in the PMF, which is, in turn, a function of the more advantageous even though it neglects the enthalpy of
pullingdirection.Inaddition,weshouldnotethatitisnearly structuralreorganization (14,16).
impossible to determine what errors, if any, are caused by AnotherobservationisthattheCoulombandGBtermsare
notsampling the lowest free-energy states of the force field numerically large, though nearly compensatory. Their sum,
ateach window. oftendenotedastheelectrostaticsofbinding,playsanimpor-
AspointedoutintheMethodssection,thesimpleoverlay tant role in the free-energy difference (16). The nonpolar
of equally spaced windows was not sufficient for broad surface area terms are uniform among different protocols.
sampling of the FK506 landscape for either solvent model. The vdW terms are quite large, especially in the FK506
Theprimaryreasonforthiswasthatthepullingdirectionthat system. It would appear that these terms must be predom-
inantly compensated by solute entropy in the final free-
energyestimation(9,41).Thetotalenthalpyestimatesforthe
differentmethodsarequitevaried.Somevariationsbetween
one-andtwo-statemodelswereexpected,giventhatthepro-
teinandligandareallowedtorelaxinthetwo-statemodel.
It seems quite unusual that the change in enthalpies
become more favorable for most cases when relaxation is
included.Tounderstandthisphenomenon,we extendedthe
GB-SA simulation times of the FKBP-BUQ system to 7 ns
in both the complexed and unbound states. After 7 ns,
DUgas1DGsolv rose to (cid:2)6.3 kcal/mol, which is higher than
the one-state model ((cid:2)14.3 kcal/mol). The enthalpy of the
complex varied little from the 2-ns value. However, the
protein energy declined by ;14 kcal/mol. This additional
calculationagreeswiththoseofotherstudies(16)inshowing
that 2-ns simulations are too short to obtain converged
DUgas1DGsolv values for two-state models.
FIGURE 2 One-dimensionalpotentialofmeanforceforpullingFK506
Because we used a GB-SA solvent model rather than
outoftheFKBPbindingpocketusingGB-SAandhybridsolventmodels.
SeeFig.1forlegendanddetails. the conventional PB in our rescoring of explicit solvent
BiophysicalJournal90(3)864–877
AbsoluteProtein-LigandBindingAffinity 871
TABLE3 SummaryofresultsforPMFcalculations(kcal/mol)ofthetwoprotein-ligandcomplexesstudiedinthiswork
Ligand Solventmodel w * DGRBTðligÞ DG (correction)y DG (calculation) DG (experimental)
min cpx np bind bind
BUQ GB-SA (cid:2)6.4(0.5) 3.5(0.3) (cid:2)0.9 (cid:2)3.8(0.6) (cid:2)4.5z
BUQ Hybrid (cid:2)7.2(0.5) 2.9(0.3) – (cid:2)4.3(0.6) (cid:2)4.5z
FK506 GB-SA (cid:2)14.6(0.5) 4.2(0.3) (cid:2)0.9 (cid:2)11.3(0.6) (cid:2)12.3§,(cid:2)12.8{
FK506 Hybrid (cid:2)13.6(0.5) 3.8(0.3) – (cid:2)9.8(0.6) (cid:2)12.3§,(cid:2)12.8{
*Errorestimatesinparentheses.
yCorrectionterminGB-SAsimulationsaccountsforsecondcavitationpointinunboundstate.
zDerivedfromexperimentallymeasuredK (29).
i
§Biereretal.(55).
{Biereretal.(56).
simulations, we employed the CHARMM PBEQ module orientational restraints to the system (35), which were in-
(see, e.g., Feig et al. (40) for calculation details) to rescore tended to contaminate the low-frequency rotational modes
the1500structuresineachoftheone-statehybridsimulation withhigh-frequencyrestraintmodes.Ther-dielectricmodes,
sets. In the BUQ set, the PB model obtained an average in contrast, were summed into the entropy by removing the
bindingDGPB¼14.1kcal/mol,whichis1.0kcal/molhigher lowest six eigenvalues in each entropy calculation, which
thantheGB-SAmodel’sresult.IntheFK506set,thebinding roughly correspond totherotationalandtranslationdegrees
DGPB ¼ 42.4 kcal/mol, which is 4.3 kcal/mol higher than of freedom. Another related problem is that strong mini-
DGGB. Note that at the 0.5-A˚ grid resolution used in the mization with the GBMV2 potential can lead to artificially
calculation, the PBEQ Poisson model is only slightly more low minima that are rotationally dependent. Moreover, the
accurate than GBMV2 compared to a benchmark 0.25-A˚ structures may distort slightly due to grid artifacts. These
resolutionPBEQPoissonmodel(40).Also,theGB-SAPMF last two problems are probably somewhat muted by error
results were in reasonable agreement with the explicit cancellation in taking the difference between two entropy
solvent results. For this reason, we attributed discrepancies terms.
withPoissonresultsassimplycontributionstooverallerror. Large differences in absolute entropy estimates exist
By way of the endpoint approach, there are two primary between NM and QH1 (Table 7) for the individual species.
methods for estimating entropy values: normal mode and Fig. 3 shows that NM has more fine structure compared to
quasiharmonic(QH).Normalmodeanalysiswasperformed QH1 (42). In Fig. 4, one can see that over a wide range of
for two potential energy functions as shown in Table 6. As frequencies, the QH1 modes provide larger entropy contri-
onemightexpect,therearesignificantdiscrepanciesbetween butionsthanNM.ItmighthavebeenexpectedthattheQH1
theentropyvaluesusingtheconventionalr-dielectric(RD4) methodwouldhavemorelow-frequencymodesthanNMas
versusGB-SA.Althoughonecancriticizethevalidityofthe seeninFig.3,asmultiplebasinsalongtheenergylandscape
entropy values based on r-dielectric, because it is a crude are explored over the MD simulation run. However, higher
treatment of solvation, the GB-SA entropy terms are not frequency modes (;100–1000 cm(cid:2)1) also seem to contrib-
without reproach. For example, the GBMV2 potential that ute to higher entropy values as seen in Fig. 4. One should
we employed has rotational variance because of the fixed- keep in mind that the QH1 method is void of the high-
orientation angular integration grid used in calculating the frequency covalent hydrogen bond stretching and bending
Bornradii(34).Thismeansthattherearebroadenergywells modes because of the SHAKE constraints imposed on the
corresponding to rigid-body rotations of each chemical system. These high-frequency modes, nonetheless, are not
species.Wetriedtolimittheeffectofthisissue,byapplying expectedtomakeanynoticeablecontributiontotheentropy
estimate.
TABLE4 Averageenergies(kcal/mol)forthedifferent
componentsoftheone-andtwo-stateMM/GB-SAsimulations
oftheFKBP-BUQcomplex TABLE5 Averageenergies(kcal/mol)forthedifferent
componentsoftheone-andtwo-stateMM/GB-SA
GB-SA GB-SA Hybrid Hybrid
simulationsoftheFKBP-FK506complex
Component* one-state two-state one-state two-state
GB-SA GB-SA Hybrid Hybrid
DUinternal 0.0 (cid:2)5.8 0.0 (cid:2)11.3
Component* one-state two-state one-state two-state
DUvdW (cid:2)9.0 (cid:2)13.3 (cid:2)11.0 (cid:2)3.2
DUCoulomb (cid:2)12.7 (cid:2)63.4 (cid:2)9.5 (cid:2)26.7 DUinternal 0.0 2.7 0.0 11.2
DGGB 10.0 67.6 13.1 19.3 DUvdW (cid:2)48.9 (cid:2)50.1 (cid:2)50.2 (cid:2)47.8
DGnp (cid:2)2.6 (cid:2)3.2 (cid:2)2.8 (cid:2)2.2 DUCoulomb (cid:2)15.3 (cid:2)46.3 (cid:2)17.3 5.9
DUgas1DGsolv (cid:2)14.3 (cid:2)18.1 (cid:2)10.1 (cid:2)24.1 DG 28.4 52.3 38.1 7.8
GB
DG (cid:2)6.0 (cid:2)6.9 (cid:2)6.3 (cid:2)6.2
*Mean6SEsfortheone-statemodelsareontheorderofafewtenthsof np
DUgas1DGsolv (cid:2)41.8 (cid:2)48.4 (cid:2)35.7 (cid:2)29.1
a kcal/mol, whereas errors for the two-state models are ;10 kcal/mol
exceptforthenonpolarterm,whichisinerrorby;1kcal/mol. *ErrorestimatesdiscussedinTable4.
BiophysicalJournal90(3)864–877
872 LeeandOlson
TABLE6 Normalmode-derivedvibrationalentropy
components(kcal/mol)fordifferentsystemsaveragedover
16equallyspacedsnapshotsofthe1.5-nsproductionrun
BUQ FK506
Component* GB-SAy RD4z GB-SAy RD4z
TS (cpx) 884.3(1.3) 1218.1(0.9) 973.0(1.4) 1260.8(1.5)
cpx
TS (prot) 876.8(1.2) 1205.9(1.0) 907.6(1.5) 1189.7(0.9)
cpx
TS (lig) 5.6(0.1) 5.4(0.0) 78.0(1.0) 73.4(0.0)
cpx
(cid:2)TDStrans(cid:2)TDSrot 17.1 17.1 21.2 21.2
free free
(cid:2)TDS 15.2(1.8) 10.3(1.4) 33.8(2.3) 23.5(1.7)
bind
*Mean6SEsareshowninparentheses.
yGBcalculationswererunwithtranslationalandrotationalrestraints.Zero
andnegativeeigenvalueswereomitted(seeResultssectionfordetails).
zIntheRD4calculations,thelowestsixvibrationalmodeswereomitted.
The results of the alternative quasiharmonic analysis FIGURE 3 Distribution of modes for normal mode and quasiharmonic
technique,QH2,arealsopresentedinTable7.Theresultsfor vibrationcalculationsfortheFKBP-BUQcomplex.Harmonicmodes(thin
line) are shown for one complex structure from the RD4 calculation.
the two variants of QH are in good agreement with each
Quasiharmonicmodes(thickline)wereobtainedfromthe1.5-nsproduction
other in the case of the BUQ complex given a particular
runoftheGB-SAcomplexsimulation.Frequencybinsstartat1cm(cid:2)1and
choice of solvent model. In contrast, the QH variants differ exponentiallyincreasewithamultiplierof1.1.
quitesubstantiallyintheFK506complex.Onepossibilityis
that FK506 is coupled quite closely with the protein and
therefore the separation of degrees of freedom between the diately obvious how one would choose a priori a particular
proteinandligandinthecomplex(17)maybetoostrongof protocol. Two protocols deserve some comment. First, the
an approximation. The QH results for BUQ with hybrid GB-SAone-stateensemble/GB-SAnormalmodeisarguably
solvent are inconsistent with the results obtained by other self-consistent and reasonably statistically reliable. On the
entropy approaches. Once again, separation of degrees of downside, this protocol lacks an estimate of the configura-
freedom of the solute and water molecules may be an im- tion entropy associated with the multiplicity of basins, and
portant factor, although, this factor alone may not be suf- the relaxation enthalpy and entropy of the free protein and
ficient toexplain suchdisagreements. ligand.Althoughtherelativefree-energyestimateisthebest
When the enthalpic and entropic components are com- ofalloftheprotocols,theabsolutefreeenergiesaretoohigh
bined,awiderangeoffree-energyestimatesemergesasseen by;5kcal/mol.Theothernotableprotocolistheoneclosest
in Table 8. Some of the possible combinations of one- and to that recommended in the article by Kollman and co-
two-state enthalpies and NM, QH1, and QH2 provide good workers (14): the hybrid one-state ensemble using RD4
absolute and relative values. Unfortunately, it is not imme- normalmodeentropy.Thisproceduresucceedsinpredicting
thebindingfreeenergyofFK506,butoverestimatesthefree
TABLE7 Quasiharmonicentropycomponents(kcal/mol)for
thedifferentsimulations
BUQ FK506
Component* GB-SA Hybrid GB-SA Hybrid
TS (cpx) 1470.4 1487.3 1579.3 1567.2
cpx
TS (prot) 1445.7 1460.6 1495.8 1472.2
cpx
TS (lig) 15.3 14.5 82.6 102.6
cpx
TS (prot) 1489.0 1424.3 1526.1 1484.0
prot
TS (lig) 11.7 17.0 93.4 114.7
lig
(cid:2)TDSðprotÞ 43.3 (cid:2)36.3 30.3 11.8
(cid:2)TDSðligÞ (cid:2)3.6 2.5 10.8 12.1
(cid:2)TDSRBTðligÞ 2.6(0.3) 2.0(0.3) 3.3(0.3) 2.9(0.3)
cpx
(cid:2)TDSRBRðligÞ 1.6(0.2) 1.2(0.3) 4.0(0.4) 3.8(0.4)
cpx
(cid:2)TDSRBTðligÞ(cid:2)TDSRBRðligÞ 4.2(0.4) 3.2(0.4) 7.3(0.5) 6.7(0.5)
cpx cpx
(cid:2)TDStrans(cid:2)TDSrot 17.1 17.1 21.2 21.2
free free
(cid:2)TDSQH1 47.4 (cid:2)28.9 61.4 52.7
bind
(cid:2)TDSQH2 43.9 (cid:2)30.6 48.4 30.6
bind
*Trajectorysnapshotsweretakenevery0.1ps.Mean6SEsareshownin FIGURE 4 Contributiontofreeentropyvaluesfromharmonic/quasihar-
parentheses. monicmodes.SeeFig.3forlegendanddetails.
BiophysicalJournal90(3)864–877