Applying Linear Interaction Energy Method for Rational Design of

A training set of 20 TG analogues was used to build a binding affinity model for ... Citation data is made available by participants in Crossref's Cit...
1 downloads 0 Views 620KB Size
Download PDF
J. Med. Chem. 2005, 48, 3005-3014

3005

Applying Linear Interaction Energy Method for Rational Design of Noncompetitive Allosteric Inhibitors of the Sarco- and Endoplasmic Reticulum Calcium-ATPase Pratap Singh,† Anastasiah M. Mhaka,‡ Soren B. Christensen,§ Jeffrey J. Gray,† Samuel R. Denmeade,‡ and John T. Isaacs*,‡ Department of Chemical and Biomolecular Engineering, Johns Hopkins University, Baltimore, Maryland 21218, Anti-Cancer Drug Discovery/Development Program, The Sidney Kimmel Comprehensive Cancer Center at Johns Hopkins, Baltimore, Maryland 21231, and Department of Medicinal Chemistry, The Danish University of Pharmaceutical Sciences, 2 Universitetsparken, 2100 Copenhagen, Denmark Received August 17, 2004

Noncompetitive inhibitors of sarco- and endoplasmic reticulum calcium-ATPase (SERCA) have important therapeutic value in the treatment of cancer, due to their ability to induce apoptosis in cancer cells in a proliferation-independent manner. Thapsigargin (TG) and its analogues are one such class of inhibitors that bind to a hydrophobic pocket located in the transmembrane region of SERCA near the biomembrane surface and interfere with calcium transport. The binding free energies of thapsigargin-based inhibitors of SERCA were computed using a novel linear interaction energy (LIE) method with a surface generalized Born (SGB) continuum solvation model. A training set of 20 TG analogues was used to build a binding affinity model for estimating the free energy of binding for 18 new inhibitors with a root-mean-square (rms) error of 1.36 kcal/mol with respect to experimental data. For 15 out of the 18 inhibitors in the test set, the rms error was 1.02 kcal/mol, which is on the order of the accuracy level achieved by highly rigorous free energy of perturbation (FEP) or thermodynamic integration (TI) methods. On the basis of the analysis of the binding cavity at the interface of the membrane surface and the cytoplasmic region, we propose that side chains of TG derivatives at the O-8 position orient toward the cytoplasmic region through a hydrophobic channel. On the basis of this insight, four analogues of varying side chain length at the O-8 position with a charged moiety at the end were designed, tested with LIE methodology, and then validated experimentally for their SERCA inhibition activity. Low levels of rms error for the majority of inhibitors establish the structure-based LIE method as an efficient tool for generating more potent and specific inhibitors of SERCA by testing rationally designed lead compounds based on thapsigargin derivatization. Introduction Ca2+-ATP-

Sarcoplasmic and endoplasmic reticulum ases (SERCA) have been discovered to be important targets for activation of apoptosis in androgen-independent prostate cancer cells for which standard androgen ablation therapy has proven to be ineffective.1-3 The SERCA belongs to the family of P-type ATPases that are responsible for active transport of cations across biomembranes.4 The SERCA uses the energy released from hydrolysis of ATP to ADP for transporting calcium ions to the lumen of sarco- and endoplasmic reticulum (ER) against the electrochemical gradient. Thapsigargin (TG), obtained from the umbelliferous plant Thapsia garganica,5 is a potent inhibitor of the SERCA.6 Inhibition of SERCA by TG triggers depletion of the ER calcium pool, leading to capacitance influx of extracellular calcium into the cell. This results in * To whom correspondence should be addressed. Address: Oncology & Urology, 1650 Orleans Street, CRB 1M44, Baltimore, MD 212311000. Tel: 410-955-7777. Fax: 410-614-8397. E-mail: [email protected]. † Johns Hopkins University. ‡ The Sidney Kimmel Comprehensive Cancer Center at Johns Hopkins. § The Danish University of Pharmaceutical Sciences.

sustained elevation of calcium concentration in the cytoplasm, forcing proliferative as well as quiescent prostate cancer cells to activate a programmed cell death pathway. The SERCA is ubiquitously expressed in nearly all cell types and hence it is necessary to selectively target the cytotoxicity of TG to prostate cancer cells. To avoid general host toxicity, we have taken advantage of the fact that prostate cancer cells secret the unique serine protease prostate-specific antigen (PSA).3 Hydrolytic activity of PSA is restricted to the extracellular fluid within sites of metastatic prostate cancer. Once in the blood, PSA is inactivated by serum protease inhibitors.7 This observation has provided a rationale for developing amine-containing analogues of TG that can be incorporated into peptide prodrugs that are inactive as SERCA inhibitors, as they cannot enter cells.3 These prodrugs can be selectively activated within the sites of metastatic prostate cancer sites by PSAcatalyzed hydrolysis to release the amine-containing TG analogues which can enter into the cells and thus selectively inhibit the SERCA.3,8,9 In a number of experimental studies during the past decade, TG analogues have been designed, synthesized, and eval-

10.1021/jm049319a CCC: $30.25 © 2005 American Chemical Society Published on Web 03/16/2005

3006 Journal of Medicinal Chemistry, 2005, Vol. 48, No. 8

Singh et al.

uated.8-12 These studies provided important insights on the structure-activity relationships between analogues and SERCA inhibition. The publication of crystal structures of nucleotidefree SERCA in calcium-bound form13 (E1.2Ca2+) and in TG-bound form14 (E2′.TG) as well as in complex with a nonhydrolyzable ATP analogue (E1.AMPPCP)15,16 and with ADP stabilized by aluminum fluoride (E1.AIFx.ADP)16 has elucidated the structures of key intermediates involved in the calcium transport cycle. The availability of such structural data has facilitated the understanding of conformational changes and dynamics involved at various steps of the transport cycle. The overall structure of SERCA consists of three cytoplasmic domains and 10 transmembrane helices (Figure 1a). The two calcium binding sites in the E1 state of SERCA are about 5 Å apart and situated in the trans-membrane region, around 4 and 7 Å below the cytoplasmic surface of the membrane. The TG binding site is located in a hydrophobic cavity formed by transmembrane helices M3, M5, and M7 (Figure 1b). The polar end of the TG molecule is located near the membrane interface between residues Phe256 and Ile829. Binding of TG to SERCA is mostly hydrophobic in nature with only one hydrogen bond formed between the Ile829 backbone and the carboxyl oxygen at the O-8 position of the TG molecule. The availability of structural information on SERCA facilitates understanding the structure-activity relationships (SAR) for SERCA inhibition and enables molecular modeling techniques to be applied for designing novel and more potent inhibitors. In this study, we have applied a newly developed structure-based linear interaction energy method implementing a surface generalized Born (SGB-LIE)17 continuum model for solvation to build a binding affinity model for estimating the free energy of binding for a diverse set of SERCA inhibitors. The LIE method18,19 has been applied on a number of protein-ligand systems with promising results,20-29 producing small errors on the order of 1 kcal/mol for free energy prediction.27,30 The magnitude of free energy changes upon binding of inhibitors to SERCA directly correlates with the experimental potency of these inhibitors; hence, fast and accurate estimation of binding free energies provides a means to screen the compound libraries for lead optimization and rational design. Materials and Methods LIE Methodology. The LIE method employs experimental data on binding free energy values for a set of ligands (referred as training set) to estimate the binding affinities for a set of novel compounds. The method is based on the linear response approximation (LRA), which dictates that binding free energy of a protein-ligand system is a function of polar and nonpolar energy components that scale linearly with the electrostatic and van der Waals interactions between a ligand and its environment. The free energy of binding for the complex is derived from considering only two states: (1) free ligand in the solvent and (2) ligand bound to the solvated protein. The conformational changes and entropic effects pertaining to unbound receptor are taken into account implicitly and only interactions between the ligand and either the protein or solvent are computed during molecular mechanics calculations. Among the various formulations of the LIE methodology developed in the past, the SGB-LIE method implementing a

Figure 1. Thapsigargin bound to SERCA. (a) The thapsigargin molecule (shown in blue and yellow spheres) binds to SERCA in a cavity formed by residues on transmembrane helices 3, 5, and 7. The ATP binding site is located near residue F487 (yellow asterisk) in the nucleotide domain (N) in purple. The phosphorylating residue (Asp351) is situated near a location marked by an orange asterisk in the phosphorylation (P) domain, which is shown in green. (b) Close up view of TG binding site formed by mostly hydrophobic residues with only one hydrogen bond between Ile829 and the oxygen atom on the O-8 side chain. surface generalized Born (SGB) model for the solvation has been shown to be 1 order of magnitude faster than the methods based on explicit solvent17 with the same order of accuracy. The SGB-LIE method also offers better accuracy in treating the long-range electrostatic interactions. The SBG-LIE method implements the original formulation proposed by Jorgensen31

Noncompetitive Inhibitors of SERCA

Journal of Medicinal Chemistry, 2005, Vol. 48, No. 8 3007

Table 1. First Subset of Inhibitors of SERCA Comprising the Training Set

for the case of continuum solvent replacing the solvent accessible surface area term by a cavity term as follows

∆G ) R(〈Uvdwb〉 - 〈Uvdwf〉) + β(〈Ueleb〉 - 〈Uelef〉) + γ(〈Ucavb〉 - 〈Ucavf〉) (1) where bracketed terms represents the ensemble average of the energy terms, such as van der Waals (Uvdw), electrostatic (Uele), or cavity (Ucav) energy. The energy terms involved can be computed using energy minimization, molecular dynamics, or Monte Carlo calculations. All the terms are evaluated for interaction between ligand, both in the free (f) and bound (b) state, and its environment. The R, β, and γ are LIE fitting parameters. The transferability and dependence of LIE parameters on force fields and protein-ligand system are still the subject of debate. In the Jorgensen formulation, LIE parameters are free coefficients that need to be determined by fitting the experimental data on the training set compounds. In the SGB model of solvation, there is no explicit van der Waals or electrostatic interaction between the solute and

solvent. The contribution for net free energy of solvation comes from two energy terms, namely, reaction field energy (Urxn) and cavity energy (Ucav):

USGB ) Urxn + Ucav

(2)

The cavity and reaction field energy terms implicitly take into account the van der Waals and the electrostatic interactions, respectively, between the ligand and solvent. The application of the SGB-LIE method for a given protein-ligand system essentially involves computing four energy components, i.e., the van der Waals and Coulombic energy between the ligand and protein and the reaction field and cavity energy between the ligand and continuum solvent. The total electrostatic energy in the SGB-LIE method is the sum of Coulombic and reaction field energy terms. Computational Details. Preparation of receptor and ligands was done using the MOE package supplied by Chemical Computing Group, Montreal, Canada. All the calculations for the SGB-LIE method were performed in the Liaison package from Schrodinger Inc.32 The Liaison module performs

3008

Journal of Medicinal Chemistry, 2005, Vol. 48, No. 8

Singh et al.

Table 2. Second Subset of Inhibitors of SERCA Comprising the Training Set

LIE calculations in the OPLS force field with a residue-based cutoff of 15 Å. The OPLS force field was also used for charge assignment and all energy calculations performed in the MOE and Liaison package. Protein Preparation. Initial coordinates for the SERCA bound with TG were taken from the RSCB Protein Data Bank (PDB ID 1IWO), which describes the structure of protein in dimer form. The coordinates for only one monomer bound with TG were used for all calculations in this study. To speed up the protein-ligand energy calculations, protein was truncated in such a way that it contained only the transmembrane domain, lumen, and part of the cytoplasmic domain near the membrane surface where the binding site for TG is located. This truncation is not expected to affect the accuracy, as residues located faraway from the binding site do not affect the protein-ligand interactions. All residues, except histidines, were kept at their default protonation states at physiological pH. The protonation states of the histidine residues were determined in such a way that chosen protonation states

favored the formation of hydrogen bonds with nearby atoms. In the cases where the choice of protonation state was not obvious, histidine residues were protonated at the δ-nitrogen. All the histidines are very far from the binding site and there are no charged residues interacting with the inhibitors in the vicinity of the hydrophobic binding site. The truncated protein contained 445 residues. The net charge on the protein was -16, but under SGB-LIE implementation, there is no need for neutralizing the protein. After the addition of hydrogen atoms, the complex structure was energy minimized in the MOE package while all atoms except hydrogen were kept fixed. The motivation behind keeping protein atoms fixed was to preserve the structure of the binding site during energy calculations. The energy-minimized SERCA structure was subsequently used for docking of TG analogues and SGB-LIE calculations. Preparation and Docking of Ligands. All the ligand molecules were built using the Builder feature in MOE and energy minimized in a vacuum. The amine moiety present in the analogues was protonated while the carboxyl group was

Noncompetitive Inhibitors of SERCA

Journal of Medicinal Chemistry, 2005, Vol. 48, No. 8 3009

Table 3. Average van der Waals (vdw), Electrostatic (ele), Reaction Field (rxn), and Cavity (cav) Energy Terms for the Free (f) and Bound State (b) Obtained from SGB-LIE Calculations ligand

Uvdwf,a kcal/mol

Uelef,a kcal/mol

Urxnf, kcal/mol

Ucavf, kcal/mol

Uvdwb, kcal/mol

Ueleb, kcal/mol

Urxnb, kcal/mol

Ucavb, kcal/mol

TG TR1 TR2 TR3 TR4 TR5 TR6 TR7 TR8 TR9 TR10 TR11 TR12 TR13 TR14 TR15 TR16 TR17 TR18 TR19

0.0 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

-28.9 -28.5 -27.5 -27.2 -26.4 -21.4 -29.0 -22.2 -28.9 -24.2 -27.0 -28.5 -28.6 -26.5 -29.7 -25.5 -31.4 -31.1 -29.2 -29.7

6.01 5.75 5.98 6.04 6.13 5.67 6.29 5.91 5.92 5.67 4.76 6.03 6.17 5.99 6.57 6.09 5.78 5.41 6.01 5.94

-73.4 -69.0 -70.2 -62.1 -75.8 -69.0 -69.5 -65.6 -54.6 -59.9 -51.9 -72.2 -76.5 -76.5 -53.6 -76.0 -67.4 -58.2 -54.3 -16.4

-12.26 -12.6 -6.8 -12.9 1.7 -6.3 -20.3 -8.4 -8.6 -7.8 -10.7 -10.3 -9.6 -6.5 -19.9 -1.3 -10.2 -7.5 -7.5 -5.7

-27.5 -27.9 -6.0 -28.5 -32.0 -10.9 -23.6 -13.5 -23.6 -14.9 -28.5 -30.6 -27.9 -14.8 272.4 -15.3 -30.4 -24.5 -25.0 -16.6

1.82 1.77 2.05 1.86 1.72 1.74 1.84 1.75 1.81 1.76 1.49 1.84 1.79 1.95 1.92 1.96 1.84 1.86 2.52 1.89

a The van der Waals and electrostatic energy terms for the free state are zero for all ligands, since there is no explicit interaction between ligand and solvent in continuum solvent.

Table 4. Binding Affinity Model Calculations for Training Set Inhibitors Using the SGB-LIE Methoda ligand

RIC50

∆Grel, kcal/mol

∆Gabs, kcal/mol

∆GLIE, kcal/mol

∆Gjk, kcal/mol

TG TR1 TR2 TR3 TR4 TR5 TR6 TR7 TR8 TR9 TR10 TR11 TR12 TR13 TR14 TR15 TR16 TR17 TR18 TR19

1.0 0.98 0.85 0.52 0.40 0.186 0.08 0.06 0.02 0.002 0.0002 0.67 0.40 0.36 0.09 0.07 0.02 0.015 0.002 0.0003

0.0 0.0116 0.0918 0.3815 0.5382 0.9852 1.4277 1.6360 2.1632 3.7234 5.1175 0.2371 0.5359 0.6022 1.4024 1.5838 2.1860 2.4503 3.5572 4.7062

-11.20 -11.19 -11.11 -10.82 -10.66 -10.21 -9.72 -9.56 -9.04 -7.48 -6.08 -10.96 -10.66 -10.60 -9.80 -9.61 -9.01 -8.75 -7.64 -6.49

-10.41 -9.828 -9.875 -9.583 -10.855 -9.785 -10.440 -9.882 -8.963 -9.11 -7.720 -10.319 -10.862 -10.447 -9.818 -10.53 -9.671 -8.533 -8.179 -6.162

-10.34 -9.72 -9.79 -9.49 -10.87 -9.78 -10.49 -9.87 -8.96 -9.18 -7.78 -10.22 -10.87 -10.43 -11.30 -10.60 -9.74 -8.61 -8.18 -4.74

a

RIC50 refers to the relative potency of a SERCA inhibitor with respect to TG and is defined as the ratio of the IC50 value of TG (TGIC50) and a SERCA inhibitor. The experimental IC50 value of TG (6.17 nM) was taken to be the average of four different measurements conducted in refs 8-11. The relative and absolute free energy values for SERCA inhibition (∆Grel and ∆Gabs, respectively) are computed using ∆Grel ∼ RT ln(RIC50) and ∆Gabs ∼ RT ln(TGIC50/RIC50). ∆GLIE and ∆Gjk refer to the absolute free energy values obtained using the SGB-LIE method and jackknife test, respectively.

deprotonated. The partial charges on the atoms were assigned using the OPLS force field. The binding mode of all the TGbased compounds is likely to be similar to that of TG, as they have very similar structures. The hydrophobic channel near the membrane surface was assumed to be the docking position for the long side chains at the O-8 position. Each ligand was superimposed with SERCA-bound TG to produce the starting configurations. Energy minimization for each analogue bound to SERCA was performed in the vacuum while all the protein atoms other than nonpolar hydrogen atoms were kept fixed. The energy-minimized complex corresponding to each analogue was transported to the LIAISON package for subsequent SGBLIE calculations.

Chart 1. Free Energy Values Estimated by the SGB-LIE Method for 20 Thapsigargin Analogues Comprising the Training Set Plotted against the Corresponding Experimental Dataa

a The rms and correlation coefficient (R2) between the two data sets are 0.834 kcal/mol and 0.71, respectively. The solid straight line represents the perfect match between experimental and predicted values.

LIE Calculations. To compute the interaction terms required in the LIE model, energy minimization was performed in the LIAISON package, and various energy terms were collected at the end of the minimization. The conjugated gradient technique was used for the minimization with a maximum of 500 steps. The root mean square (rms) gradient value of 0.01 was the criteria for the convergence of minimization. Theoretically, molecular dynamics (MD) or Monte Carlo (MC) calculations are more efficient in sampling the ensemble average values for energy terms, but sampling using MD or MC methods is computationally expensive and does not improve the accuracy by a significant amount when the initial docked positions of the compounds are close to the actual binding mode.27 The various energy terms collected by LIAISON after the energy minimization were used for building binding affinity model and free energy estimation for TG analogues. Experimental Section. All the TG compounds were synthesized as described previously.8-12,33 Measurements of IC50

3010

Journal of Medicinal Chemistry, 2005, Vol. 48, No. 8

Table 5. Inhibitors of SERCA Comprising the Test Set

Singh et al. Table 6. Estimation of Binding Free Energies (∆G, kcal/mol) for the TG Analogues Comprising the Test Seta

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

compd

RIC50

∆Grel, kcal/mol

∆Gabs, kcal/mol

∆GLIE, kcal/mol

dbTG 4a 4c 6a 5c 6c 5d 6d 6e 6f 6g 6h 6i 6j 7k 7l 7m 7n

0.190 0.220 0.003 0.590 0.520 0.680 0.540 0.590 0.600 0.01 0.011 0.059 0.333 0.384 0.004 0.294 0.833 1.250

0.971 0.886 3.398 0.308 0.382 0.226 0.360 0.309 0.299 2.722 2.665 1.678 0.650 0.565 3.352 0.724 0.108 -0.132

-10.229 -10.315 -7.803 -10.891 -10.818 -10.974 -10.840 -10.891 -10.902 -8.478 -8.535 -9.522 -10.550 -10.634 -7.848 -10.475 -11.092 -11.332

-9.005 -11.668 -8.869 -9.714 -11.474 -9.334 -11.513 -9.615 -8.936 -10.057 -9.326 -8.971 -9.539 -9.854 -11.258 -10.182 -10.202 -10.208

a The R IC50 refers to the relative potency of a SERCA inhibitor with respect to TG and is defined as the ratio of the IC50 value of TG and a SERCA inhibitor. The relative and absolute experimental free energy values for SERCA inhibition are denoted by ∆Grel and ∆Gabs, respectively. The absolute free energy values obtained using the LIE methodology is denoted by ∆GLIE.

values for SERCA inhibition for each compound were conducted as described previously.10

Results and Discussions We have applied the SGB-LIE method to a training set of 20 inhibitors of SERCA to build a binding affinity model that was then used to compute the free energy of binding for a test set of 18 new TG analogues as well as for a design set of four newly synthesized TG analogues. The training set for building the binding affinity model was comprised of two subsets of TG analogues. The first subset (Table 1) was comprised of 11 ligands designed and synthesized by Nielsen et al.,11 whereas the second set (Table 2) was comprised of nine ligands designed by Christensen et al.10 Analogues in the first subset had been synthesized by alkylating or acetylating O-11 and O-12 position in the lactol obtained by reducing thapsigargicin. Analogues comprising the second

subset are more diverse in their structure and were synthesized by derivatizing different locations in the key lactone ring of the thapsigargin molecule. For both subsets, experimental IC50 values and the relative potencies (RIC50, defined as the ratio of the IC50 value of TG with respect to the IC50 value of the analogue) for the inhibition of microsomal Ca-ATPase are available. With the 4 orders of magnitude difference between the IC50 values and the large diversity in the structures, the combined set of 20 ligands is ideal to be considered as a training set, as the set does not suffer from bias, due to the similarity of the structures. Also, the training set containing 20 analogues contains enough data points not to suffer from overparametrization by the LIE model. Training set compounds were docked into the TG binding site, and the SGB-LIE calculations were performed using the LIAISON module. The simulations were performed both for the ligand-free and ligandbound state. The various interaction energy terms described in the methods were collected by the LIAISON package and are presented here in Table 3. The largest contribution for the binding energy comes from the van der Waals (VDW) interactions, which is expected as TG and its analogues are mostly lipophilic molecules that interact favorably with a binding cavity lined with hydrophobic residues. The cavity energy term in the bound state is smaller than in the free state for all compounds, as there is less energy penalty for creating a cavity in solvent when part of the ligand is buried into the hydrophobic binding site. The reaction field energy term in the free state lies in a very narrow range for all compounds, but it varies in a wide range in the bound state as the solvent accessible surface area varies with ligand structure in the bound form. The energy values in Table 3 were used to fit eq 1 using the least-square error-fitting method. The values obtained for the three fitting parameters, R, β, and γ, are 0.072, -0.0004, and 1.228, respectively. The large value of the cavity energy term signifies the fact that binding is largely driven by

Noncompetitive Inhibitors of SERCA

Journal of Medicinal Chemistry, 2005, Vol. 48, No. 8 3011

Chart 2. Free Energy Values Estimated by the SGB-LIE Method for 18 Thapsigargin Analogues Comprising the Test Set Plotted against the Corresponding Experimental Dataa

a

The rms error for the whole set was 1.36 kcal/mol. For 15 analogues, the rms error is 1.02 kcal/mol.

the ligand’s ability to bury itself in the binding cavity, which is understandable given that most of the ligands are highly hydrophobic in nature. Even though the R value is low, VDW interactions contribute significantly toward the free energy of binding due to the large magnitude of the VDW interaction term. In Table 4, the experimental free energy values obtained from the RIC50 values and the free energy values estimated using fitting parameters are presented. The rms error between the experimental values and the values obtained by the fit was 0.834 kcal/mol, which is an indicator of the robustness of the fit. The quality of the fit can also be judged by the value of the squared correlation coefficient (R2), which was 0.71 for the training set. It is interesting to note that the estimated free energy of binding for ligand TR19, a stereoisomer of TG, comes out to be very close to the experimental value. This indicates that the SGB-LIE fit is sensitive enough to enable the distinction of the stereoisomers based on their inhibitory potencies. Chart 1 graphically shows the quality of fit. Table 4 also presents the jackknife cross-validation results. The jackknife test was performed in order to estimate the bias and prediction power of the binding affinity model developed using the training set. The jackknife test involves leaving one compound out of the training set and then using fitting parameters obtained from the rest of the training set to predict the binding free energy of the ligand left out. The rms error between the experimental and the predicted free energy values for the jackknife test is 1.02 kcal/mol, which suggests a reasonable prediction capability of the SGB-LIE method for estimation of free energy values. The biggest contribution for error is coming from the ligands TR9, TR10, and TR19. The jackknife test produces overestimated (more negative than the experimental value) values for free energy of binding for TR9 and TR10, but the opposite is true for ligand TR19. The large value of error (∼1.7 kcal/mol) for TR19 is linked to the fact that the jackknife test produces quite different values of fitted LIE parameters (R ) 0.1, β ) -0.0015, and γ ) 0.76) when TR19 is left out. Interestingly, TR19 is the only compound in the training set that has different

stereochemistry at the O-8 position. By including this compound in the training set, we train the LIE to make a close prediction for similar analogues having different stereochemistry at the O-8 position. When all three outliers were left out from the training set, the values of the LIE parameters did not differ much (R ) 0.076, β ) -0.0002, and γ ) 1.19), but the quality of fit is improved (rms ) 0.687 kcal/mol) as expected from the fact that these three compounds were the source of largest error in the training set. Satisfied with the robustness of the binding affinity model developed using the training set, we applied the LIE model to the TG analogues comprising the test set. The analogues comprising the test set (Table 5) were synthesized by Christensen et al.8 and by Jakobsen et al.9 by replacing the butanoyl residue at the O-8 position with the residue containing aromatic amines or a side chain containing an amino acid suitable for conjugating to a peptide substrate. The main aim of these synthesis studies was to develop a drug moiety to be incorporated into a prodrug for the treatment of prostate cancer. Since the experimental values of IC50 for these inhibitors are already available, this set of molecules provides an excellent data set for testing the prediction power of the SGB-LIE method for new ligands. Table 6 presents the free energy values estimated for the 18 test compounds for which experimental IC50 values were available to enable the accuracy check. The overall rms error between the experimental and predicted free energy values was 1.36 kcal/mol. This rms error is reasonably low in comparison with the results of previous LIE studies.27,34 For compounds 6e, 6f, and 7k, the error is more than 1.36 kcal/mol. The rms error for the rest of the 15 compounds is 1.02 kcal/mol, which means that the LIE modeling was able to predict the free energy of the 15 compounds within 1.02 kcal/mol, which is comparable to the level of accuracy achieved by the most accurate methods, such as free energy perturbation. The squared correlation coefficient for the test set compounds is low (R2 ) 0.22), but the level of rms error (1.36 and 1.02 kcal/mol, with and without the three outliers, respectively) is low enough to discriminate milli-, micro-, and nanomolar inhibitors.

3012

Journal of Medicinal Chemistry, 2005, Vol. 48, No. 8

Singh et al.

Figure 2. The channel formed by an opening between the M3 and M5 helices near the O-8 position of TG. (a) The front and side view of the thapsigargin binding site depicting the location of the channel that makes possible for the long side chains at the O-8 position in thapsigargin analogues to get exposed to the solvent on the cytoplasmic side (side view). (b) Close-up view depicting residues on M3 and M5 helices lining the channel.

The SGB-LIE model estimates for the free energy values for the analogues 7l, 7m, and 7n are quite close to the experimental one. These analogues were synthesized by substituting the O-8 butyric acid with an 11carbon-long alkyl chain conjugated with the R-amino acid at the end. The LIE results confirm with the experimental evidence that these compounds have SERCA inhibition potencies very close to TG. The estimated free energy values for the test set ligands are plotted against the experimental data in Chart 2. Apart from three poorly estimated values, there is a close match between the experimental and LIE free energy values for the rest of the ligands in the test set. The poor predictions for 6e, 6f, and 7k may be due to the novel binding modes for these compounds, which are different from TG binding. Investigation of the TG binding site in the SERCA pump reveals that the O-8 position of TG is aligned toward a tight hydrophobic channel near the biomembrane interface formed by an opening between the M3 and M5 helices. Capped by loop L67 (the loop between

M6 and M7) on the front end, the channel leads to a solvent-accessible region near Glu309 (Figure 2). We propose this hydrophobic channel to be the docking position of side chains at O-8 position of TG analogues. The side chains at the O-8 position of ligands 6e, 6f, and 7k are of intermediate in size with a charged amine at the end. Due to insufficient length, side chains at the O-8 position seem to be unable to reach out to the solvent through the hydrophobic channel. Unfavorable interactions of charged amine with the hydrophobic residues lining the channel make it likely for these ligands to adopt a binding orientation different than the one we have used for LIE calculations. The proposed binding orientation of the long side chains at O-8 through the hydrophobic channel is consistent with the experimental observation that long derivatization at the O-8 position is possible without significant loss of inhibitory activity. Interestingly, even charged groups can be placed at the end of the side chain, provided the hydrophobic portion of the side chain is long enough to expose the charged moieties to the

Noncompetitive Inhibitors of SERCA

Journal of Medicinal Chemistry, 2005, Vol. 48, No. 8 3013

Table 7. Free Energy Values (∆G, kcal/mol) and Relative Potencies (RIC50), Obtained from the SGB-LIE Method and Experimental Data for Design Set Inhibitors

high level of accuracy for a range of compounds with varying inhibition potencies. Despite the limitation imposed by the insufficient sampling inherent in the energy minimization protocol, the method has reproduced experimental data with reasonably small error for the majority of TG analogues. Using LIE methodology, we have been able to verify the experimental observation that derivatized thapsigargin compounds, with their O-8 butyric acid side chain replaced with R-amino acid conjugated substrates, have inhibition potencies close to that of thapsigargin. Moreover, the SGB-LIE method is able to predict the relative potencies of rationally designed ligands with relative success. On the basis of careful observation of binding cavity, we propose a binding model for TG analogues derivatized at the O-8 position in which long side chains at the O-8 position dock toward the hydrophobic channel in such a way that the polar groups at the side chain ends are exposed to solvent near the cytoplasmic surface. A detailed study on the structure-activity relationships for SERCA inhibitors can throw light on the moieties and functional groups important in determining the inhibition potency. This study is currently being undertaken in our lab. We are also studying the effect of sampling protocols, i.e., energy minimization, molecular dynamics, hybrid Monte Carlo methods, on the quality of predictions. The close estimation of inhibition potencies of a wide range of compounds establishes the LIE methodology as an efficient tool for screening novel compounds with very different structures. Compared to the empirical methods, such as scoring function approaches, the LIE method is more accurate due to the semiempirical approach adopted in which experimental data are used to build the binding affinity model. The SGB-LIE method seems promising when compared to the FEP or TI methods in achieving comparable accuracy with must faster speed even for structurally very different ligands.

solvent on the cytoplasmic side. This is in contrast to the fact that most of the interactions between TG and SERCA are hydrophobic in nature. To explore the effect of side chain length and position of the charged moiety on inhibitory activity, four novel compounds with different side chain length at the O-8 position were synthesized (Table 7). All four compounds carry a charged moiety at the end of the side chain. The SGBLIE methodology was used to estimate the free energy of binding and SERCA inhibitory potencies for these compounds. To judge the accuracy of the SGB-LIE estimation of inhibitory potencies, we synthesized and tested these compounds in our lab for SERCA inhibitory potency using an ATPase assay described in the method section. Table 7 presents the free energy values and corresponding IC50 values for the four newly designed ligands. Experimentally determined relative potencies of the drugs are also provided in order to evaluate the accuracy of LIE predictions. For all compounds, LIE predictions produce exactly the same trend for relative potencies, even though the exact magnitudes of these values do not match very well. The difference in the exact magnitudes of estimated vs experimental free energy of binding for compounds in the design set and a few ligands in the test set may be a contribution of two factors in addition to the limitations imposed by inadequate sampling and force field parametrization. The first has to do with the fact that calculation of absolute free energy of binding from experimental IC50 values is only an approximation due to variation in experimental conditions. The second factor contributing to error is the limitation of the continuum solvent model for SERCA inhibitors. Even though the TG binding cavity is located on the interface of membrane and cytoplasm, there are subtle interactions between alkyl side chains of TG analogues and lipid molecules that have not been accounted for in our SGB-LIE calculations. Conclusions We have demonstrated that the SGB-LIE method can be applied to estimate the free energy of binding with a

Acknowledgment. This work was supported by NIH SPORE grant P50CA58236. References (1) Furuya, Y.; Lundmo, P.; Short, A. D.; Gill, D. L.; Isaacs, J. T. The Role of Calcium, pH, and Cell-Proliferation in the Programmed (Apoptotic) Death of Androgen-Independent ProstaticCancer Cells Induced by Thapsigargin. Cancer Res. 1994, 54, 6167-6175. (2) Tombal, B.; Denmeade, S. R.; Gillis, J. M.; Isaacs, J. T. A supramicromolar elevation of intracellular free calcium ([Ca2+](i)) is consistently required to induce the execution phase of apoptosis. Cell Death Differentiation 2002, 9, 561-573. (3) Denmeade, S. R.; Jakobsen, C. M.; Janssen, S.; Khan, S. R.; Garrett, E. S.; Lilja, H.; Christensen, S. B.; Isaacs, J. T. Prostate Specific Antigen-activated thapsigargin prodrug as targeted therapy for prostate cancer. J. Natl. Cancer Inst. 2003, 95, 9901000. (4) Kuhlbrandt, W. Biology, structure and mechanism of P-type ATPases. Nat. Rev. Mol. Cell Biol. 2004, 5, 282-295. (5) Rasmussen, U.; Christensen, S. B. Thapsigargine and Thapsigargicine, 2 New Histamine Liberators from Thapsia-Garganica L. Acta Pharm. Suecica 1978, 15, 133-140. (6) Lytton, J.; Westlin, M.; Hanley, M. R. Thapsigargin Inhibits the Sarcoplasmic or Endoplasmic-Reticulum Ca-ATPase Family of Calcium Pumps. J. Biol. Chem. 1991, 266, 17067-17071. (7) Christensson, A.; Laurell, C. B.; Lilja, H. Enzymatic-Activity of Prostate-Specific Antigen and Its Reactions with Extracellular Serine Proteinase-Inhibitors. Eur. J. Biochem. 1990, 29, 755763. (8) Christensen, S. B.; Andersen, A.; Kromann, H.; Treiman, M.; Tombal, B.; Denmeade, S.; Isaacs, J. T. Thapsigargin analogues for targeting programmed death of androgen-independent prostate cancer cells. Bioorg. Med. Chem. 1999, 7, 1273-1280.

3014

Journal of Medicinal Chemistry, 2005, Vol. 48, No. 8

(9) Jakobsen, C. M.; Denmeade, S. R.; Isaacs, J. T.; Gady, A.; Olsen, C. E.; Christensen, S. B. Design, synthesis, and pharmacological evaluation of thapsigargin analogues for targeting apoptosis to prostatic cancer cells. J. Med. Chem. 2001, 44, 4696-4703. (10) Christensen, S. B.; Andersen, A.; Poulsen, J. C. J.; Treiman, M. Derivatives of Thapsigargin as Probes of Its Binding-Site on Endoplasmic-Reticulum Ca2+ ATPasesStereoselectivity and Important Functional-Groups. FEBS Lett. 1993, 335, 345-348. (11) Nielsen, S. F.; Thastrup, O.; Pedersen, R.; Olsen, C. E.; Christensen, S. B. Structure-Activity-Relationships of Analogues of Thapsigargin Modified at O-11 and O-12. J. Med. Chem. 1995, 38, 272-276. (12) Carsten M. Jakobsen, A. M., John T. Isaacs, Soren Brogger Christensen, Carl E. Olsen, and Samuel R. Denmeade Thapsigargin Analogues Applicable for Prodrugs Activated by ProstateSpecific Membrane Antigen. In press 2004. (13) Toyoshima, C.; Nakasako, M.; Nomura, H.; Ogawa, H. Crystal structure of the calcium pump of sarcoplasmic reticulum at 2.6 angstrom resolution. Nature 2000, 405, 647-655. (14) Toyoshima, C.; Nomura, H. Structural changes in the calcium pump accompanying the dissociation of calcium. Nature 2002, 418, 605-611. (15) Toyoshima, C.; Mizutani, T. Crystal structure of the calcium pump with a bound ATP analogue. Nature 2004, 430, 529-535. (16) Sorensen, T. L. M.; Moller, J. V.; Nissen, P. Phosphoryl transfer and calcium ion occlusion in the calcium pump. Science 2004, 304, 1672-1675. (17) Zhou, R. H.; Friesner, R. A.; Ghosh, A.; Rizzo, R. C.; Jorgensen, W. L.; Levy, R. M. New linear interaction method for binding affinity calculations using a continuum solvent model. J. Phys. Chem. B 2001, 105, 10388-10397. (18) Aqvist, J.; Medina, C.; Samuelsson, J. E. New Method for Predicting Binding-Affinity in Computer-Aided Drug Design. Protein Eng. 1994, 7, 385-391. (19) Aqvist, J.; Marelius, J. The linear interaction energy method for predicting ligand binding free energies. Combin. Chem. High Throughput Screening 2001, 4, 613-626. (20) Tominaga, Y.; Jorgensen, W. L. General model for estimation of the inhibition of protein kinases using Monte Carlo simulations. J. Med. Chem. 2004, 47, 2534-2549. (21) Leiros, H. K. S.; Brandsdal, B. O.; Andersen, O. A.; Os, V.; Leiros, I.; Helland, R.; Otlewski, J.; Willassen, N. P.; Smalas, A. O. Trypsin specificity as elucidated by LIE calculations, X-ray structures, and association constant measurements. Protein Sci. 2004, 13, 1056-1070. (22) Ostrovsky, D.; Udier-Blagovic, M.; Jorgensen, W. L. Analyses of activity for factor Xa inhibitors based on Monte Carlo simulations. J. Med. Chem. 2003, 46, 5691-5699.

Singh et al. (23) Tounge, B. A.; Reynolds, C. H. Calculation of the binding affinity of beta-secretase inhibitors using the linear interaction energy method. J. Med. Chem. 2003, 46, 2074-2082. (24) Wesolowski, S. S.; Jorgensen, W. L. Estimation of binding affinities for celecoxib analogues with COX-2 via Monte CarloExtended linear response. Bioorg. Med. Chem. Lett. 2002, 12, 267-270. (25) Chen, J. G.; Wang, R. X.; Taussig, M.; Houk, K. N. Quantitative calculations of antibody-antigen binding: Steroid-DB3 binding energies by the linear interaction energy method. J. Org. Chem. 2001, 66, 3021-3026. (26) Pierce, A. C.; Jorgensen, W. L. Estimation of binding affinities for selective thrombin inhibitors via Monte Carlo simulations. J. Med. Chem. 2001, 44, 1043-1050. (27) Wall, I. D.; Leach, A. R.; Salt, D. W.; Ford, M. G.; Essex, J. W. Binding constants of neuraminidase inhibitors: An investigation of the linear interaction energy method. J. Med. Chem. 1999, 42, 5142-5152. (28) Lamb, M. L.; Tirado-Rives, J.; Jorgensen, W. L. Estimation of the binding affinities of FKBP12 inhibitors using a linear response method. Bioorg. Med. Chem. 1999, 7, 851-860. (29) Smith, R. H.; Jorgensen, W. L.; Tirado-Rives, J.; Lamb, M. L.; Janssen, P. A. J.; Michejda, C. J.; Smith, M. B. K. Prediction of binding affinities for TIBO inhibitors of HIV-1 reverse transcriptase using Monte Carlo simulations in a linear response method. J. Med. Chem. 1998, 41, 5272-5286. (30) van Lipzig, M. M. H.; Vermeulen, N. P. E.; Wamelink, M.; Geerke, D.; Jongejan, A.; ter Laak, T. M.; Meerman, J. H. N. Molecular dynamics simulations with the estrogen receptor: Prediction of ligand binding affinity and orientation by the linear interaction energy method. Drug Metab. Rev. 2003, 35, 107107. (31) Carlson, H. A.; Jorgensen, W. L. An Extended Linear-Response Method for Determining Free-Energies of Hydration. J. Phys. Chem. 1995, 99, 10667-10673. (32) FirstDiscovery2.7, 2.7 ed.; Schrodinger Inc.: Portland, 2004. (33) Andersen, A.; Cornett, C.; Lauridsen, A.; Olsen, C. E.; Christensen, S. B. Selective Transformations of the Ca-2+ Pump Inhibitor Thapsigargin. Acta Chem. Scand. 1994, 48, 340-346. (34) Hou, T. J.; Zhang, W.; Xu, X. J. Binding affinities for a series of selective inhibitors of gelatinase-A using molecular dynamics with a linear interaction energy approach. J. Phys. Chem. B 2001, 105, 5304-5315.

JM049319A