DASH: A Novel Analysis Method for Molecular Dynamics Simulation

Mar 18, 2005 - A novel molecular dynamics (MD) analysis algorithm, DASH, is introduced in this paper. DASH has been developed to utilize the sequentia...
7 downloads 19 Views 356KB Size
Download PDF
3214

J. Med. Chem. 2005, 48, 3214-3220

DASH: A Novel Analysis Method for Molecular Dynamics Simulation Data. Analysis of Ligands of PPAR-γ David W. Salt,† Brian D. Hudson,‡ Lee Banting,‡ Matthew J. Ellis,*,‡ and Martyn G. Ford‡ Department of Mathematics, Lion Terrace, Buckingham Building, University of Portsmouth, Portsmouth, PO1 3HE, and Centre for Molecular Design, Institute of Biomedical and Biomolecular Science, University of Portsmouth, King Henry I Street, Portsmouth PO1 2DY Received September 27, 2004

A novel molecular dynamics (MD) analysis algorithm, DASH, is introduced in this paper. DASH has been developed to utilize the sequential nature of MD simulation data. By adjusting a set of parameters, the sensitivity of DASH can be controlled, allowing molecular motions of varying magnitudes to be detected or ignored as desired, with no knowledge of the number of conformations required being prerequisite. MD simulations of three synthetic ligands of the orphan nuclear receptor PPARγ were generated in vacuo using Tripos’s SYBYL and used as the training set for DASH. Two X-ray crystal structures of PPARγ complexed with Rosiglitazone were compared to gain knowledge of the pharmacophoric conformation; this showed that the conformation of the ligand is significantly different between the two structures, indicating that there is no distinct conformation in which rosiglitazone binds to PPARγ but multiple binding modes. An investigation into simulation length was carried out. A simulation of 5 ns was found to give highly variable results, whereas a simulation of 25 ns gave a representative window of motion for molecules of this size. DASH was compared with Ward’s hierarchical cluster analysis method. The results show that DASH analysis is as good as Ward analysis in some areas (e.g. conformation identification) and is superior in others (e.g. speed and input size). Introduction Determination of solid-state configurations of drugs and their biomolecular complexes, by X-ray crystallography, is often seen as the ‘holy grail′ of drug design. However, the dynamic nature of the molecules when they are not in crystal form is not always considered fully, i.e., a NMR study of PPAR-γ has shown pronounced molecular flexibility prior to ligand binding, which is then lost on complex formation.1 There are also examples where the solid-state drug or drug-receptor complex may not reflect the bioactive conformations in solution.2-5 There appears to be a gradual move away from solid-state structure-based design6 toward solution-based screening, not only due to the rapidity of screening, but also the increased success of potency from leads discovered this way.7 Computational approaches to the elucidation of the conformations of drug molecules and their complexes, and their role in drug design, are now well proven.8 The structures derived from this type of modeling form the basis for the multitude of QSAR approaches used.9 Often, however, the methods used to obtain these structures lack both molecular mechanical and dynamical information.10,11 This paper takes an approach suitable for the rigorous analysis of very flexible bioactive molecules, applies novel data sampling techniques to extended molecular dynamics simulations of three ligands of PPAR-γ (Figure 1), and is intended to extend the existing pharmacophoric picture.12 * To whom correspondence should be addressed. E-mail: [email protected]. Phone: 023 9284 3623. Facsimile: 023 9284 3722 † Department of Mathematics, University of Portsmouth. ‡ Centre for Molecular Design, Institute of Biomedical and Biomolecular Science, University of Portsmouth.

The peroxisome proliferator-activated receptors (PPARs) are a subset of a ligand-activated nuclear transcription factor (NTF) superfamily that includes the estrogen (ER), retinoid (RXR) and farnesoid (FXR) receptors. The PPARs include a pharmacologically important member, PPAR-γ, which serves as a key regulator of adipocyte differentiation13 and glucose homeostasis,14 the modulation of which has a proven role in the treatment of type 2 diabetes mellitus.15 Cellular transcription ensues after (i) apoPPAR-γ binds its endogenous ligand, (ii) this ‘activated′ complex forms a heterodimer with a ligand ‘activated′ RXR receptor,16 (iii) displacement of co-repressor proteins from the PPAR and (iv) recruitment of partner transcriptional coactivators facilitate the successful targeting of responsive gene elements. Step (i) requires the PPAR-γ-ligand complex to undergo a conformational change, common to the members of this NTF superfamily, by which the ligand binding (LBD) and activation function 2 (AF-2) domains move. Much is known about the features of this process for the estrogen receptor (ER), and an increasing amount is being discovered about the molecular features of the ligands necessary for agonism (full and partial) and antagonism of estrogen-responsive element (ERE) transcription.17,18 However, less is known about the PPAR-γ-ligand complex interactions. The nature of the endogenous ligand of PPAR-γ is still proving elusive, although synthetic full/partial agonists and covalently modifying antagonists are known. The structural nature, in the solid state, of some of the PPAR-γ synthetic ligand complexes have been determined,19,20 and the dynamical nature of the solution phase interaction has been studied.21 The torsion angles provide a direct means of analyzing the conformations visited by a compound throughout a

10.1021/jm049216s CCC: $30.25 © 2005 American Chemical Society Published on Web 03/18/2005

DASH: A Novel MD Analysis Algorithm

Journal of Medicinal Chemistry, 2005, Vol. 48, No. 9 3215

MD simulation.10 Analytical programs that facilitate placing ensembles of torsion angles into separate groups, thus distinguishing between distinct conformations, are termed cluster analysis programs. A widely accepted method of cluster analysis for MD simulations is the Ward method, a hierarchical system that has wellknown limitations, not least of which is the inability to analyze large datasets.22 DASH is a novel algorithm for the analysis of MD simulation data and is a method that automatically identifies the number of distinct conformations (states) occupied during a MD simulation and classifies each sampled conformation into one of these states. DASH is introduced in this study as a superior method for conformational analysis of small druglike molecules. The DASH Algorithm (Dynamics Analysis by Salt and Hudson). Step 1. For each torsion angle separately, a categorical smoothing technique is applied23 that partitions the angles into a number (p) of disjoint and exhaustive regions. This is achieved by binning the values into class intervals of user-defined size and automatically searching for the positions of the localized modes. The cut points are then identified as being midway between these modal positions. The precision of this step in the algorithm can be improved by applying the categorical smoothing to a moving average representation of the angles that sharpens the positions of the modes. The order of the moving average is selected to achieve the required degree of separation between neighboring localized ranges. At the end of this step the cut points cutk (k ) 1, 2,......, p + 1) have been generated where -180 ) cut1 < cut2 < .... < cutp+1 ) 180. Because of the circular nature of the data, cut1 ) cutp+1. Step 2. Each disjoint range of angles is then assigned a range number Rk, such that all angles ai in this range are coded as Rk i.e., code ai ) Rk if cutk < ai < cutk+1. If the algorithm detects ‘wrap-around’ at ( 180, i.e., a localized range artificially split, then the right-most range is assigned the same code number as the first i.e., R1 ) Rp. Step 3. The number of distinct conformations for a particular molecule is calculated by combining the torsion angle codes from Step 2 in such a way that each distinct combination of torsion angle codes at time j produces a unique identifier that is termed the state code (state). For time point j the state code can be expressed as

statej ) cj1 + (cj2 - 1) × c1 + (cj3 - 1) × c1 × c2 +‚‚‚+ (cjNt - 1) × c1 ׂ‚‚cNt - 1 Nt

) cj1 +

(cji - 1) × Πki-1 ∑ ) 1ck i)2

(eqn.1)

where Nt is the number of torsion angles being considered and ci is the number of different codes for torsion angle i. Step 4. The final step involves pruning out the low occupancy, low probability states; it is intrinsic in the nature of MD simulations that as the length of a simulation increases, the chances of these states i) being visited and ii) then being recognized by DASH as a true state increases. The current procedure used to ensure

Figure 1. (a) The chemical structure of TZD01 (rosiglitazone); T1-T8 indicates the torsion angles that are present throughout all of the compounds in this study. The atoms marked with * were used for the distance measurements. (b and c) The chemical structure of TZD02 and TZD06, respectively.

Figure 2. Graph showing the effect the ‘smooth′ function of has DASH on the number of states that are identified for a 25 ns simulation of TZD01.

that all of these states are removed is to increase the ‘smooth function’ of DASH until a leveling off is observed. This can be seen in Figure 2; the plateau, which is observed between 50 and 90, shows the high occupancy states. For all of the simulations, a smooth value of ∼50 was found to be appropriate; this removes/ merges in all of the potential states which do not occur in bouts of more than 0.2% of a 25 ns simulation or 1% of 5 ns simulations. This leads to a consistent number of states through a series of simulations of the same compound (data not shown). DASH first splits each torsion angle into a parameterdefined group; thus all individual torsion angles that fall within a certain range will be given a single identifier (Rk) (Figure 3a-c). DASH then considers each

3216

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

Figure 3. (a-c) Torsion angle against sample conformation index for TZD01 with DASH torsion categorization shown in red; (d-f) frequency distributions for the representative conformational classes.

of these torsion angle identifiers with respect to what has occurred previously (Markov process), e.g. for DASH to identify a state change there must have been x changes in the torsion angle in the same direction of a predetermined magnitude, with the new state being occupied for a minimum of n counts. This enables DASH to disregard the high-speed movements of individual torsion angles while clearly identifying where significant movement has occurred (Figure 3a-c). The final process performed by DASH is to recombine these single torsion descriptors back into ensembles; these are then considered by DASH, again as a time series, using the ‘smooth′ function (Figure 2) to remove low occupancy states (a ‘smooth′ value of 50 will remove any states which occur in bouts of less than 50). Results and Discussion DASH Training Set. A small compound set of three PPAR-γ ligands (TZD02, TZD01 and TZD06 (Figure 1)) was chosen to represent a range of activities of synthetic PPAR-γ ligands (low, medium and high respectively). Details of the protocol used to generate MD simulations of these are given in Methods. Analysis of MD Trajectories. The analysis of torsion angle data was undertaken with the hope of identifying distinct conformations that a molecule may adopt as the MD simulation progresses. The distinct conformations are those that, except for relatively small

Salt et al.

changes in the torsion angles, occupy discrete regions of torsion angle space. A popular method utilized for identifying distinct conformations is cluster analysis performed on the observed torsion angles at each sampled time point, with their associated geometry being described by their centroids. There are three considerable drawbacks to this approach: (i) the method is very CPU intensive (several hours on a PC) particularly when the simulations are allowed to run for sufficient time to allow all important conformations to be represented; (ii) there is the inevitable problem of having to decide on how many clusters there are and finally; (iii) angular statistics such as torsion angles have the obvious feature that they are cyclical and thus give circular statistics. Methods of analysis such as classical cluster analysis, which is based on linear mathematics, are likely to produce misleading results if this property is not handled appropriately. An approach that (i) determines the number of distinct conformers visited during a MD simulation automatically, (ii) is computationally efficient and (iii) takes into account the circular nature of the data would be desirable. Sampled torsion angle values for TZD01 (Figure 1a), which have been plotted against the time at which the sample was taken, are shown in Figure 3a-c. It is clear from these time series that as MD simulations progress, the torsion angles at a given instant are residing in one of possibly several preferred localized ranges. Ignoring the time variable and viewing this information in frequency diagrams (Figure 3d-f), it can be seen that where these ranges are mutually exclusive they are characterized by unique modality. The situation is rather more complex when these localized ranges overlap with the merged frequency counts giving rise to larger localized ranges of torsion angle values that are either multimodal or exhibit tendencies toward being so. This is a significant problem for any conformation analysis method, which operates only on the frequency distributions. The combinations of localized ranges occupied determine the geometrical configurations or conformers. In Figure 3a-c we see three torsion angles plotted as time series: Torsion 3 (a) has only two localized ranges that are mutually exclusive, the modes being at approximately -90 and +90. Torsion 6 (b) has three localized ranges with modes at (again approximately) -180, -80 and +80 but has the added complication that the localized range centered on (180 is artificially split by the linear scale; a further complication can be seen with this torsion angle at approximately 4000 ps when there is observed a fourth localized range centered at about -50°. This short run of values can be seen to widen the range of values associated with those centered on -80°. The behavior of Torsion 7 (c) is more complex with several examples of these overlapping ranges. As can be seen from these plots, the torsion angles fluctuate about certain defined levels, moving from one level to another as the simulation progresses. In Figure 3a-c, the readout for DASH is shown by the black lines; the arbitrary state numbering employed by DASH has been replaced by the DASH centroid values. The reason for the superimposition is to illustrate that the correct number of codes have been identified, the change points separating the regions of

DASH: A Novel MD Analysis Algorithm

different codes occur at the correct times and the centroid values relate to the actual torsion data. As mentioned above, Torsion 3 (a) adopts two distinct preferences and therefore is given two categorical codes, i.e., two levels. In (b) and (c) we see that Torsions 6 and 7 respectively, which have increasingly more complex sets of change points, have both been appropriately coded. Since DASH has been designed with a range of ‘tuning′ parameters that allow for fine movements in torsion angles either to be recognized or ignored, it is easy to test out several parameters for the best result. Simulation Length. Originally a series of 5 ns MD simulations were performed on the 3 molecules (TZD01 TZD02 and TZD06; Figure 1a, b, and c respectively) to generate a test set for DASH. On repetition of these 5 ns MD simulations a markedly different set of results was gained where the numbers of conformations bore little resemblance between sister simulations. To investigate whether this was a problem with DASH or with the data itself, a distance monitor between atoms indicated in Figure 1a was used to perform a rough analysis of the simulations. This was a tertiary split based on the X-ray crystal structure of rosiglitazone (TZD01) bound to PPAR-γ (PDB entry: 1FM6),24 which exhibits a distance of 11.25 Å; the crystal structure has a resolution of 2.1 Å; therefore, distances of (1.05 Å were used (10.2 and 12.3 Å) giving the groups of conformations 12.3. Binary, tertiary and quaternary splits were all investigated; the tertiary splits showing the most promising results. The initial results from this distance analysis indicated that any two 5 ns simulations could exhibit a wide range of variation in the time spent in any of the distance groups. A major problem with MD simulations, especially where the time scales are short and/or the rotational barriers are not insignificant, is that the simulations can become trapped in low energy local minimum conformations for prolonged periods. This would make any MD analysis instantly flawed, as not all of the conformations that the compound is able to visit would necessarily be represented. Thus studies such as QSAR, which rely on accurate temporal structural information, would be impossible. This is suggested in the 5 ns simulation of TZD06 where the molecule exists in a range of folded conformations inconsistent with that of the crystal structure (PDB entry: 1FM9).24 A further simulation of TZD06, the most active compound in the series, was performed to investigate the effect of the length of simulation. This additional run was of 50 ns duration and showed that although TZD06 does spend ∼97.7% of its time in the folded conformations, the ‘crystal′ and extended conformations are in fact visited, if somewhat sporadically (Figure 4). This result is inconsistent with the hypothesis that the compounds have become trapped in local minima, as there is clearly free conformational movement even if folded conformations are dominant; since TZD06 is the most active of the molecules, the hypothesis that the pharmacophore is in the range of folded conformations can thus be proposed. The hypothesis that a 5 ns simulation is insufficient to sample all of the conformations visited by a molecule of this size is also supported. Longer time scale simulations (of 50 ns duration) were therefore performed for all three molecules, and the

Journal of Medicinal Chemistry, 2005, Vol. 48, No. 9 3217

Figure 4. 25 ns MD simulation of TZD06. The results are split into the three distance groups: (1) folded (>10.2 Å); (2) crystal (10.2-12.3 Å); (3) extended (10.2 Å) during 50 ns simulations. Each simulation is split into 10 × 5 ns samples, 5 × 10 ns samples and 2 × 25 ns samples. TZD01(a) is the original 50 ns simulation. TZD01(b) is two 25 ns simulations combined.

41.14% to 10.63% and 89.32% to 9.05% variance respectively; the higher variance seen for TZD01 and TZD02 compared with TZD06 is due to the degree of molecular flexibility inherent in the individual compound’s structure. Of course, all of these time lengths and occupation of conformations are critically dependent on the nature of the force field used, in particular the relative energies. The nature of the MD formalism ensures that the lowest energy conformations of a system will be more highly populated in any stochastic sampling scheme, such as the MD methods employed here, or other methods such as Monte Carlo. Other energy evaluation methodologies such as QM and DFT could be used and may give different populations dependent on the relative energies. The arguments above do not attempt to resolve this issue. Rather, within the parameters of the energy evaluation we have attempted to provide an internally consistent methodology that is independent of the particular energy function in use. It would be interesting to compare the MM energy function utilized here with similar dynamical simulations using higher order energy theories. Comparison with Ward’s Clustering. In order for an analysis method to prove the equal or superior of an already established and trusted method there are certain benchmarks that must be met. These are (i) the ability to generate results that are as relevant as, or better than, those generated by the old method, (ii) increased productivity, and (iii) a broader scope of employment. With the development of DASH, all of these criteria have been met when compared to the commonly used hierarchical cluster method Ward. The results generated by both DASH and Ward for TZD01 are shown in Figure 6a and b, respectively. The number of conformations into which the programs divide the data is dependent on a number of factors that differ between the two methods. With Ward analysis there is the choice of either directly setting the number of conformations required, or setting a similarity level. The

Figure 6. Plot of the states/conformations identified by (a) DASH and (b) Ward, respectively. All results have been numbered in order of their first appearance.

problem with the first of these approaches is that a certain level of accuracy must be attained in order for the results to be reliable, i.e. if the data is split into too many groups, conformations that are similar enough to belong in one group would be divided to form two or more smaller groups; conversely, if the data is split into too few groups, conformations which would usually be considered two or more classifications could be pulled together to create one larger group. A major advantage of DASH is that this arbitrary choice of the number of clusters expected is unnecessary once the relevant parameters have been set. However, computational speed is of no value if the results are inconsistent. It is therefore necessary to compare the quality of the results between the two methods. To achieve this, Ward clustering was performed using, as the number of clusters, the number of states identified by the DASH method. Thus for TZD01, where DASH predicts 23 discrete states, the number of Ward clusters was set to 23. This allows reasonable comparisons of the results of the two methods to be made. As DASH assigns an arbitrary state label to each of the states identified and Ward numbers them in order of first appearance, all of the DASH states were renumbered according to their order of appearance in the simulation to allow an easier comparison. Inspection of Figure 6a and b shows that similar processes are being carried out both by the DASH and the Ward algorithms. (i) At the time points between ∼1500 and ∼2500 ps a long occupancy is seen as state 6 (a) and conformation 13 (b); this shows that both programs are capable of not only identifying unique structures but accurately determining their transition points. (ii) There are similar progressions between states/conformations, for example in the last 1000ps of

DASH: A Novel MD Analysis Algorithm

the simulation; this steady increase in state number, mirrored in (a) and (b), indicates that both DASH and Ward recognize similar changes in structure. (iii) The recurrence of states that have been previously visited are also very similar between the two methods. This can most clearly be seen by looking at the repeated conformations which occur between ∼2500 and 3000 ps. Although there are clearly more states reported by Ward’s analysis, the repetition of structures first seen in the initial 500 ps of the simulation is clearly mirrored by both Ward and DASH. The differences between the two programs are equally apparent from Figure 6, primarily the splitting of single states reported by DASH into multiple conformations reported by Ward. The fundamental reason for this difference in splitting is that DASH considers the MD simulation data as a time series. Ward projects each of the ensembles of torsion angles into a matrix and then splits these according to the desired similarity level, disregarding any temporal information contained within the results. In light of the localized ranges of torsion angles observed in Figure 3a-c, and the similarities and differences observed between DASH and Ward’s analysis in Figure 6a and b, it can be proposed that DASH is not only as good at analyzing MD data as Ward, but improves upon it by taking into account, and then removing, the structurally insignificant, high-frequency movements. Although the receptor bound conformation of these ligands is unknown, the results from this study indicate that the active conformations of the PPARγ agonists are folded and have a distance measurement (Figure 1a) of under 10.3 Å. This does not agree with the ligand-bound conformation of rosiglitazone that has a distance measurement of 11.25 Å.23 There are two possible explanations for this: 1. The active conformation of these compounds may lie within a different distance group; to clarify this, more information is needed about the activities of the compounds. QSAR studies could then be performed which would identify the active conformation of the molecules. Another course of action would be to use docking search programs that simulate the binding of the ligand to the receptor. A rigid structure of PPARγ without rosiglitazone or any other bound ligand would be needed for this, as it has been shown that the structure of PPARγ is changed when ligand bound.1 2. The use of X-ray crystallography to visualize molecules presents its own problems. The picture that is obtained is an average structure, not a direct representation of the molecule. This means the crystal structure, although useful, is not an ideal representation. This view is supported by comparison of the conformations of TZD01 extracted from the X-ray structures found under PDB entry 1FM624 (Figure 7a and b) and CSD entry PONXUL (Figure 7c). The pronounced difference in conformation between the two structures suggests that rosiglitazone undergoes a conformational change on binding to PPARγ. It is important therefore not to rely on the single X-ray structure but to investigate the range of accessible conformations, which is the purpose of DASH. To investigate this further, a series of 24 extra compounds, all containing the backbone present in the three molecules used for this study (Figure 1a-c), are

Journal of Medicinal Chemistry, 2005, Vol. 48, No. 9 3219

Figure 7. The conformations of rosiglitazone extracted from the X-ray crystal structures 1FM6 (a and b, green) and PONXUL (c, yellow). The atoms used previously, which identify the ends of the torsion angle measurements, are colored red and blue. For all of the conformations, the thiozole group is oriented the same.

currently under investigation. 25 ns MD simulations have been generated for all 27 of these compounds and applied to DASH analysis. These dataset will be used for future QSAR studies. Conclusions MD offers an easy and effective way of characterizing the movements of molecules over time. The novel analysis program DASH, presented herein, has been shown to distinguish between the main conformational changes made by small flexible molecules in a fashion comparable to other analysis methods such as Ward hierarchical clustering. Its advantages are as follows: (i) it takes only seconds to run, whereas cluster analysis using the Ward method22 took 1.6 h; (ii) due to the treatment of the data as a time series, DASH can be used on larger data series than Ward; (iii) no information concerning the number of clusters is required prior to application of the algorithm. It can be concluded that the process of analyzing each of the torsion angles separately, classifying the main angles adopted and those surrounding them into single torsion classifications and combining these to give the molecular conformations is a superior approach for the analysis of MD data since only significant conformation changes that involve torsion angle movements that last for above a set duration are considered of structural importance. It is also suggested by this study that the X-ray crystal structures of PPARγ represent various bound conformations of rosiglitazone (TZD01),24 though the actual conformation that PPARγ is in when ligand binding occurs may be very different. This finding is supported by an NMR study carried out on ligand-bound PPARγ and apo-PPARγ,1 which shows decreased protein flexibility upon binding of rosiglitazone. Other studies have also indicated that the preconception of X-ray crystallography as providing a definitive representation of the binding mode of a ligand-substrate complex can

3220

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

be misleading. For example, the consideration of multiple binding modes in docking studies have been reported to yield a higher success rate than consideration of the crystal structure alone,25 and increasingly studies are reporting an induced fit or stabilization of proteins upon ligand binding.26,27 MD simulations, combined with DASH, provide a useful tool in the study of molecules the size of rosiglitazone, the results of which can be used to generate conformational information. The dataset of MD simulations for the TZDs has value in the development of analytical techniques and in further investigations such as QSAR studies. Methods Molecular Dynamics Simulations. Models of the compounds were constructed using standard fragments found within the 2D and 3D molecular databases of the Tripos Sybyl 6.9 package.28 Gasteiger-Huckel partial atomic charges were applied and the resultant structures were energy minimized using 100 steps of steepest descent minimization. Molecular dynamics (MD) simulations were performed in vacuo using (i) 6 ps, beginning at 50 K and increasing incrementally by 50 K every 1 ps to 300 K as a heating protocol, (ii) a 20 ps period at 300 K was employed for equilibration and (iii) 5 ns, 25 ns and 50 ns periods at 300 K were used in the simulations. Conformations were calculated every 1 fs of simulation and representative structures ‘captured′ every 1 ps. The standard SYBYL 6.9 force field was used for all simulations. Following the MD simulation, torsion angles T1-T8 and distance measurements (defined as for TZD01 in Figure 1a) were calculated for each series and exported for further manipulation. The torsion angle and time columns were converted to a space-delimited format upon which Ward’s cluster analysis was performed using a function in the MATLAB programming environment.29 DASH analysis was performed using an inhouse Visual Basic program.

Acknowledgment. We thank the EPSRC for their funding (M.J.E.) and Tripos Associates Ltd. for the contribution of their software, SYBYL 6.9. References (1) Johnson, B. A.; Wilson, E. M.; Li, Y.; Moller, D. E.; Smith, R. G.; Zhou, G. Ligand-induced stabilization of PPARγ monitored by NMR spectroscopy: Implications for nuclear receptor activation. J. Mol. Biol. 2000, 298, 187-194. (2) Gangjee, A.; Yu, J., McGuire, J. J.; Cody, V.; Galitsky, N.; Kisliuk, R. L.; Queener, S. F. Design, synthesis, and X-ray crystal structure of a potent dual inhibitor of thymidylate synthase and dihydrofolate reductase as an antitumor agent. J. Med. Chem. 2000, 43, 3837-51. (3) Li, D.; Levy, L. A.; Gabel, S. A.; Lebetkin, M. S.; DeRose, E. F.; Wall, M. J.; Howell, E. E.; London, R. E. Interligand Overhauser effects in type II dihydrofolate reductase. Biochemistry. 2001, 40, 4242-52. (4) Nugiel D. A.; Voss, M. E.; Brittelli, D. R.; Calabrese, J. C. An approach to the design of novel cognitive enhancers using molecular modeling and X-ray crystallography. Drug Des. Discovery 1995, 12, 289-95. (5) Sanders, W. J.; Nienaber, V. L.; Lerner, C. G.; McCall, J. O.; Merrick, S. M.; Swanson, S. J.; Harlan, J. E.; Stoll, V. S.; Stamper, G. F.; Betz, S. F.; Condroski, K. R.; Meadows, R. P.; Severin, J. M.; Walter, K. A.; Magdalinos, P.; Jakob, C. G.; Wagner, R.; Beutel B. A. Discovery of potent inhibitors of dihydroneopterin aldolase using CrystaLEAD high-throughput X-ray crystallographic screening and structure-directed lead optimization. J Med. Chem. 2004, 47, 1709-18. (6) Nienaber, V. Is X-ray crystallography a viable tool for lead discovery? Drug Discovery Today 2002, 7, 650. (7) van Dongen, M.; Weigelt, J.; Uppenberg, J.; Schultz, J.; Wikstro¨m, M. Structure-based screening and design in drug discovery. Drug Discovery Today 2002, 7, 471-478.

Salt et al. (8) Holtje, H. D.; Sippl, W.; Rognan, D.; Folkers, G. Molecular Modelling - Basic principles and applications, 2nd ed.; WileyVCH: Weinheim, 2003. (9) Devillers, J., Ed. Comparative QSAR; Taylor & Francis, Washington, D.C., 1998. (10) Ford, M. G.; Hoare, N. E.; Hudson, B. D.; Nevell, T. G.; Banting, L. QSAR studies of the pyrethroid insecticides. Part 3. A putative pharmacophore derived using methodology based on molecular dynamics and hierarchical cluster analysis. J. Mol. Graph. Model. 2002, 21, 29-36. (11) Hudson, B. D.; George, A. R.; Ford, M. G.; Livingstone D. J. Structure-activity relationships of pyrethroid insecticides. Part 2. The use of molecular dynamics for conformation searching and average parameter calculation. J. Comput.-Aided Mol. Des. 1992, 6, 191-201. (12) Kurogi, Y. Three-dimensional quantitative structure-activity relationships (3D-QSAR) of antidiabetic thiazolidinediones. Drug. Des. Discovery 1999,16, 109-18. (13) Grimaldi, P. A. The roles of PPARs in adipocyte differentiation. Progress in Lipid Res. 2001, 40, 269-281. (14) Willson, T. M.; Brown, P. J.; Sternbach, D. D.; Henke B. R. The PPARs: from orphan receptors to drug discovery. J. Med. Chem. 2000, 43, 527-50. (15) Parker, J. C. Troglitazone: the discovery and development of a novel therapy for the treatment of Type 2 diabetes mellitus. Adv. Drug Delivery Rev. 2002, 1173-1197. (16) Lee, M. S.; Kliewer, S. A.; Provencal, J.; Wright, P. E.; Evans, R. M. Structure of the retinoid X receptor alpha DNA binding domain: a helix required for homodimeric DNA binding. Science 1993, 260, 1117-21. (17) Mueller, S. O.; Hall, J. M.; Swope, D. L.; Pedersen, L. C.; Korach, K. S. Molecular determinants of the stereoselectivity of agonist activity of estrogen receptors (ER) alpha and beta. J. Biol. Chem. 2003, 278, 12255-62. (18) Rodriguez, A. L.; Tamrazi, A.; Collins, M. L.; Katzenellenbogen, J. A. Design, synthesis, and in vitro biological evaluation of small molecule inhibitors of estrogen receptor alpha coactivator binding. J. Med. Chem. 2004, 47, 600-11. (19) Sauerberg, P.; Pettersson, I.; Jeppesen, L.; Bury, P. S.; Mogensen, J. P.; Wassermann, K.; Brand, C. L.; Sturis, J.; Woldike, H. F.; Fleckner, J.; Andersen, A. S.; Mortensen, S. B.; Svensson, L. A.; Rasmussen, H. B.; Lehmann, S. V.; Polivka, Z.; Sindelar, K.; Panajotova, V.; Ynddal, L.; Wulff, E. M. Novel tricyclic-alphaalkyloxyphenylpropionic acids: dual PPARalpha/gamma agonists with hypolipidemic and antidiabetic activity. J. Med. Chem. 2002, 45, 789-804. (20) Cronet, P.; Petersen, J. F.; Folmer, R.; Blomberg, N.; Sjoblom, K.; Karlsson, U.; Lindstedt, E. L.; Bamberg, K. Structure of the PPAR-alpha and -gamma ligand binding domain in complex with AZ 242; ligand selectivity and agonist activation in the PPAR family. Structure (Camb) 2001, 9, 699-706. (21) Johnson, B. A.; Wilson, E. M.; Li, Y.; Moller, D. E.; Smith, R. G.; Zhou, G. Ligand-induced stabilization of PPAR-gamma monitored by NMR spectroscopy: implications for nuclear receptor activation. J. Mol. Biol. 2000, 298, 187-94. (22) Leach, A. Molecular Modelling; 2nd ed.; Prentice Hall: Upper Saddle River, NJ, 2001. (23) Hastie, T. J.; Tibshirani, R. J. Generalised Additive Models; Chapman and Hall: London, 1990. (24) Gampe, R. T., Jr.; Montana, V. G.; Lambert, M. H.; Miller, A. B.; Bledsoe, R. K.; Milburn, M. V.; Kliewer, S. A.; Willson, T. M.; Xu, H. E. Asymmetry in the Ppargamma/Rxralpha Crystal Structure Reveals the Molecular Basis of Heterodimerization Among Nuclear Receptors. Mol. Cell 2000, 5, 545. (25) Kallblad, P.; Mancera, R. L.; Todorov, N. P. Assessment of multiple binding modes in ligand-protein docking. J. Med. Chem. 2004, 47, 3334-3337. (26) Sotriffer, C. A.; Kramer, O.; Klebe, G. Probing flexibility and “induced-fit” phenomena in aldose reductase by comparative crystal structure analysis and molecular dynamics simulations. Proteins 2004, 56, 52-66. (27) de Olveira, C. A. F.; Guimaraes, C. R. W.; Barreiro, G.; de Alencastro, R. B. Investigation of the induced-fit mechanism and catalytic activity of the human cytomegalovirus protease homodimer via molecular dynamics simulations. Proteins 2003, 52, 483-491. (28) SYBYL 6.9; Tripos Inc.: St. Louis, MO, 2003. (29) MATLAB is available from The MathWorks Inc.: Natick, MA.

JM049216S