Vol. 47 No. 1/2000 1–9 QUARTERLY

A theoretical model for predicting the free energy of binding between anthracycline antibiotics and DNA was developed using the electron density functional (DFT) and molecular mechanics (MM) methods. Partial DFT-ESP charges were used in calculating the MM binding energies for complexes formed between anthracycline antibiotics and oligodeoxynucleotides. These energies were then compared with experimental binding free energies. The good correlation between the experimental and theoretical energies allowed us to propose a model for predicting the binding free energy for derivatives of anthracycline antibiotics and for quickly screening new anthracycline derivatives.

Cassinelli, 1998).In addition, the anthracyclines are arguably the best-characterized group of DNA intercalators.Numerous high-resolution crystallographic and NMR structures of daunorubicin and other anthracyclines bound to DNA have been reported (Chaires, 1996).The DNA sequence preference of daunorubicin binding to DNA has been characterized by thermodynamic, kinetic, and footprinting data (for reviews see Chaires, 1990;1995a;1995b).Yet, despite the intense efforts being made to synthesize more efficacious anthracyclines with fewer toxic side effects, doxorubicin remain the best chemotherapeutic agent of this class.In fact, it is becoming clear that the rational design of new anthracycline derivatives will require, in addition to structural information (see, e.g., Chaires et al., 1997;Lampidis et al., 1997; Priebe & Perez-Soler, 1993; Priebe et al., 1995), the development of new theoretical models.Some groups have already begun working on such models.One recent study devoted to the estimation of the free energy of binding of anthracyclines with DNA (Bagiñski et al., 1997) produced two interesting results: (1) the design of a thermodynamic cycle involving unbound and bound states of anthracyclines and (2) the determination of the free energy of binding as the sum of free energy changes along the path of the cycle.Although the theoretical and experimental results were in very good agreement, the reported approach is computationally intensive and very sensitive to the selection of some parameters of the model (see Bagiñski et al., 1997); in short, the total binding free energy is the sum of more than 10 individual contributions of both signs whose absolute value reaches 30-50 kcal/ mole, which in turn is almost 1 order of magnitude higher than the final result.More recently, a free energy component analysis has been used to model protein-DNA complexes (Jayaram et al., 1999).However, while this approach provides a good physical background for theoretical computations of the binding free energy, its practical application would be limited when comparing larger number of ligands.
In this light, the purpose of the present study was to formulate a method for a quick, practical analysis of anthracyclines derived from a common lead compound and of strongly related ligands having similar binding modes.In conducting our study, we followed Hermans et al. (1999) in assuming that the free energy of binding DG can be conceptually divided into the following four components: DG = DG i + DG s + DG c + DG a , where DG i -represents direct ligand-macromolecule interactions involving Van der Waals and Coulomb interaction energy of ligand atoms with atoms of the macromolecule; DG s -is the difference between the solvation energy of the macromolecule and the ligand in complex and in isolation; DG c -represents the confinement contribution in the bound state, the ligand being confined to a limited range of orientations and a small volume, DG a -represents conformational adjustments (i.e., adjustments in the conformations of the ligand and the macromolecule made to meet requirements of the bound state).
Although this decomposition is approximate, it is useful in aiding the understanding of the differences in binding among various compounds.The first and the fourth component can be determined using MM and/or MD methods.In the case of closely related compounds, the second and third components can be accounted for indirectly on the basis of the following qualitative analysis.The solvation component consists of two terms -the first one is proportional to the change of the total solvent accessible area, and the second one arises from the reorganization of the ionic density around the target molecule.Any change in the total solvent-accessible area is proportional to the area of contact between the target and the ligand, and the shape of the ionic cloud should not be to a large extent dis-turbed by the change of the ligand details.Since the non-bonded component will be roughly proportional to the surface of molecular contacts, it is reasonable to assume that the solvation component energy of binding will be correlated with the non-bonded energy.The confinement term can be neglected since it is very similar for closely related ligands having the same binding mode.
If all these assumptions are valid, then the free energy of binding of anthracyclines should be linearly correlated with the energy of binding computed by the MM method.With this in mind, we devised an MM-based model for calculating the free energy of binding between anthracycline antibiotics and DNA, and used it to calculate the free energies of binding for a series of previously reported as well as novel anthracycline derivatives.

MATERIALS AND METHODS
Using the MM method, we computed the binding energies of a DNA octamer with molecules from a set of anthracycline antibiotics for which the experimental free energies of binding were previously reported (Chaires et al., 1995).We also computed binding energies with this DNA octamer for a set of novel mono-and bis-intercalating compounds based on daunorubicin and doxorubicin (Figs. 1  and 2, respectively).Three-dimensional (3-D) models of all the anthracyclines used were built using an Insight II visual interface (MSI, 1997a) (Fig. 3) to manipulate the experimental structure of bis-daunorubicin (WP631) bound to 5¢-d(ACGTACGT)-3¢ (Robinson et al., 1997), PDB ID: 1AL9 (Bernstein et al., 1977;Sussman et al., 1998).The same octamer was used in the computations.In the case of monointercalating compounds, 3-D structures were generated, first by modeling the WP631 structure and then combining two symmetrically bound molecules.To ensure consistency of the computed energies, the sin-gle molecular structure was used as the reference starting structure.
The geometries of all structures were optimized using density functional theory (DFT) as implemented in the Dmol program (MSI, 1996).The atomic partial charges on the anthracycline antibiotics, reproducing the molecular electrostatic potential, so-called ESP-charges, were computed.All calculations were performed on a fine mesh, with the local VWN functional using the DNP basis set without frozen orbitals.
All MM calculations were carried out using the Amber force-field (Weiner et al., 1984;1986) as implemented in InsightII release 97.0.An implicit solvent model with the distance dependent "dielectric constant" e = r and reduced charges on the DNA phosphate groups was applied (see Rudnicki & Lesyng, 1995;1997).The Discover 2.97 program (MSI, 1997b) was used to calculate the energy minimizations.For an overview of the theoretical methods applied see Lesyng & McCammon (1993).In the case of monointercalating molecules, the binding energy per molecule was computed using the formula: ( ) The coefficient 1/2 appears in this equation because there are two antibiotics bound to each DNA molecule.In the case of bisintercalating molecules, the binding energy per molecule was computed using the formula: All MM energies were obtained for the optimized structures, with all energy gradients lower than 0.1 kcal/(mole Å).
Because DNA without an explicit solvent is too flexible, we therefore applied an auxiliary potential restraining DNA heavy atoms close to their initial, experimental positions.We then tested each molecule at four different force constants: 100 kcal/(mole Å 2 ), 10 kcal/ (mole Å 2 ), 1 kcal/(mole Å 2 ), and 0 kcal/(mole Å 2 ), the latter constant being one of the unconstrained optimization.The theoretical binding energies were then compared with the experimentally derived free energies of binding using a linear regression procedure.We compared the root mean square deviation (RMSD) between the initial and optimized structures, the correlation between the theoretical and experimental energy values, and the difference between the free energies of binding obtained experimentally and from our model.
The straightforward procedure described above required only one computationally intensive step -the computation of the atomic ESP charges of the relatively large drug molecules.But it was desirable to simplify the procedure even further by either computing the ESP charges more quickly (i.e., by minimizing the geometry using a fast force field minimizer and performing a single-point DFT energy calculation to find the ESP charges) or using a simpler procedure to obtain partial charges (i.e., by computing atomic charges using a heuristic algorithm).Opting for the second, simpler approach, we applied the charge equilibration method (Rappe & Goddard, 1991) implemented in Cerius2 release 3.8 (MSI, 1998) and then compared the results with the results obtained using the DFT-ESP charges.

RESULTS AND DISCUSSION
Energy minimization using different force constant values, k, of the potential restraining the DNA structure showed that the largest k leads to the best correlation between the experimental free energies of binding and the theoretical ones (Table 1).Conversely, All derivatives exhibit the same chiral properties as the original compounds.
smaller restraining potentials lead to a significantly higher RMSD between the minimized and experimental coordinates and a weaker correlation between the predicted and the experimental energy.It is worth noting, however, that building a model with good correlation between the computed binding energy and the experimental free energy of binding does not require very good description of the electrostatic interactions.Looking at calculations carried out using the partial charges, obtained by the heuristic procedure (Table 1, column E), one may notice that the correlation between the experimental and the computed energy obtained using this charge model and the strong restraining potential is comparable to the correlation obtained using the DFT-ESP charges and the weaker restraining potential.When one analyzes the force field binding energies and experimental free energies of binding for a selected set of previously reported anthracycline antibiotics (Table 2), the correlation between the theoretical and experimental values of the binding energies for k = 1 (1a) and k = 100 kcal/(mole Å 2 ) (1b) (Fig. 4), and the free energies of binding for several novel anthracycline derivatives calculated using our model (Table 3), it becomes clear that the binding energies computed using the force field method are much higher than the free en-  ergy of binding for the same components.However, the very clear correlation between these values, which certainly arises from a deeper relation between them, allows reasonably accurate predictions to be made about the free energy of binding.For example, one may notice that the difference between correlations obtained with the DFT-ESP partial charges and with the charges obtained by the equilibration method, is rather small.This shows that the details of the electrostatic interactions are relatively less important for the differentiating interactions between DNA and various anthracycline antibiotics than Van der Waals interactions and solvation.This significant dependence on the restraining potential suggests that a more rigid geometry is required to get a good correlation between the computed and experimental free energies of binding.It is therefore reasonable to assume that the disparity between the computed energy and the experimental free energy of binding for the anthracycline antibiotics can be attributed to interactions with the solvent.In fact, it has already been shown that some details of the solvation shell are strongly dependent on the structure of the solute (Rudnicki & Pettitt, 1997).This dependence can in turn justify the dependence of the correlation between the computed and experimental free energy of binding on the force constant of the restraining potential.

Table1. Correlations between the experimental free energies of binding (DG) and the theoretical energy of binding (DE) for various models
In conclusion, our proposed MM-based approach of linear correlation allows fast, effective, and reasonably accurate estimation of the free energy of binding of new anthracycline derivatives to DNA.This makes our model applicable to fast screening of derivatives of a common lead compound, given a good initial model of interactions between the lead compound and the target molecule, and allows it to be used to select potential drugs on Vol. 47 Binding of anthracyclines to DNA 7  the basis of binding affinity.If the binding of a ligand is comparable to or stronger than the one obtained for daunorubicin and doxorubicin (DE » -9 kcal/mol), then the molecule should be selected for further studies.There is no guarantee, however, that such a ligand will be a better drug, since there are many other factors involved and their influence can only be determin in further studies.

Figure 1 .
Figure 1.Schematic representation of monointercalating derivatives of daunorubicin a) and doxorubicin b).

Figure 2 .
Figure 2. Schematic representation of bisintercalating derivatives of anthracycline antibiotics.WP 631 is the basic compund described in (Robinson et al., 1997) and the derivatives represent the same chiral symmetry.Only a half of the symmetrical structure is presented.