Assessment of a Single Decoupling Alchemical Approach for the Calculation of the Absolute

16 Pages·2017·2.61 MB·English

by  Emilio Gallicchio

ORIGINALRESEARCH
published:08March2018
doi:10.3389/fmolb.2018.00022
Assessment of a Single Decoupling
Alchemical Approach for the
Calculation of the Absolute Binding
Free Energies of Protein-Peptide
Complexes
DeniseKilburg1,2andEmilioGallicchio1,2,3*
1DepartmentofChemistry,BrooklynCollege,Brooklyn,NY,UnitedStates,2Ph.D.PrograminChemistry,TheGraduate
Center,CityUniversityofNewYork,NewYork,NY,UnitedStates,3Ph.D.PrograminBiochemistry,TheGraduateCenter,
CityUniversityofNewYork,NewYork,NY,UnitedStates
The computational modeling of peptide inhibitors to target protein-protein binding
interfaces is growing in interest as these are often too large, too shallow, and too
feature-less for conventional small molecule compounds. Here, we present a rare
successful application of an alchemical binding free energy method for the calculation
of converged absolute binding free energies of a series of protein-peptide complexes.
Specifically, we report the binding free energies of a series of cyclic peptides derived
fromtheLEDGF/p75proteintotheintegrasereceptoroftheHIV1virus.Thesimulations
Editedby: recapitulatetheeffectofmutationsrelativetothewild-typebindingmotifofLEDGF/p75,
EmilAlexov, providingstructural,energeticanddynamicalinterpretationsoftheobservedtrends.The
ClemsonUniversity,UnitedStates
equilibrationandconvergenceofthecalculationsarecarefullyanalyzed.Convergenceis
Reviewedby:
aidedbytheadoptionofasingle-decouplingalchemicalapproachwithimplicitsolvation,
NobuyukiMatubayasi,
OsakaUniversity,Japan which circumvents the convergence difficulties of conventional double-decoupling
JanaShen,
protocols.Weherebypresentthesingle-decouplingmethodologyandcriticallyevaluate
UniversityofMaryland,Baltimore,
UnitedStates its advantages and limitations. We also discuss some of the challenges and potential
*Correspondence: pitfalls of binding free energy calculations for complex molecular systems which have
EmilioGallicchio generally limited their applicability to the quantitative study of protein-peptide binding
[email protected]
equilibria.
Specialtysection:
Keywords:protein-peptidecomplexes,bindingfreeenergy,alchemicalcalculations,single-decoupling,implicit
Thisarticlewassubmittedto
solvent,LEDGF/p75,HIVintegrase
MolecularRecognition,
asectionofthejournal
FrontiersinMolecularBiosciences 1. INTRODUCTION
Received:05January2018
Accepted:21February2018 Protein-protein interactions are pervasive in biological systems as they drive and regulate
Published:08March2018 critical functions within the cell (Kastritis and Bonvin, 2013). Concomitantly, there is a strong
Citation: interest in developing therapeutic drugs targeting protein-protein interactions (Higueruelo
KilburgDandGallicchioE(2018) et al., 2009; Basse et al., 2013; Labb et al., 2013; Arkin et al., 2014). The rationale
AssessmentofaSingleDecoupling of using small-molecules or peptides to influence protein-protein binding rests on the
AlchemicalApproachforthe
observation that, even though they may involve large and complex molecular assemblies
CalculationoftheAbsoluteBinding
consisting of thousands of atoms, often protein-protein interactions are commanded by a
FreeEnergiesofProtein-Peptide
Complexes.Front.Mol.Biosci.5:22. relatively small number of interfacial amino acid residues that are part of short linear
doi:10.3389/fmolb.2018.00022 recognition motifs (London et al., 2010). Standard drug design approaches employed
FrontiersinMolecularBiosciences|www.frontiersin.org 1 March2018|Volume5|Article22
KilburgandGallicchio SingleDecouplingforPeptideBinding
for the development of small molecule inhibitors have been obtainedasthedifferenceoftwomuchlargervalueswhichgrow
found generally insufficient compared to peptide constructs in with ligand size and ligand charge and whose statistical and
disrupting protein-protein interactions. Peptides, being larger systematicerrorscombineadditively(DengandRoux,2009).For
thanmostsmallmoleculeinhibitors,cantargetlargeandshallow example, for a protein-peptide complex, the double-decoupling
binding sites with a small ratio of binding energy per atom approachwouldrequirethechallengingcalculationofthelarge
(Kastritis and Bonvin, 2013), which are often categorized as solvationfreeenergyofthepeptidewithinanuncertaintysmall
“undruggable”bystandarddefinition. in comparison to the binding free energy (ideally a fraction
A large variety of approaches are available to model small of a kilocalorie per mole). Furthermore, the free energy of
molecule-protein interactions. These range from ligand-based decoupling from the receptor environment, which arguably
cheminformatic and pharmacophore models (Lavecchia and involvesmuchslowerandcomplexreorganizationprocessesthan
DiGiovanni,2013;Yanetal.,2016),tostructure-baseddocking thedecouplingfromthesolution,wouldhavetobeconvergedto
andscoring(Ewingetal.,2001;Grayetal.,2003;Verdonketal., thesamedegree(Gumbartetal.,2012).
2003; Friesner et al., 2004; Kozakov et al., 2006; Zhou et al., To begin to address these challenges, in this work we
2007;Perrymanetal.,2014;Pierceetal.,2014)andphysics-based assesstheapplicabilityofasingle-decouplingalchemicalmethod
atomisticmodels(Gilsonetal.,1997;Wangetal.,2015;Belletal., (Gallicchioetal.,2010)(referredtohereafteras“SDM”forshort)
2016; Ellis et al., 2016; Zuckerman and Chong, 2017). Which to the calculation of protein-peptide binding free energies. As
are the focus of this work. The ultimate goal of physics-based outlined below, the method employs a λ-dependent potential
computationalmodelsofbindingisthequantitativepredictionof energy function which interpolates between the dissociated
equilibriumconstantsofbinding(or,equivalently,standardfree and associated states of the complex in solution, avoiding the
energies of binding) from statistical mechanics principles and intermediategas-phasestaterequiredbythedouble-decoupling
modelsofinteratomicinteractions(Changetal.,2007;Gallicchio method.Becauseitcomputesthebindingfreeenergydirectly,the
andLevy,2011).Pathway-basedmodelsofbindingmeasurethe method is expected to exhibit favorable convergence properties
free energy changes along a thermodynamic path linking the sufficienttohandlechallengingprotein-peptidesystems.Indeed
unboundandboundstatesofthemolecularcomplex.Thesetake themethodhasbeenalreadytestedwithrelativelylargeligands
theformofphysicalpathwaysinwhichligandandreceptorare includingsometothesamereceptorsite(discussedbelow)which
broughttogetheralongaspatialcoordinate(Héninetal.,2005; istheobjectofthiswork(Lapelosaetal.,2012;Gallicchioetal.,
Gumbart et al., 2013; Jo et al., 2015; Lapelosa, 2017), as well as 2014).
pathwaysinso-calledalchemicalspaceinwhichligand-receptor HIV-1 is a RNA retrovirus that must integrate a reverse
interactions are progressively dialed-in (Chodera et al., 2011; transcribed copy of its RNA genome into host DNA to
Dengetal.,2017). permanently infect a cell (Smith and Daniel, 2006). Essential
Models targeting protein-peptide binding are generally not viral catalytic enzymes include reverse transcriptase (RT) and
asestablished(KilburgandGallicchio,2016).Bioinformaticand integrase(IN).TheformerconvertstheviralRNAgenomeinto
knowledge-based approaches based, for example, on sequence doublestrandedcDNA.Onceinsidethenucleus,INcatalyzesthe
and structure-based homology and phylogenetic profiles, are insertionofviralcDNAintothehostchromosome(Cherepanov
commonlyemployedtosearchforlikelyinteractingproteinpairs et al., 2005a). Integrase catalyzes this reaction in two ways:
(ShoemakerandPanchenko,2007;Zhangetal.,2012).Theability (1) it removes the 3 terminal GT nucleotides from both ends
′
of predicting binding affinities of protein-peptide interactions of the viral cDNA and (2) inserts the newly processed 3
′
fromstructuralmodelsisacrucialchallengethatwouldenhance terminiintohostchromosomalDNA(Cherepanovetal.,2005b).
our understanding of biological regulatory systems and be However, HIV1-IN itself cannot associate with chromatin and
highlybeneficialfordrugdesignanddevelopment.Eventhough requires an endogenous transcriptional coactivator, p75, more
estimating binding free energies from physical principles is a commonly known as LEDGF (lens epithelium-derived growth
dauntingcomputationalproblem,worthwhileattemptsshouldbe factor)(Cherepanovetal.,2005b).Inrecentyears,severalsmall
made to solve it. A main obstacle is that the size and flexibility molecule inhibitors of the HIV1-IN interaction with LEDGF
of peptides make them difficult systems to model (You et al., have been designed by focusing on interactions with selected
2016). Following the work of characterizing physical pathways residues of LEDGF, most notably Ile365, Asp366, and Leu368
approachesforthisproblem(Héninetal.,2005;Gumbartetal., (Cherepanovetal.,2005b;Rhodesetal.,2011).However,known
2013; Jo et al., 2015; Kilburg and Gallicchio, 2016; Lapelosa, mutations of HIV1-integrase have shown high resistance to
2017), the purpose of this work is to assess the applicability of these inhibitors, prompting further study into other potential
alchemicalmethodsfortheestimationofprotein-peptidebinding interactionsbetweenLEDGFandintegrase(Rhodesetal.,2011).
freeenergies. AnextensivestudybyTsiangetal.(2009)reportedthebinding
The double-decoupling method (Gilson et al., 1997), is the affinitiesofnativeLEDGFtoHIV1-INaswellasmanyLEDGF
leading approach for the alchemical calculation of absolute derivedpeptides,includingbothlinearandcyclizedmotifs.One
binding free energies from first principles. It is based on a of the goals of their study was to find the smallest LEDGF-
thermodynamic cycle involving the free energies of decoupling peptidesequencecapableofinhibitingintegrase.Theyfoundthat
the ligand from the solution with and without the presence of small peptides are indeed capable of displacing LEDGF from
the receptor. One of the difficulties of this approach is that HIV1-IN and concluded that cyclizing the peptides resulted in
the binding free energy, which in reversible binding processes greater potency. Building on this, Rhodes et al. (2011) studied
is relatively small and weakly dependent on ligand size, is the interactions between HIV1-IN and multiple linear and
FrontiersinMolecularBiosciences|www.frontiersin.org 2 March2018|Volume5|Article22
KilburgandGallicchio SingleDecouplingforPeptideBinding
boundligands,andcrystallizationionswereremovedandprotein
sidechain protonation states were assigned assuming pH 7 and
absenceofsignificantpKashifts(GluandAspdeprotonatedand
Lys and Arg protonated). A key residue in the LEDGF binding
site, His171, was protonated at the Nδ position as previously
investigated (Gallicchio et al., 2014). The cyclic SLKIDNLD
peptidewasadaptedfromthecrystalstructureofHIVintegrase
from which it was bound, 3AVB, as formerly reported (Rhodes
et al., 2011). Four cyclic mutants (ALKIDNLD, ALKIDNMD,
SLKINNLD,andSLKADNLD)inadditionwithalinearpeptide
H-SLKIDNLD-OH were constructed by modifying the cyclic
SLKIDNLD peptide using the Maestro program (Schrödinger,
LLC).
2.2. Single-Decoupling Binding Free
Energy Protocol
Absolute binding free energies of the HIV-IN with the cyclic
andlinearH-SLKIDNLD-OHpeptidecomplexeswerecomputed
with implicit solvation using the Single Decoupling (SDM)
(Gallicchio et al., 2010) version of Double Decoupling Method
FIGURE1|Representationofthewild-typecyclicpeptide(SLKIDNLD)bound (DDM) (Gilson et al., 1997). SDM was originally introduced
totheLEDGFbindingdomainofHIV1integrase(PDBid:3AVB).Thecyclic astheBindingEnergyDistributionAnalysisMethod(BEDAM)
peptide(green)isshownintuberepresentation.Theintegrasereceptoris (Gallicchio et al., 2010) when it was based on the Weighted
showninsurfacerepresentationcoloredbyelectrostaticcharge.
Analysis Histogram Method (WHAM) (Kumar et al., 1992).
RecentapplicationsofSDMhavebeenbasedonmulti-statebin-
cyclicpeptidesderivedfromtheSLKIDNLD(residues362-369) lessfreeenergyinferencemethods(ShirtsandChodera,2008;Tan
binding motif of LEDGF (Figure1). This motif is located on a etal.,2012).Similarlytodouble-decoupling,inSDMthestandard
loopregion,betweentheα1andα2helicesofLEDGF,whichis bindingfreeenergy1G◦bbetweenareceptorRandtheligandLis
consideredthemaininteractionsite(Cherepanovetal.,2005b). expressedas:
Rhodes et al. found that, in agreement with previous studies
(Tsiangetal.,2009),thelinearpeptide,H-SLKIDNLD-OH,gave 1G◦b =−kBTlnC◦Vsite+1Gb (1)
no activity up to 1 mM, while a cyclized version of the same
where C is the standard concentration of ligand molecules
◦
sequence produced an IC50 value of 70.0 µM using a strand (1M or 1,668 Å 3), V is the volume of the binding site (see
− site
transferinhibitionassay.
below), and 1G is the excess binding free energy, defined as
b
In this work we conduct alchemical single-decoupling
the free energy difference between the coupled state, in which
calculations of the binding free energies of a set of LEDGF-
thereceptorandligandarefullyinteracting,andtheuncoupled
derived peptides to the HIV1-IN receptor. The results
state,inwhichthereceptorandligandareonlyinteractingwith
recapitulate the effect of mutations relative to the wild-type
the solvent and not with each other. In both the coupled and
bindingmotifofLEDGF,particularlythoseatpositionsthathave
uncoupled states the ligand is sequestered in the binding site
beenfoundtobecriticalforbinding,andtheyprovidestructural,
regionasdefinedbelow.
energetic,anddynamicalinterpretationsoftheobservedbinding
Unlike the Double-Decoupling Method (DDM), which
affinities.Theequilibrationofthesystemsandtheconvergence
requires two alchemical calculations going through an
oftheresultsarecarefullyanalyzed.Themethodisfoundtoyield
intermediate “vacuum” state (Gilson et al., 1997; Boresch
reliable binding free energies of cyclic peptides while it fails to
etal.,2003),SDMemploysaλ/temperature-dependentreduced
produceaconvergedestimateforalinearpeptide.Wediscussthe
effective potential energy function which allows for a direct
SDM methodology and the procedure of choosing appropriate
alchemical thermodynamic path between the uncoupled and
simulation parameters. We also discuss some of the challenges
coupledstatesofthecomplex,asfollows:
and potential pitfalls of binding free energy calculations of
complex molecular systems which have generally limited their Uλ(r) β[U0(r) λu(r)], (2)
= +
applicabilitytothequantitativestudyofprotein-peptidebinding
equilibria. whereβ 1/kBT,
=
U (r) U (r ) U (r ) (3)
0 R R L L
2. METHODS = +
is the effective potential energy function corresponding to the
2.1. System Preparation
uncoupledstate,and
TheLEDGFbindingsiteconstructofHIVintegrasewasprepared
aspreviouslyreported(Gallicchioetal.,2014).Watermolecules, u(r ,r ) U(r ,r ) U(r ) U(r ) (4)
R L R L R L
= − −
FrontiersinMolecularBiosciences|www.frontiersin.org 3 March2018|Volume5|Article22
KilburgandGallicchio SingleDecouplingforPeptideBinding
isthebindingenergyfunctiondefinedastheeffectivepotential and
energy difference between the coupled and uncoupled states
∂f (b mx y
i i
of the complex in conformation r = (rr,rL), obtained by ∂b = + σ2 − =0 (10)
translating the ligand from the solvent and into the receptor i (cid:18) i (cid:19)
X
without changing the internal coordinates of either molecule.
Solvingthesetwosimultaneousequationsyields
It is straightforward to show that the λ-dependent potential
energyfunctioninEquation(2)linearlyinterpolatesbetweenthe
(1/σ2) (xy/σ2) (x/σ2) (y/σ2)
u0ntoco1u.pledandcoupledstatesofthecomplexasλisvariedfrom 1S◦b = P i(1/Pσi2)i i(xii2/σ−i2)P− i xii/Pσi2 2i i (11)
In lieu of simulating each λ and/or temperature state
P P (cid:0)P (cid:1)
independently and to improve conformational sampling Theerrorontheentropywascalculatedusingtheerrorsonthe
efficiency,withSDMweutilizeaHamiltonianreplicaexchange individualpoints.
(HREM) λ/temperature-hopping approach where λ and
temperature values are swapped periodically according to a 1 (xi <x>)2
− (12)
Monte Carlo procedure with acceptable exchanges adhering to σ2 = σ2
(cid:18) i (cid:19)
theMetropolisalgorithm(Paletal.,2016). X
To furtherimprove convergenceof thefree energyatvalues Finally,thebindingreorganizationenergy1E wascomputed
reorg
aroundλ 0,weapplyinthisworka“softcore”bindingenergy astheresidualofthereorganizationfreeenergyofbindingafter
=
function: subtractionoftheentropiccomponent:
u′(r) umaxtanh uum(ra)x ifu(r)>0 (5) 1Ereorg =1Greorg+T1S◦b (13)
=(u(r) (cid:16) (cid:17) ifu(r) 0
≤ The errors reported for: 1E and 1G were calculated
reorg reorg
usingstandarderrorpropagation.For1E weusedthestandard
b
whereu issomelargepositivevaluecomparedtothethermal
max error of the mean. All statistical uncertainties, including error
energy. In this work we have set u 1000 kcal/mol. The
max = barsinplots,arereportedastwicethestandarddeviation(96%
energy function modeled in Equation (5) replaces the binding
confidenceinterval).
energy function in Equation (4), when applicable, by limiting
the maximum value of the binding energy for very small λ 2.3. Computational Details
valueswhileleavingthevaluesofthefavorablebindingenergies In this work we employ an effective potential energy function
unchanged. The resulting data is subsequently processed using based on the AGBNP2 implicit solvent model (Gallicchio and
UWHAManalysis(Tanetal.,2012)tocalculatethebindingfree Levy,2004;Gallicchioetal.,2009),andtheOPLS-AAforcefield
energyofthecomplex. (Jorgensen et al., 1996; Kaminski et al., 2001) for the covalent
Thereorganizationbindingfreeenergyiscalculatedbytaking and non-covalent interactions. Parallel molecular dynamics
the difference between the binding free energy and the average simulationswereperformedwiththeIMPACTprogram(Banks
interactionenergyas: et al., 2005). Replica exchange conformational sampling was
conducted for all combinations of eight temperature spanning
1G◦reorg =1G◦b−1Eb (6) 300 to 379 K, and 26 intermediate λ steps at λ = 0.0, 0.002,
0.0048,0.006,0.008,0.01,0.015,0.02,0.0225,0.025,0.03,0.0325,
Additionally,anadvantageofcalculatingthebindingfreeenergy 0.035,0.04,0.07,0.1,0.25,0.35,0.45,0.55,0.65,0.71,0.78,0.85,
multiple temperatures is the availability of the conformational 0.92,and1,foratotalof208replicas.Thebindingsitevolume
entropyofbinding1S was defined as any conformation in which the peptide center
◦b
of mass is within 6 Å of the center of mass of the Cα atoms of
∂ residues 168-174 and 178 of chain A and residues 95-99, 102,
1S◦b =−∂T1G◦b (7) 125,128,129,and132ofchainB(residueandchaindesignations
according to 3NFB crystal structure) of HIV1-IN. The peptide
As we have only several binding free energies at discrete wassequesteredwithinthebindingsitebymeansofaflat-bottom
temperatures,theabovedifferentialwassolvedbyusingalinear harmonic potential with a force constant of 3.0 kcal/mol Å2
leastsquaresapproximation. applied to atoms with distances greater than 6 Å. The volume
of the binding site is calculated to be 904 Å corresponding to
b mxi yi 2 kBTlnC◦Vsite = 0.69kcal/mol.Asinpreviouswork(Gallicchio
f = + σ − (8) et al., 2014), the Cα atoms of the residues of the CCD domain
i (cid:18) i (cid:19) ofHIV1-INwererestrainedtotheoverallstructurebymeansof
X
sphericalharmonicrestraintswithaforceconstantof1kcal/mol
Å2, excluding the residues closest to the binding site (residues
∂f (b mxi yi)xi 167-178 of chain A and 95-102 and 124-132 of chain B) which
+ − 0 (9)
∂m = σ2 = wereleftunrestrained.
i (cid:18) i (cid:19)
X
FrontiersinMolecularBiosciences|www.frontiersin.org 4 March2018|Volume5|Article22
KilburgandGallicchio SingleDecouplingforPeptideBinding
We simulated the complexes of HIV1-IN with six peptides:
cyclicSLKIDNLD,linearH-SLKIDNKL-OH,cyclicALKIDNLD,
cyclicALKIDNMD,cyclicSLKINNLD,andcyclicSLKADNLD.
Multi-dimensionalreplicaexchangecalculationswereperformed
forabout3.3nsofmoleculardynamicsperreplica,or867nsof
averagesimulationtimeforeachcomplexandapproximately5.2
µs of simulation time in total. Binding energies were sampled
withthefrequencyof25psforatotalof35,000bindingenergy
samples per complex. All binding free energy calculations were
conducted on Louisiana State University’s XSEDE SuperMic
supercomputerusingtheASyncREjobdistributionmiddleware
(Gallicchioetal.,2015).
2.4. Error Analysis and Determination of
Equilibration and Convergence
FIGURE2|Representativereversecumulativeconvergenceplotforthe
Statistical uncertainties of binding free energy estimates were
bindingfreeenergyforthecomplexwithpeptideMutant2(S362A/L368M)
calculated using UWHAM’s built-in Fisher’s analysis. The highlightingthedifferentmethodsforchoosingequilibrationtime.
UWHAM formalism, based on mathematical distribution (1)representstheequilibrationtime,teq,returnedbythequantitativeprotocol
theory, uses the curvature of the likelihood function to weadoptedinthisworkasdescribedinthetext,(2)representstheteqchosen
qualitativelybasedonthefirstobservedinflectionpointoftheplotand(3)
calculate the variance of free energy estimates. A detailed
representsthesmallestteqcorrespondingtoafreeenergyestimate
derivation is provided by Tan et al. (2012) Fisher’s analysis statisticallyindistinguishablefromthoseatlongerequilibrationtimes.The
assumes uncorrelated data, we therefore assessed the statistical resultant1G0b(teq)are−7.7±0.2,−7.5±0.2,and−7.1±0.4kcal/molwith
inefficiencyofthebindingenergytimeseriesdatabycomputing thethreemethods,respectively.Thethreemethodsyieldstatisticallyequivalent
uncertaintiesbasedonastandard1/√N erroronprogressively resultsinthiscase.
smaller, random sets of the data and then subsequently
plotted those errors against the error provided by UWHAM
that contained the same number of data points. From these A suitable equilibration time can be chosen qualitatively as,
comparisons, as well as the data obtained from running auto- for example, the smallest equilibration time that gives a free
correlation analysis, we have concluded that binding energies energy estimate statistically indistinguishable from those at
were collected with sufficiently small frequency so as to make longer equilibration times (see Figure2) (Yang and Karplus,
statisticalcorrelationsnegligible. 2003).
Most molecular dynamics simulations are initiated with In this work we have explored the quantitative protocol
structures that are atypical of equilibrium conformations. As for the choice of the equilibration time recently proposed for
binding energy calculations are sensitive to small perturbations averagesofcorrelatedtimeseries(Chodera,2016).Whenapplied
inconfiguration,itistypicalpracticetoremoveaninitialportion to a generic time series xi x(ti), this method consists of
=
ofthetrajectoryinwhichthesystemisapproachingequilibration picking the equilibration time so as to maximize the effective
so as not to adulterate the equilibrium result (Klimovich et al., number, n , of statistically independent samples of the time
eff
2015). To determine the amount of initial data to eliminate, in series(AllenandTildesley,1993)whichremainsafterremoving
this work we employ a method similar in the spirit of reverse thefirstneqdatapointscollectedpriortothechosenequilibration
cumulative averaging from Yang and Karplus (2003) and the time (see Figure3 for an example). The number of statistically
autocorrelation analysis discussed recently by Chodera (2016) independentsamplesisgivenby
In this approach, we examine the time series of binding free
n n
energyestimatesasafunctionofincreasingequilibrationtimeteq n − eq , (14)
eff
measuredfromthebeginningofthesimulation.Specifically,we = g(neq)
define1G [t ]asthebindingfreeenergyestimateobtainedby
◦b eq
wherenisthetotalnumberofpointsinthetimeseries,n n is
discarding initial data up to simulation time t . The sequence eq
eq −
thenumberofpointsremaininginthetimeseriesafterremoving
of binding free energy estimates so obtained is referred to as
thefirstn datapoints,and
the reverse cumulative profile (Yang and Karplus, 2003). In eq
these profiles (see Figure2 as an example) binding free energy
g(n ) 1 2τ (15)
estimates to the right, at long equilibration times, are those eq = + eq
leastaffectedbybiasintroducedbyunequilibrateddataatearly
isthestatisticallyinefficiencyofthetimeseries,whereτ isthe
times. Conversely, estimates to the left, at short equilibration eq
correlationlengthgivenby:
times,areexpectedtobethemostbiased.Thedecisionofwhich
equilibrationtimetopickisnotobviousbecausetheleastbiased
estimates,sincetheycorrespondtothesmallestbindingenergy tmax−1 t
τ 1 C (16)
eq t
samples, are also those with the largest statistical uncertainties. = t 1 (cid:18) − tmax(cid:19)
X=
FrontiersinMolecularBiosciences|www.frontiersin.org 5 March2018|Volume5|Article22
KilburgandGallicchio SingleDecouplingforPeptideBinding
FIGURE3|(A)Reversecumulativebindingfreeenergyconvergenceplotfor
FIGURE4|(A)Reversecumulativeconvergenceplotand(A)corresponding
thecomplexwithpeptideMutant4,I365A.Thearrowpointstothecomputed
neffplotforthebindingfreeenergyofthecomplexofHIVintegrasewiththe
equilibrationtime.(B)Plotoftheeffectivenumberofindependentsamples,
SLKIDNLDlinearpeptide.Bothfiguresshowthatneitheraqualitative(A)nora
neff,asafunctionofequilibrationtimeforthesamecomplex.Themaximum quantitativeapproach(B)providesavalidequilibrationtimeforthissystem.
indicatestheestimatedequilibrationtimeof1.9ns.
Theconclusionisthatthebindingfreeenergyforthissystemhasnotreached
equilibration.
wheretheautocorrelationfunction,C isdefinedas:
t
caseofmutant4(I365A),forexample,thisprocedurereturnsan
equilibrationtimeof1.9nsperreplica.
x x x 2
Ct = h nhxn2n+it−i−hxnhi2ni (17) smaInll egqenuielribalr,atwioenotbimserevsefotrhastysttheemosutthliantewdopurlodcbedeucroenrseidtuerrends
to have equilibrated quickly based on qualitative arguments,
It was found that this approach incorporates the trade-off, statisticaluncertainties,andtheshapeofthereversecumulative
discussed above, between minimizing the variance while at the profiles. Conversely, the procedure returned relatively long
same time minimizing the bias in the estimate. The variance equilibration times (see Results) for systems that qualitatively
obviously decreases as the number of samples included in the appear to equilibrate slowly, for example those that display
estimateincreases,thatisasn decreases.Incontrast,Chodera persistent drift of the reverse cumulative profile for short
eq
observed that the correlation length, and, consequently, the equilibrationtimes(Figures3,4).Inthecaseofalinearpeptide,
statisticalinefficiency,growasatypicalsamplesatthebeginning considered from multiple perspective to have failed to achieve
of the simulation are included (Chodera, 2016), presumably equilibration, this procedure correctly indicated the lack of a
becausetheseintroduceanoveralldriftinthetimeseries. reasonableequilibrationtimeasthen curvedidnotdisplaya
eff
We obtained reasonable estimates of equilibration times by clearmaximumasshowninFigure4.
applyingtheapproachabovetothereversecumulativeprofileof Having determined the equilibration time, we then assessed
thebindingfreeenergy1G [t ].Werationalizethisobservation convergence of the free energy calculations by examining the
◦b eq
bynotingthat,whilethereversecumulativeprofileisnotatime conventional forward cumulative profiles (see Figure5 and
series of an observable, n still captures the trade-off between Figure S1 for an example) obtained considering only binding
eff
maximization of the number of samples and minimization of energy samples collected at times following the equilibration
the statistical inefficiency, which is larger at small equilibration time until some maximum time t, which is varied. Systems
timesduetheoveralldriftofthereversecumulativeprofileinthat that are considered converged display a forward cumulative
region (see Figure3). Specifically, we considered the family of profile with monotonically decreasing uncertainties, ideally
reversecumulativefreeenergyprofiles1G [t>t ],thatisthose randomfluctuations,or,moreoften,smalldownwardorupward
◦b eq
derivedbythefullreversecumulativeprofileafterremovalofthe driftscontainedwithincomputeduncertainties.Tobetterassess
valuesfort<t .Foreachprofileweperformedautocorrelation convergence, it has also been helpful to consider the quantity
eq
analysistocomputethecorrespondingcorrelationlengthτ(t ), (δ1G0/δt), which measures the differential effect of adding
eq b
the statistical inefficiency g(t ) 1 2τ(t ) and effective samplescontainedbetweentandt 1t.Asshownforexample
eq eq
= + +
number of statistically independent samples n (t inFigure5,inconvergedcalculationsthisfunctionprogressively
eff max
∝ −
t )/g(t ). We then selected the optimal equilibration time as decays to zero with smaller and smaller fluctuations with
eq eq
that one that maximizes n . As illustrated in Figure3, in the increasingsimulationtime.
eff
FrontiersinMolecularBiosciences|www.frontiersin.org 6 March2018|Volume5|Article22
KilburgandGallicchio SingleDecouplingforPeptideBinding
3. RESULTS Consistent with experimental findings (Tsiang et al., 2009;
Rhodes et al., 2011), mutation of the critical D366 residue
3.1. Thermodynamic Parameters
(M3 in Table1) resulted in a significant reduction of binding
Computed thermodynamic parameters for the binding of five
affinity. This residue is known to form a strong hydrogen
cyclicLEDGF-derivedpeptidesconsideredinthisworkarelisted
bonding interaction motif with the backbone atoms of His171
inTable1.Thecomplexeswithwildtype(SLKIDNLD),mutant1
andGlu170ofHIVintegrase(Figure8),whichisalsofoundin
(ALKIDNLD), and double-mutant 2 (ALKIDNMD) have been
manysmall-moleculeinhibitors(Gallicchioetal.,2014).Theloss
studied by Rhodes et al. (2011) who obtained IC values and
50 ofthisinteractionisobservedtoleadtomajordisruptionofthe
crystalstructures.Mutants3(SLKINNLN)and4(SLKADNLD)
peptide-receptor binding interface and loss of binding affinity.
have been included in this study to probe the thermodynamic
Analogously,mutationofI365toalanine(M4inTable1)causes
andstructuralconsequencesofmutationsinresiduesknownto
disruption of favorable hydrophobic packing peptide-receptor
be critical for successful binding (Tsiang et al., 2009). We have
interactions and loss of binding affinity, although to a smaller
been unsuccessful at computing the binding free energy of the
degreethantheD366Nmutation.
linear wild typepeptide (see below), which showed weak or no
Table1 also reports the decomposition of the computed
activity in experimental tests (Tsiang et al., 2009; Rhodes et al.,
binding free energies into average binding energies and
2011).
reorganization free energies, and the further decomposition of
Thecomputedbindingfreeenergiesforwildtypeandmutant
thelatterintoreorganizationenergyandconformationalentropy
1cyclicpeptidesarequalitativelyconsistentwithµM-rangeIC
50 of binding. Average binding energies (1E , third column in
b
values(Rhodesetal.,2011)measuredforthesepeptides(85.2
± Table1) reflect the strength of direct interatomic interactions
24.5 µM and 39.7 7.1 µM, respectively). Also in agreement
± betweenthepeptideandtheproteinreceptor.Thesearerelatively
withinhibitionexperiments,mutationoftheS366serineresidue
constant around 58 kcal/mol across the peptide set, with the
toalaninecausesminorchangesinbindingaffinity.Thedouble −
exceptionofmutant3,forwhich1E isabout20kcal/molless
b
mutant S366A/L368M (M2 in Table1), while undergoing a
favorable than the others. This behavior confirms that the loss
greaterconformationalchangerelativetowildtypeandmutant
ofbindingaffinityformutant3ismainlyduetothelossofkey
1 (Figures7, 8), is predicted to have similar binding affinity to
hydrogen bonding interactions between D366 and a histidine
the wild type and the single mutant (Figure6B). Experimental
residue in the binding pocket of the receptor (Figure8). The
values for this double-mutant are contradictory. Rhodes et al.
magnitude of the interaction energy loss (20 kcal/mol) is also
(2011) report no detectable inhibition activity for this peptide,
consistent with the loss of two strong hydrogen bonds and the
while they were able to obtain a crystal structure of it bound
correspondingfavorableelectrostaticinteractions.
to HIV integrase (PDB id 3AVJ) (Rhodes et al., 2011). Taken
In contrast, the weakening of the binding of mutant 4 (M4
together, this evidence suggests that, while mutant 2 (M2) is
in Table1), in which a bulky isoleucine residue is replaced
capableofbindingtoHIVintegrase,somehowitcannotinhibit
by alanine, displays a completely different thermodynamic
iteffectively.
signature.Theaveragebindingenergyforthispeptidecomplex
is of similar magnitude than the higher affinity peptides, and
it is actually predicted to be greater than that of mutant
1 even though their predicted binding constants differ by
at least two orders of magnitude. The weaker binding of
mutant 4 is predicted to be mainly due to a large unfavorable
reorganization free energy (fourth column in Table1), which,
in turn, is caused by unusually large intramolecular energy
strain(1E )onlypartiallycompensatedbyasmallerentropic
reorg
penalty (T1S0). As further discussed below, this behavior
b
is consistent with the observed structural response of the
peptide and receptor attempting to compensate for the loss
of favorable packing interactions (Figure8D). In contrast, the
thermodynamic parameters collected in Table1 indicate that
the loss of hydrogen bonding interactions in mutant 3 results,
effectively, in the unbinding of the peptide with little or no
tendency for conformational reorganization aimed at regaining
bindingasshowninFigure8C(seebelowforadiscussionofthe
observedstructuralchanges).
3.2. Structural Considerations
The equilibration times of the complexes are found to be
FIGURE5|(A)ForwardcumulativeconvergenceplotforMutant4afterthe highly correlated with the extent of their conformational
initialbiashasbeenremoved.(B)Plotofthefluctuationsof1G0 thatshowsa
B reorganization to go from the starting structure (which in all
steadydecaytowardzero.
cases was the crystal structure of the wildtype peptide) to
FrontiersinMolecularBiosciences|www.frontiersin.org 7 March2018|Volume5|Article22
KilburgandGallicchio SingleDecouplingforPeptideBinding
the equilibrated conformational ensemble. Thus, as expected, by an hydrogen bond with the hydroxyl group of Thr174 and
the binding free energy of the wildtype peptide converged withasaltbridgewiththeprotonatednitrogenatomofHis171.
relatively quickly (Figure9B). Conversely, the I365A mutation BothSer362andLys364ofthepeptideformsstronginteractions
causedextensiveconformationalrearrangementswhichrequired withthesidechainofGlu170.Ile365forms backbone-backbone
alongerequilibration(Figure6D). interactionswithGln168andThr174.ThesidechainofIle365is
Consistentwiththecrystalstructure(3AVB),residuesSer362, additionallynestledintoanhydrophobicpocketofthereceptor
Lys364, Ile365, Asp366, and Asn367 of the cyclic WT peptide lined by residues Ala128, Trp131, Trp132, and Met178. The
form the primary intermolecular interactions with the receptor solvent-exposed residues of the peptide (Leu363, Leu368, and
(Figure7).ThecarboxylategroupofAsp366ofthepeptideforms Asp369),donotinteractconsistentlywiththereceptoralthough
the strong hydrogen bonds with the NH backbone groups of they make occasional intramolecular interactions with each
Glu170 and His171 of the receptor which are the hallmark of other.
theLEDGF-HIVINcomplex(Cherepanovetal.,2005a)aswell A serine to alanine mutation of residue 362 (M1 mutant)
as those with peptido-mimetic synthetic inhibitors (Peat et al., causes a slight, but statistically significant decrease in the
2014). The carboxylate group of Asp366 is further anchored predictedbindingaffinity(Table1).Auto-correlationanalysisof
TABLE1|ComputedthermodynamicparametersforthebindingbetweenHIV1-integraseandaseriesofcyclicpeptidesderivedfromtheLEDGFprotein.
Peptide 1G◦ba 1Eba 1Greorga 1Ereorga T1S◦ba
SLKIDNLD(WT) 8.6 0.4 59.3 0.5 50.6 0.5 23.6 1.9 27.0 1.8
− ± − ± ± ± ±
ALKIDNLD(M1) 8.1 0.2 56.9 0.4 48.8 0.4 22.1 0.4 26.7 0.1
− ± − ± ± ± ±
ALKIDNMD(M2) 7.6 0.2 59.5 0.5 51.9 0.5 24.6 0.8 27.3 0.6
− ± − ± ± ± ±
SLKINNLD(M3) 2.7 0.2 36.1 0.6 38.3 0.6 12.4 0.6 25.9 0.6
± − ± ± ± ±
SLKADNLD(M4) 3.7 0.2 57.6 0.5 53.9 0.5 33.2 0.6 20.7 0.4
− ± − ± ± ± ±
aInkcal/mol.
FIGURE6|Reversecumulativebindingfreeenergyconvergenceplotsfor(A)Mutant1(S362A),(B)Mutant2(S362A/L368M),(C)Mutant3(D366N),and
(D)Mutant4(I365A).Arrowsmarktheestimatedequilibrationtimesbasedonautocorrelationanalysis:(A)0.20ns,(B)0.25ns,(C)0.20ns,and(D)1.90ns.
FrontiersinMolecularBiosciences|www.frontiersin.org 8 March2018|Volume5|Article22
KilburgandGallicchio SingleDecouplingforPeptideBinding
mutants,themostsignificantstructuralrearrangementoccurred
for the receptor. The void in the hydrophobic binding pocket
left from the isoleucine to alanine mutation slowly collapsed
around the alanine sidechain (Figure8D). Concomitantly, the
peptide moves deeper into the binding cavity causing the
solvent-exposed residues (Leu363, Leu368, and Asp369) to
make increased interactions with the receptor while slightly
decreasing interactions between Asp366 and Glu170/His171.
These rearrangements recover most of the interaction energy
of binding (Table1, third column) at the expense of higher
conformationalstrain(fifthcolumn),resultinginabindingfree
energy5kcal/mollessfavorablethanwildtypeinagreementwith
experimental evidence (Cherepanov et al., 2005b; Tsiang et al.,
2009).
FIGURE7|Representativestructureextractedfromthesimulationofthe
wildtypecyclicpeptide(WTinTable1)boundtotheLEDGF-bindingdomain 3.3. The Potential Pitfalls of
ofHIVintegrase.Thefollowingintermolecularinteractionsarehighlighted: Protein-Peptide Binding
Lys364-Glu170,Asp366-Thr174,His171,andAsn367-Gln95.Alsoshownin
theintramolecularinteractionbetweenAsp168andGln169ofthepeptideand Alchemical binding free energy calculations require careful
thekeybidentatedinteractionbetweenAsp366andthebackbonenitrogen planningandassessmentoftheresults(Klimovichetal.,2015).
atomsofGlu170andHis171. Thisisparticularlysoforprotein-peptidecomplexesduetotheir
largeconformationalvariability.Herewediscussthechallenges
weencounteredinapplyingourcalculationprotocols,tunedfor
smallmoleculebinding,topeptides.
thereversecumulativebindingfreeenergyprofiles(Figure6A) The first challenge we encountered turned out to be related
revealed only a small equilibration bias of about 1 kcal/mol, to our initial choice of the λ schedule, which was based on the
consistent with the relatively minor conformational changes one that had been used with success on large ligands for the
relativetothewildtypepeptide(Figure8A).Themostnoticeable samereceptorsite(Gallicchioetal.,2014).Withthisλschedule
structural difference is the shift of the interaction between we collected a substantial amount of binding free energy data
Glu170andSer362,nowabsent,toAsn367andtheconcomitant that, when displayed using a traditional forward cumulative
weakening of the Lys364–Glu170salt bridge which now occurs plots showed no particular disturbance (Figure9B) beyond a
less often. These subtle changes are the cause of the weakening persistent drift signaling slow convergence. It was only when
of the binding energy and the slight less favorable binding free we analyzed the reverse cumulative plot (Figure9A) that we
energyrelativetothewildtype(Table1). discoveredthattheslowconvergencewasasymptomofaserious
ThedoublemutationS362A:L368M(M2inTable1)resulted numericalartifact.
in no significant change of binding free energy relative to Thereversecumulativeplot(Figure9A)showsanabruptshift
the S362A single-mutant. However, the methionine sidechain in the binding free energy corresponding to an equilibration
at position 368 now affords additional polar interactions with timeofapproximately1ns.Subsequentexaminationoftheerror
Thr125 and it induces the recruitment of Asp369 (previously bounds returned bytheUWHAM(Tanetal.,2012)procedure,
completelysolventexposed)inthesameinteraction(Figure8B). revealed that without the equilibration data bias introduced by
Asaresult,thepeptideshiftsoverallpositiontowardthe124-131 the binding energy data collected at times prior to 1 ns, the
α-helixofthereceptor. uncertaintyinthecomputedbindingfreeenergywassolarge(on
TheD366Nmutationresulted,essentially,inthedissociation theorderof1000kcal/mol)astomaketheresultsmeaningless.
of the complex (Figure8C), consistent with the unfavorable It was eventually established that the failure in obtaining a
computed binding free energy (Table1, M3 mutant) and reasonable binding free energy estimate was due in this case to
previous experiments (Tsiang et al., 2009). The peptide was thechoiceoftheλschedule,whichlackedsufficientintermediate
kept in proximity of the pocket only by the applied tether (see λ values between 0.1 and 0.4 to ensure sufficient ensemble
Methods). The conformations of the receptor and the peptide “overlap” of adjacent alchemical states (Klimovich et al., 2015).
rapidly reorganized to screen the unmet positive charge of Specifically, after some equilibration, for λ < 0.1 all of the
His171whiletheAsn366sidechainshiftedawayfromtheGlu170 binding energy values clustered around the maximum value
backbone. The rapid equilibration of the peptide, freed from u [seeEquation(5)],whichwerenotobservedinanyofthe
max
strong intermolecular interactions, is the cause of the smaller other λ states. As a result, the UWHAM multi-state maximum
equilibration time observed (Figure6C) for the M3 mutant likelihood function (Tan et al., 2012) lacked a well-defined
relativetotheothermutants. maximumcorrespondingtothebindingfreeenergyprofile.This
The I365A mutant (M4 in Table1) displayed the most problemhadnotbeenapparentwhenunequilibrateddataatthe
significant and lengthy equilibration of the binding free energy beginningofthesimulation,whichcontainedasmallamountof
among the peptides studied (Figure6D). Unlike the other binding energy data bridging alchemical states across λ 0.1,
=
FrontiersinMolecularBiosciences|www.frontiersin.org 9 March2018|Volume5|Article22
KilburgandGallicchio SingleDecouplingforPeptideBinding
FIGURE8|Representativestructuresextractedfromthesimulationsof(A)theS362Amutantcyclicpeptide(M1inTable1)and(B)theS362A:L368Mdouble
mutantcyclicpeptide(M2),(C)theD366Ncyclicpeptide(M3),and(D)theI365Amutantcyclicpeptide(M4)withtheHIVintegraseLEDGF-bindingdomain.All
interactionsshownwerepreviouslydescribedinFigure7.Noticein(C)thattheD366Nmutationcausesseveredistortionofthepeptideaswellasthereceptor
relativetotheboundconformationofthecomplex.
wasincludedintheestimate,resultinginaforwardcumulative linear form of the wild-type peptide (H-SKIDNLD-OH). As
convergenceplotnotclearlyoutoftheordinary. showninFigure4,thereversecumulativeplotforthecomplex
As a result of this assessment, we modified the λ scheduled with this peptide does not show a plateau, indicating that the
to insert intermediate λ states between 0.1 and 0.4 more likely system was still in the equilibration phase up until the longest
to sample both large and smaller binding energy values. This simulationtimeconsidered.Lackofconvergenceisconfirmedby
resulted in binding free energy estimates with reasonable error theforwardcumulativeplotincludingallofthebindingenergy
boundswhenanalyzedwithbothforwardandreversecumulative samples(Figure11A).
plots(Figures9C,D). Analysis of molecular dynamics trajectories revealed that
This experience highlights the need for careful examination the lack of convergence was due to the slow equilibration
of equilibration artifacts and error bounds when assessing the of the relative populations of two competing conformational
resultsofalchemicalbindingfreeenergycalculations.Clearly,in states of the complex (Figure11B). The first state corresponds
thiscase,thechoicesofsimulationlengthandλscheduletuned to a compact conformation of the linear peptide tethered by
forsmall-moleculebindingwerenotsuitableforthemodelingof an intramolecular salt bridge similar to the conformation of
the binding of the larger peptides studied here. The optimized the cyclic peptide. The other, more extended, conformation
λ schedule addressed the requirement of a suitable, unbroken, resultedfromthereplacementoftheintramolecularsaltbridge
thermodynamicpathbetweenthecoupledanduncoupledstates with intermolecular salt bridges between the termini of the
of the complex and enabled the assessment of equilibration peptide and residues lying on the rim of the protein receptor
bias in the binding free energy estimates. However it did not pocket (Figure10). This extended bound conformation, whose
guaranteeconvergenceofthebindingfreeenergyestimateswhen formation was found to be reproducible across independent
these are affected by challenges of a different nature, as the runs,isstronglyfavoredbybindingenergy(approximately 120
−
followingexampleillustrates. kcal/mol,comparedwithbindingenergiesofhalfthismagnitude
One of aims in this work had been to study the forcyclicpeptidesandforthecompactconformationofthelinear
thermodynamic factors that favor the binding of the cyclic peptides–seeTable1andFigure11–andstronglydisfavoredby
wild-typepeptideoverthelinearone,H-SKIDNLD-OH(Tsiang reorganizationfreeenergy.
et al., 2009; Rhodes et al., 2011) (Figure10). However, we Becauseofthetrade-offbetweenenergeticandreorganization/
were unable to obtain a converged binding free energy for the entropic opposing thermodynamic driving forces, the compact
FrontiersinMolecularBiosciences|www.frontiersin.org 10 March2018|Volume5|Article22