Bioorganic & Medicinal Chemistry 7 (1999) 851–860
`
`Estimation of the Binding A(cid:129)nities of FKBP12 Inhibitors Using
`a Linear Response Method
`
`Michelle L. Lamb,y Julian Tirado-Rives and William L. Jorgensen*
`
`Department of Chemistry, Yale University, New Haven, CT 06520-8107, USA
`
`Received 15 September 1998; accepted 2 November 1998
`
`Abstract—A series of non-immunosuppressive inhibitors of FK506 binding protein (FKBP12) are investigated using Monte Carlo
`statistical mechanics simulations. These small molecules may serve as sca€olds for chemical inducers of protein dimerization, and
`have recently been found to have FKBP12-dependent neurotrophic activity. A linear response model was developed for estimation
`of absolute binding free energies based on changes in electrostatic and van der Waals energies and solvent-accessible surface areas,
`which are accumulated during simulations of bound and unbound ligands. With average errors of 0.5 kcal/mol, this method pro-
`vides a relatively rapid way to screen the binding of ligands while retaining the structural information content of more rigorous free
`energy calculations. # 1999 Elsevier Science Ltd. All rights reserved.
`
`Introduction
`
`The binding protein of the immunosuppressant natural
`product FK506 (Fig. 1) has been the target of extensive
`investigation by both biochemical and theoretical tech-
`niques during the last 10 years. Discovery of the cis-
`trans peptidyl-prolyl
`isomerase (PPIase or rotamase)
`activity of FKBP12 (MW=12 kDa) led to dissection of
`the rotamase mechanism1–4 and hopes for the rapid
`design of low molecular weight immunosuppressant
`molecules that inhibited this activity. The crystal struc-
`ture of FK506 bound to FKBP12 was crucial to this
`endeavor, as it demonstrated that the peptidomimetic a-
`ketoamide and pipecolyl portions of the ligand were
`buried but
`that much of
`the macrocycle remained
`exposed to solvent.5 However, in later studies rotamase
`inhibition was found to be insu(cid:129)cient for immunosup-
`pression; the FK506–FKBP12 complex associates with
`the surface of calcineurin (CN), a serine/threonine
`phosphatase, and hinders binding of subsequent pro-
`teins in the T-cell signaling pathway.6 The ‘‘e€ector’’
`region of FK506, which contacts CN, is opposite the a-
`ketoamide-pipecolic acid moiety.7,8 Thus, the PPIase
`inhibitors developed through structure-based design
`e€orts (e.g. compounds 1–7 in Table 1) formed a set of
`potential
`sca€olds
`for
`immunosuppressive
`e€ector
`components.9,10
`
`Key words: Monte Carlo; linear response; FKBP12; rotamase inhibi-
`tors.
`*Corresponding author. Fax: +203-432-6299; e-mail: bill@adrik.chem.
`yale.edu
`y Current address: Dept. of Pharmaceutical Chemistry, University of
`California, San Francisco, Box 0446, San Francisco, CA 94143-0446.
`
`these non-immunosuppressive FKBP12
`Many of
`ligands are also able to promote neuronal growth in
`vitro and in vivo through binding to FKBP12.11–13
`Dose-response studies of neurite outgrowth in chick
`dorsal root ganglia resulted in an ED50 of 0.058 nM for
`for example.11 In addition,
`tethered
`compound 10,
`dimers of molecules similar to 4 have been used to
`‘‘chemically-induce’’ dimerization14 of targeted cellular
`proteins that had been adapted to include FKBP12
`domains.15 Modification of targets in appropriate sig-
`naling pathways can result in apoptosis or the induction
`of transcription; the prospects of this technology for
`gene therapy are intriguing.
`
`Consequently, we have used theoretical techniques to
`probe this class of small molecules (Table 1) to gain
`insight
`into the physical basis
`for di€erences
`in
`FKBP12-binding activity and to evaluate new methods
`for the calculation of protein–ligand binding a(cid:129)nities.
`In previous work,16 relative free energies of binding
`((cid:1)(cid:1)Gb) for compounds 2, 4, and 8–10 were calculated
`by the free energy perturbation (FEP) technique using a
`Metropolis Monte Carlo (MC) algorithm for config-
`urational sampling.17 The parameters and geometry of
`one ligand were transformed into those of another with
`these simulations, both in solution and while bound to
`the protein. At each step of the transformation, the
`system was brought to equilibrium and the free energy
`di€erence relative to the previous step was computed.
`The di€erence in binding free energy for the two ligands
`was then found from the di€erence in the total free
`energy change for each (bound and unbound) transfor-
`mation. While theoretically rigorous and accurate, cal-
`culations of this type are computationally demanding
`
`0968-0896/99/$ - see front matter # 1999 Elsevier Science Ltd. All rights reserved.
`P I I : S 0 9 6 8 - 0 8 9 6 ( 9 9 ) 0 0 0 1 5 - 2
`
`

`
`852
`
`M. L. Lamb et al. / Bioorg. Med. Chem. 7 (1999) 851–860
`
`Table 1. FKBP12 inhibitorsa
`
`Compound
`
`Ki,app, nM
`
`Figure 1. The immunosuppresent drug FK506.
`
`and impractical for a large set of structurally diverse
`ligands.
`
`An attractive alternative has emerged in linear interac-
`tion energy or linear response techniques (LR), as
`recently reviewed.18 Unlike most FEP calculations,
`absolute free energies of binding ((cid:1)Gb) for protein–
`ligand systems are estimated, and simulations of non-
`physical states (intermediate steps) are not required. In
`other respects, the molecular dynamics (MD) or MC
`simulation protocols are the same.
`
`
`
`(cid:11) (cid:135) a (cid:1)ELennard(cid:255)Jones
`
`
`
`
`As first proposed, the linear interaction energy method
`employs eq (1) for the estimation of binding free ener-
`gies.19
`(cid:1)Gb (cid:136) b (cid:1)ECoulomb
`(cid:133)1(cid:134)
`These terms represent the change in interaction between
`a ligand and its environment (solvent and/or protein)
`upon binding. The electrostatic energy di€erences are
`obtained from average Coulombic interaction energies
`between the ligand and solvent (E Coulomb
`l(cid:255)water ) and the ligand
`and protein (E Coulomb
`l(cid:255)FKBP ) accumulated during simulations
`of the bound and unbound ligands. The Lennard–Jones
`(van der Waals) energy di€erences are found in an ana-
`logous manner. The original value of b=0.5 is con-
`sistent with analyses of the response of polar solutions
`to changes in electric fields, such as the charging of an
`ion in water. Observed linear correlations between the
`molecular size (surface area, chain length) and solvation
`free energies of hydrocarbons suggested van der Waals
`interactions might respond linearly as well. A scale fac-
`tor a=0.161 was obtained empirically from fitting to
`experimental binding data for a small set of endothia-
`pepsin inhibitors.19 Recently, A˚ qvist and co-workers
`have advocated the assignment of one of four reduced
`values of b according to the charge-state or polarity of
`the ligand and have obtained an appropriate value for a
`by fitting to data for a larger set of protein–ligand
`pairs.20–22 An extended linear response equation,23,24 in
`which a cavitation term based on solvent-accessible sur-
`face areas (SASA) is added and all parameters are
`empirically determined, is given in eq (2) below. This
`model was first applied to the calculation of free energies
`of hydration for organic solutes; however, optimization of
`
`(cid:11)
`
`aRotamase inhibition data taken from refs 9 and 10 unless otherwise
`noted.
`bData from ref 11.
`cKi reported for a mixture of R and S isomers.
`
`

`
`M. L. Lamb et al. / Bioorg. Med. Chem. 7 (1999) 851–860
`
`853
`
`scale factors to b=0.131, a=0.131 and g=0.014 yielded
`acceptable estimates of
`thrombin-binding a(cid:129)nity as
`well.25 In addition, the extended linear response equation
`has potential for predictions of other properties important
`for pharmacological activity, such as the estimation of
`ligand lipophilicity.24 It has been observed that the con-
`tribution of the (cid:1)SASA term is often nearly constant, so
`
`
`the case in which the simple addition of a constant
`improves the fit has also been considered [eq (3)].20,25
`h
`(cid:1)Gb (cid:136) b (cid:1)ECoulomb
`
`
`(cid:1)Gb (cid:136) b (cid:1)ECoulomb
`
`(cid:11) (cid:135) g (cid:1)SASA
`(cid:11) (cid:135) a (cid:1)ELennard(cid:255)Jones
`
`
`(cid:11) (cid:135) a (cid:1)ELennard(cid:255)Jones
`
`
`(cid:11) (cid:135) g
`
`i
`(cid:133)2(cid:134)
`
`(cid:133)3(cid:134)
`
`In contrast to most proteins studied previously with this
`technique, FKBP12 has a distinctly hydrophobic binding
`pocket lined with aromatic residues, and only two inter-
`molecular hydrogen bonds are observed in the crystal
`structure of inhibitor 4 with the protein.9 Consequently, it
`was expected that the electrostatic behavior of these mole-
`cules might deviate from linear response and that the van
`der Waals contributions to binding might be larger than
`previously observed for other systems. To determine the
`suitability of the linear response approximation for binding
`to FKBP12, MC simulations of the bound and unbound
`states for all 11 inhibitors in Table 1 were performed.
`
`Computational Details
`
`The modeling strategy employed here was consistent
`with the FEP study described previously,16 based on the
`4-FKBP12 crystal structure (1FKG)9 and the OPLS
`force field.26–28 (A full listing of parameters for these
`molecules may be found in ref 29) The first five inhibi-
`tors were easily built from 4. Atoms of the 1-phenyl
`substituent were simply removed to obtain 2, and for
`compound 1, a united-atom cyclohexyl group30 was
`positioned in the plane of the original 3-phenyl ring.
`The vinyl group of 3 was also represented with united-
`atom parameters and was positioned at the minimum of
`the CH3(cid:255)C(cid:255)CH(cid:136)CH2 torsional energy profile (180.0(cid:14))
`and aligned with an edge of the 1-phenyl ring of 4.29
`One side of this ring was within 3.2 A˚ of His87 and
`Tyr82, while the other was positioned more than 3.5 A˚
`from Phe46 and Glu54. Accordingly, the orientation
`which maximized hydrophobic contact between the
`protein and 3 was chosen. Inhibitor 5 incorporated the
`crystal structure orientations of the cyclohexyl and tert-
`pentyl groups from 5-FKBP12 (1FKH).9 In 1 and 5, the
`cyclohexyl groups were treated as rigid units, as were all
`ligand and protein aromatic groups. It was thought that
`the conformational flexibility within the rings of these
`substituents would be less important to binding than the
`overall flexibility of the ligand, and thus sampling was
`focused accordingly. Starting geometries
`for
`the
`unbound ligand simulations were taken from these
`initial bound conformations.
`
`All simulations were performed using the MCPRO
`program and Monte Carlo configurational sampling.31
`
`The ligands and protein–ligand complexes were solvated
`with 22 A˚
`spheres of TIP4P water molecules. First,
`1(cid:2)106 (1 M) configurations of water-only equilibration
`with preferential sampling of molecules close to the inhi-
`bitor was performed, followed by 16 M configurations of
`sampling of the entire system with all solvent molecules
`sampled uniformly. Next, data was collected for 8 M
`configurations averaged in blocks of 2(cid:2)105 configura-
`tions. To ensure convergence, averaging for all unbound
`ligands was extended to 16 M configurations.
`
`To obtain protein-bound structures of 6 and 7, the final
`conformations of 4 and 5 above were epimerized within
`the FKBP12 binding pocket via a slow perturbation
`protocol. Eight sequential double-wide windows with
`(cid:1)l=10.0625 were used. In each window, 4 M config-
`urations were sampled to slowly transform between the
`two ligands in an energetically reasonable way. This
`procedure was repeated with the free ligands in solution.
`The simulations of ligands 8–11 were started from the
`final
`structures of FEP simulations
`reported pre-
`viously.16,29 In each case, energy components were
`averaged over 8 M (bound) and 16 M (unbound)
`configurations of the system.
`
`As was mentioned above, the initial structure for 2 had
`been generated from the bound conformation of 4 with the
`1-phenyl atoms removed. Following an FEP calculation
`from 4 to 2, further sampling was carried out within the
`first windows of both a phenyl!pyridyl FEP16 and a car-
`bonyl!hydroxyl FEP.29 The final conformations of these
`windows were then used to start two additional linear
`response MC simulations, 2a and 2b, respectively. These
`additional data points provide one gauge of precision.
`
`Solvent-accessible surface areas for the ligands in both
`aqueous and protein environments were calculated
`using the SAVOL2 program.32 This algorithm has been
`incorporated into MCPRO, and the necessary atomic
`radii are calculated from the corresponding OPLS s
`parameters via 1/2 (21/6s). Using the standard solvent
`radius for water, 1.4 A˚ , the SASAs of the ligands were
`calculated for the structure at the end of each block of
`MC configurations and were averaged.
`
`Average energy and solvent-accessible surface area dif-
`ferences were fit
`to the experimental binding data
`(Table 1) to obtain linear response parameters a, b, and
`g according to eqs (1)–(3). This procedure was per-
`formed with a Simplex-based algorithm. As inhibition
`data for the majority of the ligands in this set has been
`reported by Holt and co-workers,9,10 in cases where two
`values have been measured the values from these
`authors were used for fitting.
`
`Results and Discussion
`
`Average intermolecular energy components and SASAs
`
`Each ligand and protein–ligand complex was solvated
`with a sphere of explicit water molecules and sampled as
`described above. Average energy components and
`
`

`
`854
`
`M. L. Lamb et al. / Bioorg. Med. Chem. 7 (1999) 851–860
`
`ligand SASAs accumulated during the simulations are
`reported in Table 2. As has been observed with throm-
`bin inhibitors, the ligand–solvent Coulombic energy
`components, E Coulomb
`l(cid:255)water , fluctuate the most during the
`simulations.25 However, both the magnitude of the
`energy contribution and the magnitude of the fluctua-
`tion are reduced compared to the positively charged or
`zwitterionic protease inhibitors, as expected for the
`neutral ligands of this set. All of the ligands in solution
`have undergone hydrophobic collapse relative to their
`bound conformations but to di€ering extents, as sig-
`nificant flexibility is observed in the propyl side chain.
`The Lennard–Jones interactions for the a-ketoamide
`molecules in solution scale generally with molecular
`size; 1 and 2 (and 8–10) have the least favorable average
`energies, followed by compound 3, and finally there is a
`small range of energies for 4–7. The solvent accessible
`surface area of amide carbonyl O3 is generally ca. 10 A˚ 2
`less than that of the keto (O4) or ester (O2) atoms, and
`usually only one water molecule is found to interact
`with this atom. This is consistent with molecular
`dynamics results for trans-FK506.29 Often, O4 interacts
`strongly with one water molecule (H-O4 distance, ca.
`1.8 A˚ ), while a second molecule hovers 0.5 A˚
`further
`away. Also, water molecules are found to interact with
`the faces of aromatic rings. In the case of 2, an unusual
`long-lived solvation pattern is noted. The water mole-
`cule that hydrogen bonds to O3 further associates with
`one directed into the center of the less-accessible face of
`the phenyl ring (Fig. 2).
`
`Typical of many FKBP12-ligand complexes,33 the crys-
`tal structure of 4-FKBP12 contains contacts between
`
`Figure 2. The hydrogen bonding network in 2, with waters that bridge
`between the amide O3 and phenyl ring.
`
`O4 and aromatic hydrogens of Tyr26, Phe36, and Phe99,
`with the pipecolyl ring sitting over Trp59.9 The 3-phe-
`nylpropyl moiety binds in the solvent-exposed FK506-
`cyclohexyl groove of FKBP12 between Ile56 and Tyr82,
`and these residues form hydrogen bonds with the ester
`and amide carbonyl oxygens of the ligand. The 1-phenyl
`substituent interacts with Phe46 and the tertiary pentyl
`group of the inhibitor. The two crystallographic inter-
`molecular hydrogen bonds are maintained throughout all
`of the FKBP12-ligand simulations, and none of the
`ligands deviates significantly from the original orientation
`
`Table 2. Average interaction energies (kcal/mol) and solvent-accessible surface areas of the inhibitors (A˚ 2) from the aqueous and FKBP12 MC
`simulationsa
`
`Compound
`
`1 aq
`1FKBP
`2 aq
`2a aq
`2b aq
`2FKBP
`2aFKBP
`2bFKBP
`3 aq
`3FBKP
`4 aq
`4FKBP
`5 aq
`5FKBP
`6 aq
`6FBKP
`7 aq
`7FKBP
`8 aq
`8FBKP
`9 aq
`9FBKP
`10 aq
`10FKBP
`11 aq
`11FKBP
`
`ECoulomb
`l(cid:255)water
`(cid:255)23.16(0.32)
`(cid:255)0.24(0.15)
`(cid:255)31.28(0.46)
`(cid:255)29.91(0.35)
`(cid:255)27.11(0.38)
`(cid:255)5.15(0.31)
`(cid:255)2.14(0.19)
`(cid:255)1.62(0.19)
`(cid:255)27.68(0.36)
`(cid:255)1.41(0.20)
`(cid:255)31.89(0.42)
`(cid:255)5.30(0.38)
`(cid:255)28.43(0.50)
`(cid:255)6.49(0.19)
`(cid:255)32.84(0.55)
`(cid:255)8.00(0.19)
`(cid:255)29.64(0.34)
`(cid:255)3.92(0.23)
`(cid:255)30.53(0.43)
`(cid:255)13.51(0.48)
`(cid:255)29.88(0.50)
`(cid:255)1.43(0.19)
`(cid:255)33.45(0.40)
`(cid:255)11.95(0.31)
`(cid:255)42.57(0.43)
`(cid:255)2.72(0.34)
`
`EL(cid:255)J
`l(cid:255)water
`(cid:255)35.96(0.15)
`(cid:255)15.09(0.12)
`(cid:255)34.59(0.21)
`(cid:255)33.90(0.16)
`(cid:255)33.26(0.16)
`(cid:255)13.74(0.09)
`(cid:255)13.99(0.12)
`(cid:255)13.95(0.14)
`(cid:255)36.52(0.16)
`(cid:255)15.31(0.13)
`(cid:255)39.63(0.22)
`(cid:255)17.84(0.11)
`(cid:255)40.12(0.19)
`(cid:255)16.54(0.10)
`(cid:255)38.35(0.12)
`(cid:255)18.10(0.13)
`(cid:255)38.48(0.15)
`(cid:255)16.10(0.15)
`(cid:255)32.75(0.18)
`(cid:255)12.03(0.15)
`(cid:255)31.67(0.17)
`(cid:255)14.04(0.09)
`(cid:255)32.89(0.22)
`(cid:255)12.82(0.10)
`(cid:255)31.24(0.21)
`(cid:255)15.19(0.10)
`
`ECoulomb
`l(cid:255)FKBP
`
`(cid:255)19.14(0.16)
`
`(cid:255)19.73(0.21)
`(cid:255)21.56(0.16)
`(cid:255)20.97(0.20)
`(cid:255)21.15(0.20)
`(cid:255)21.25(0.15)
`(cid:255)20.58(0.14)
`(cid:255)17.26(0.15)
`(cid:255)17.09(0.18)
`(cid:255)18.03(0.16)
`(cid:255)23.59(0.16)
`(cid:255)20.02(0.18)
`(cid:255)18.23(0.15)
`
`EL(cid:255)J
`l(cid:255)FKBP
`
`(cid:255)39.53(0.14)
`
`(cid:255)40.29(0.13)
`(cid:255)38.01(0.17)
`(cid:255)37.92(0.21)
`(cid:255)44.06(0.15)
`(cid:255)42.73(0.19)
`(cid:255)44.00(0.16)
`(cid:255)38.46(0.20)
`(cid:255)43.65(0.21)
`(cid:255)38.80(0.19)
`(cid:255)39.40(0.17)
`(cid:255)38.77(0.15)
`(cid:255)38.61(0.13)
`
`SASAb
`
`665.4(1.2)
`202.1(1.6)
`653.9(1.9)
`629.3(2.3)
`620.6(2.3)
`202.9(5.1)
`222.4(6.8)
`188.6(7.0)
`680.9(1.2)
`221.2(3.0)
`715.1(1.4)
`248.8(2.5)
`728.3(1.8)
`244.9(2.1)
`712.2(1.6)
`232.4(3.0)
`688.3(1.9)
`282.3(7.1)
`618.4(1.8)
`176.0(4.7)
`624.9(2.2)
`179.3(6.4)
`630.8(2.2)
`183.2(6.5)
`626.7(1.4)
`185.5(4.4)
`
`aThe standard error of the means is given in parentheses.
`bCalculated from structures saved every 2(cid:2)105 configurations, N=40 (8 M) or 80 (16 M). (Standard deviations range from 10–50 A˚ 2.)
`
`

`
`M. L. Lamb et al. / Bioorg. Med. Chem. 7 (1999) 851–860
`
`855
`
`of 4. In general, the aromatic rings of the inhibitors are
`oriented perpendicularly to the ring of Tyr82, although a
`clearly T-shaped interaction is not dominant. Consistent
`with the hydrophobic binding pocket, the largest con-
`tribution to the protein–ligand interaction energy in all
`cases comes from the Lennard–Jones terms. In addition,
`the most favorable van der Waals interactions with
`FKBP12 are observed for the more hydrophobic 3, 5,
`and 7 compared to their more aromatic or smaller
`counterparts.
`
`An estimate of the precision of the simulations was
`determined from the three simulations of 2 in solution.
`One model was derived directly from the 4-FKBP12
`crystal structure, while the others (2a and 2b) were
`obtained through earlier FEP simulations from 4. In
`solution, there is a 4 kcal/mol range in the E Coulomb
`1(cid:255)water
`components and a 2 kcal/mol range in E L(cid:255)J
`l(cid:255)water. Varia-
`tions among the average energy components reported
`for all simulations of 2-FKBP12 are on the order of 2–
`3 kcal/mol. It must be noted that the terms from simu-
`lations 2a and 2b are more similar to each other than to
`those from the original 2 simulation.
`
`As energy and SASA di€erences between the protein
`and aqueous environments are the quantities that are
`scaled to estimate binding a(cid:129)nity, these values are
`recorded in Table 3. Two of the lowest a(cid:129)nity inhibi-
`tors, 3 and 11 have the largest van der Waals energy
`di€erences between bound and unbound ligands.
`Ligands 6, 7, and 11 have the most unfavorable
`(cid:1)ECoulomb, while only 8 has a net favorable electrostatic
`interaction upon binding. Approximately 450 A˚ 2 of each
`ligand is buried upon binding to FKBP12. This repre-
`sents 65% of the surface of inhibitor 4, for example. The
`largest change is noted for inhibitors 5 and 6 (ca. 480 A˚ 2
`buried), while the smallest di€erence is found for 7 (ca.
`410 A˚ 2).
`
`Optimization of the LR equations
`
`The energies and surface areas of Table 2 were used first
`with previously determined scaling parameters to compute
`
`FKBP12-binding a(cid:129)nities. The original parameters of
`A˚ qvist19 do not describe binding to this protein well at
`all; the free energies of binding are underestimated with
`an average unsigned error (<jerrorj>) of 9.3 kcal/mol.
`Results with the thrombin-derived parameters are of
`more appropriate magnitude ((cid:255)9.3 (5) to (cid:255)6.3 (11) kcal/
`mol, <jerrorj>=1.4 kcal/mol), although both the set
`of atomic radii for SASAs and the all-atom representa-
`tion of the thrombin inhibitors di€ered from the study
`presented here. However, to improve the correspon-
`dence with experiment for FKBP12, new values for a, b,
`and g were found by fitting the average energy and sur-
`face area di€erences to the experimental binding free
`energies. The results for the various models investigated
`are summarized in Table 4.
`
`The first step was to employ eq (1) with b=0.5 and
`derive an appropriate value for a, as was done originally
`for endothiapepsin.19 This model resulted in a=0.626
`and yielded free energies which deviated significantly
`from experiment. In particular, the maximum unsigned
`error for model 1 was 4.4 kcal/mol, which diminishes its
`predictive value given that the range of experimental
`binding a(cid:129)nities is 3.4 kcal/mol. Furthermore, com-
`pound 3, a poor inhibitor, was predicted to bind as well
`as the highest a(cid:129)nity ligand, 8.
`
`As expected, the linear response assumption for elec-
`trostatic energies, b=0.5, does not appear to hold well
`for binding to FKBP12. When the value of b is set
`based on ligand composition20 and a derived empiri-
`cally with eq (1) (model 2) or both b and a treated as
`free parameters (model 3), average unsigned error and
`RMS to experiment were improved but the maximum
`errors are still greater than 2kcal/mol. With model 2, the
`a-ketoamide ligands required b=0.43; 11 called for
`b=0.37 due to its single hydroxyl substituent.20
`
`Fitting the data with three parameters yielded an RMS
`deviation of 0.7 kcal/mol, whether or not the SASA dif-
`ference was included explicitly (eq (2) versus eq (3)).
`Values of b=0.139, a=0.194, and g=0.0145 (model 4)
`ranked the ligands in a qualitatively reasonable way
`
`Table 3. Calculated energy and surface-area di€erencesa with representative binding a(cid:129)nities
`(cid:1)EL(cid:255)J
`
`(cid:1)SASA
`
`Compound
`
`(cid:1)ECoulomb
`
`1
`2
`2a
`2b
`3
`4
`5
`6
`7
`8
`9
`10
`11
`
`3.8
`6.4
`6.2
`4.5
`5.1
`5.3
`1.4
`7.6
`8.6
`(cid:255)1.0
`4.8
`1.5
`21.6
`
`(cid:255)18.7
`(cid:255)19.4
`(cid:255)18.1
`(cid:255)18.6
`(cid:255)22.7
`(cid:255)20.9
`(cid:255)20.4
`(cid:255)18.2
`(cid:255)21.3
`(cid:255)18.1
`(cid:255)21.8
`(cid:255)18.7
`(cid:255)22.6
`
`(cid:255)463.3
`(cid:255)451.0
`(cid:255)406.9
`(cid:255)432.0
`(cid:255)459.7
`(cid:255)466.3
`(cid:255)483.4
`(cid:255)479.8
`(cid:255)406.0
`(cid:255)442.4
`(cid:255)445.6
`(cid:255)447.6
`(cid:255)441.2
`
`model 2
`(cid:255)9.6
`(cid:255)8.9
`(cid:255)8.2
`(cid:255)9.2
`(cid:255)11.5
`(cid:255)10.2
`(cid:255)11.6
`(cid:255)7.6
`(cid:255)9.0
`(cid:255)11.3
`(cid:255)11.0
`(cid:255)10.6
`(cid:255)5.5
`
`model 4
`(cid:255)9.8
`(cid:255)9.4
`(cid:255)8.5
`(cid:255)9.2
`(cid:255)10.4
`(cid:255)10.1
`(cid:255)10.8
`(cid:255)9.4
`(cid:255)8.8
`(cid:255)10.1
`(cid:255)10.0
`(cid:255)9.9
`(cid:255)7.8
`
`aUnits: kcal/mol and A˚ 2, respectively.
`bReferences given in notes to Table 1, (cid:1)G=RTlnKi.
`cCompound 3 not included in the derivation of this model.
`dFree energy estimated using Ki=450 nM.
`
`(cid:1)Gb, kcal/mol
`
`model 6
`
`9.7
`(cid:255)9.4
`(cid:255)8.6
`(cid:255)9.3
`((cid:255)10.9)c
`(cid:255)10.3
`(cid:255)10.9
`(cid:255)9.0
`(cid:255)9.3
`(cid:255)10.2
`(cid:255)10.5
`(cid:255)10.0
`(cid:255)7.8
`
`exptb
`(cid:255)9.2((cid:255)9.2)
`(cid:255)9.5((cid:255)9.0)
`(cid:255)9.5((cid:255)9.0)
`(cid:255)9.5((cid:255)9.0)
`(cid:255)9.0
`(cid:255)10.9((cid:255)10.6)
`(cid:255)11.1
`(cid:255)8.7d
`(cid:255)8.9
`(cid:255)9.2((cid:255)9.4)
`(cid:255)10.1
`(cid:255)11.1
`(cid:255)7.7
`
`

`
`856
`
`M. L. Lamb et al. / Bioorg. Med. Chem. 7 (1999) 851–860
`
`Table 4. Summary of the a, b, and g parameters determined by fitting to the experimental (cid:1)Gb dataa
`max. jerrorjc
`
`gb
`
`<jerrorj>c
`
`RMS to exptc
`
`Model
`
`eq
`
`b
`
`a
`
`1
`2
`3
`4
`
`5
`6
`
`7
`
`1
`1
`1
`2
`
`3
`2
`
`3
`
`0.50d
`0.43d,e
`0.201
`0.139
`0.138(cid:139)0.016
`0.145
`0.176
`0.174(cid:139)0.013
`0.180
`
`0.626
`0.599
`0.536
`0.194
`0.191(cid:139)0.006
`0.134
`0.348
`0.348(cid:139)0.041
`0.328
`
`0.0145
`0.0146(cid:139)0.0003
`(cid:255)7.75
`0.0084
`0.0085(cid:139)0.0002
`(cid:255)4.21
`
`4.39
`2.49
`2.22
`1.39
`
`1.12
`1.09
`
`1.12
`
`1.24
`1.01
`0.71
`0.57
`
`0.55
`0.47
`
`0.47
`
`1.75
`1.25
`0.90
`0.72
`
`0.69
`0.58
`
`0.59
`
`aEnergy components and SASAs from MC simulations of 1, 2, 2a, 2b, and 3–11 were fit unless otherwise noted in the text. Cross-validation values
`with uncertainties given in italics.
`beq (2), kcal/mol.A˚ 2, eq (3), kcal/mol.
`cMaximum unsigned error, average unsigned error, and RMS deviation to experiment in kcal/mol.
`dFixed value.
`eb=0.37 for hydroxyl-containing ligand 11 (see text).
`
`(Table 3); again 3 was predicted to bind too well and
`was responsible for the largest deviation from experi-
`ment (1.4 kcal/mol). In performing the ‘‘leave-one-out’’
`cross-validation of these scaling parameters, it became
`clear that compound 3 was the notable outlier. Thus,
`the fitting was repeated without including data for this
`compound, resulting in optimal parameters b=0.176,
`a=0.348, and g=0.0084 (model 6). The predicted order
`of a(cid:129)nities for the remaining ligands improved slightly,
`and there was a corresponding decrease in RMS devia-
`tion to 0.6 kcal/mol. The average unsigned error in this
`case improves by 0.1 kcal/mol to 0.5 kcal/mol. A plot of
`calculated versus experimental binding free energies based
`on model 6 is shown in Figure 3. Additional justification
`
`for setting aside compound 3 may be found when g is a
`constant term in the linear response model. When 3 is
`included in the data set (model 5), the predicted binding
`energy results largely from the constant rather than the
`scaled energy terms (Table 4). However, when the fitting
`the value of g
`is performed without 3 (model 7),
`((cid:255)4.21 kcal/mol) accurately reflects the nearly constant
`value of g*(cid:1)SASA found with model 6.
`
`Overall, with or without 3, the LR calculations are
`ordering the binding a(cid:129)nities well. In view of the sta-
`tistical noise in the simulations and the experimental
`uncertainties, average errors of 0.5 kcal/mol are as good
`as can be hoped for.
`
`Figure 3. Calculated versus experimental (cid:1)Gb plotted for all ligands. The binding free energies were estimated using optimal parameters b=0.176,
`a=0.348, and g=0.0084 and eq (2). Compound 3, not included in this model, is not shown.
`
`

`
`M. L. Lamb et al. / Bioorg. Med. Chem. 7 (1999) 851–860
`
`857
`
`Structural observations
`
`With reasonable reproduction of the observed binding
`free energies, we now turn to structures obtained during
`the simulations for possible insights into the origins of
`the predicted a(cid:129)nities. The benchmark example in our
`earlier work was the di€erence between the 3-phenyl-
`propyl compound, 2, and the (R)-1,3-diphenylpropyl, 4.
`The reduced binding for 2 was attributed to a loss of
`hydrophobic contacts and specific aryl C(cid:255)H. . .N,O
`interactions with the protein upon removal of the 1-
`phenyl substituent.16 The relative free energy of binding
`computed with the LR method for any simulation of 2
`and 4 ((cid:1)(cid:1)Gb=0.9–1.8 kcal/mol) is in excellent agree-
`ment with the experimental and FEP-calculated values
`of 1.4 kcal/mol. Appropriately, the conformations of 2
`and 4 bound to FKBP12 resemble those reported pre-
`viously, although at the end of the simulations, the
`phenyl rings in these complexes are found nearer to
`Tyr82 than Glu54 (not shown). The mobility of the
`ligands within that region of
`the protein is seen
`throughout the MC simulations.
`
`Among the final structures of the bound simulations of
`compound 2, the two from the FEP-generated struc-
`tures (2a and 2b) are more similar to one another than
`to that of 2, especially with regard to the dicarbonyl
`O(cid:136)C(cid:255)C(cid:136)O dihedral angle, which is (cid:255)108(cid:14) in 2a and 2b
`and (cid:255)125(cid:14) in 2. However, the electrostatic energy dif-
`ferences for 2 and 2a are very similar (6.4 and 6.2 kcal/
`mol, respectively), while 2b yields 4.5 kcal/mol. There is
`also variation in the (cid:1)SASA terms for each model of 2
`that parallels the trend in contribution to binding from
`(cid:1)EL(cid:255)J for these ligands, in that 2 is most favorable, 2b
`is intermediate, and 2a is the least favorable. As a result,
`the predicted binding free energies for 2 and 2b agree
`well with the experimental values (Table 3), while 2a is
`estimated to bind 1 kcal/mol less well.
`
`As one would expect, compound 1, which contains no
`aromatic ring, has the least favorable average electro-
`static interaction of all ligands in solution ((cid:255)23.16 kcal/
`mol); however, its van der Waals interactions with sol-
`vent are comparable to those of 2. The flexible propyl
`side chain allows the cyclohexyl ring of 1 to pack
`against the side chain of Ile56 (Fig. 4). It is also observed
`that much less water is found within 3 A˚ of the cyclo-
`hexyl group than the phenyl group,
`in both bound
`(Fig. 4) and unbound structures.
`
`One of the best inhibitors, 5, is notable for its packing
`among aromatic side chains of the protein. This is
`found both in its crystal structure with FKBP129 and in
`the structures modeled here. One He of Tyr82 is 2.8 A˚
`from the ligand phenyl ring center and this protein side
`chain fills the space between the tert-pentyl and phenyl
`groups (Fig. 5). Furthermore, the pipecolyl and cyclo-
`hexyl rings are parallel to one another and separated by
`Phe46. The cyclohexyl moiety has more contact with
`Phe46 than does the analogous aromatic group of 4.
`Without these contacts in solution, the cyclohexyl ring
`is free to rotate ca. 30(cid:14) into a position still parallel to
`but displaced from the pipecolyl ring.
`
`Figure 4. Final simulated structures of 1-FKBP12 and 2-FKBP12.
`Protein residues and solvent molecules within 3 A˚ of the ligands are
`displayed.
`
`Figure 5. Snapshots from bound and unbound simulations of 5.
`
`The experimental and calculated binding a(cid:129)nities of 6
`and 7 are much reduced from their (R)- counterparts, 4
`and 5; the scaled electrostatic and van der Waals energy
`di€erences from model 6 each contribute ca. 1 kcal/mol
`to the less favorable (cid:1)Gb for these molecules. Thus, it
`appears that the bound conformations found for these
`ligands through epimerization of
`their counterparts
`(Fig. 6) are reasonable models. When bound to the
`protein, the two phenyl rings of 6 form a ‘‘parallel-
`stacked and displaced’’ structure with one another; and
`the distance between ring centers is 4.7 A˚ . This align-
`ment of aromatic rings has also been noted previously
`for benzene dimer in solution27 but is not seen in the
`simulations of 4-FKBP12. In 7-FKBP12, the Tyr82HZ-
`O3 hydrogen bond is the shortest of all complexes,
`1.7 A˚ , and the interaction between Ile56H and O2 is the
`longest at 2.6 A˚ . Unlike the 1-phenyl group of 6, the
`cyclohexyl ring in 7 interacts with Phe46 more closely
`than with its own tert-pentyl moiety.
`
`

`
`858
`
`M. L. Lamb et al. / Bioorg. Med. Chem. 7 (1999) 851–860
`
`the protein. This ‘‘bump-hole’’ design could be used to
`avoid unproductive binding of
`ligand by wild-type
`FKBP12 in the cytoplasm in a chemically-induced
`dimerization scheme.34,35 The electrostatic solute–sol-
`vent energy is 12–15 kcal/mol more favorable for the
`hydroxyl-containing 11 than its a-ketoamide parent,
`which is consistent with the preferential solvation of
`isopropanol compared with acetone36 and the ability of
`alcohols to be incorporated into the hydrogen bonding
`networks of water.37 The last structure in the unbound
`simulation has the hydroxyl hydrogen pointing towards
`ring (H(cid:255)O4(cid:255)C9(cid:255)C8 dihedral angle,
`the pipecolyl
`(cid:255)22(cid:14); H(cid:255)O3 distance, 2.4 A˚ ). It does make the expected
`two hydrogen bonds to water molecules, one as a
`donor, one as an acceptor. When bound to the protein,
`however there are no solvent molecules directly inter-
`acting with the ligand in the O4 binding pocket. The
`hydroxyl hydrogen of the ligand is positioned directly
`between and below the Hd and He of Phe36; these aro-
`matic hydrogens
`interact with O4 (Fig