Toward an Understanding of the Drug−DNA Recognition Mechanism

11480

J. Phys. Chem. 1996, 100, 11480-11487

Toward an Understanding of the Drug-DNA Recognition Mechanism. Hydrogen-Bond Strength in Netropsin-DNA Complexes Carlos Alema´ n,*,† M. Cristina Vega,‡ Lydia Tabernero,§ and Jordi Bella⊥ Departament d’Enginyerı´a Qı´mica, Escola Te` cnica Superior d’Enginyers Industrials de Barcelona, UniVersitat Polite` cnica de Catalunya, AVinguda Diagonal 647, E-08028 Barcelona, Spain, Departament de Biologia Molecular i Cellular, Centre d’InVestigacio´ i DesenVolupament, C.S.I.C., Jordi Girona 18, E-08034 Barcelona, Spain, Center for AdVanced Biotechnology and Medicine, 679 Hoes Lane, Piscataway, New Jersey 08854, and Department of Chemistry, Rutgers, The State UniVersity of New Jersey, Piscataway, New Jersey 08855 ReceiVed: January 10, 1996; In Final Form: April 12, 1996X

An ab initio quantum mechanical study of the binding of basic minor-groove drugs to DNA has been undertaken considering three interacting model systems, two concerning hydrogen bonds between amide moieties and one involving an interaction with a charged group. These are thymine‚‚‚formamide, thymine‚‚‚Nmethylacetamide, and thymine‚‚‚thyleneiminium. The effect of the solvent on the interaction energy has been explored by using a self-consistent reaction field (SCRF) method based on the high-level MiertusScrocco-Tomasi algorithm. Furthermore, we have made a comparison of the intermolecular geometries of drug-DNA interactions by extracting information from the Nucleic Acid Data Base. The results indicate that interactions involving a charged group are about 5 times stronger than hydrogen bonds between noncharged groups in a gas-phase environment. However, both types of interactions are greatly modulated by the solvent. Thus, whereas a hydrogen bond between noncharged groups is clearly a hydrophobic interaction, the strong polarization effect induced by the charged group would eliminate the unfavorable effect of the solvent if a small variation of the intermolecular geometry is considered. These results suggest that interactions involving charged groups play a crucial role in the drug-DNA recognition and binding mechanism.

Introduction Detailed knowledge of the interactions of DNA-binding ligands with DNA should allow the rational design of ligands with improved pharamacological properties or as usefulness research tools. Detailed three-dimensional structures of drugDNA complexes have been reported during the last years.1-16 The results have provided important information about both sitespecific and nonspecific nucleic acid recognition. However, a more detailed analysis of the different binding sites found in ligand-DNA structures is required for the understanding of many biological properties. Netropsin (Figure 1) is a DNA minor-groove-binding drug widely studied by multiple techniques.17-21 It is a long, flexible molecule with positively charged ends and a number of proton donor and acceptor groups for possible hydrogen bonds spaced between. There are several single-crystal X-ray diffraction studies of this drug complexed with different DNA sequences.3,4,7,12 The different structures indicate four important features: (i) the binding preference for stretches of A‚T base pairs over G‚C ones, (ii) the crucial role of hydrogen bonding interactions with O2 of thymines and N3 of adenines, (iii) the strong interactions of the positive charged ends with double helix, and (iv) the importance of the van der Waals interactions between the drug and the sugar rings that form the walls of the minor groove. Indeed, hydrogen bonds involving both charged and noncharged groups are arguably the most important interactions in drug-DNA complexes. * All the correspondence to this author. † Universitat Polite ` cnica de Catalunya. ‡ Centre d’Investigacio ´ i Desenvolupament, C.S.I.C. § Center for Advanced Biotechnology and Medicine. ⊥ Rutgers, The State University of New Jersey. X Abstract published in AdVance ACS Abstracts, June 1, 1996.

S0022-3654(96)00126-8 CCC: $12.00

Figure 1. Structure of netropsin.

Several models22-25 have been developed on the basis of experimental studies to elucidate the exact roles played by charged ends of netropsin and by hydrogen bonds. However, they still remain unclear at present time. The charged ends appear to be important since they are a commom feature of many minor-groove-binding drugs, e.g., Hoechst 33258 and distamycin. However, a netropsin derivative with both ends removed still binds to poly(dA)‚poly(dT), although of course less efficiently than does netropsin.26 On the other hand, the importance of the hydrogen bond between noncharged groups on drug-DNA recognition mechanisms has been thought to be essential. However, there are some experimental data that do not give support to this hypothesis and make it remain open. Thus, bisquaternary ammonium heterocycle, SN 18071, which © 1996 American Chemical Society

Drug-DNA Recognition Mechanism

J. Phys. Chem., Vol. 100, No. 27, 1996 11481

a

b

c

Figure 2. Structure of the energy minima found for the complexes thymine‚‚‚formamide (a), thymine‚‚‚dimethylacetamide (b), and thymine‚‚‚ methylenimidium cation (c), at the HF/6-31G*, HF/3-21G, and HF/6-31G* levels, respectively.

cannot form hydrogen bonds binds also to DNA, showing a similar A‚T minor groove specificity.27,28 In this work we present an energetic analysis of the drugDNA interactions. We investigate the relative energies of hydrogen bonds involving noncharged and charged groups both in vacuo and in solution using quantum chemical calculations. Three complexes were investigated: thymine‚‚‚formamide, thymine‚‚‚N-methylacetamide, and thymine‚‚‚methyleneiminium cation. In order to compare and explain these results we perform an analysis of the interactions found in X-ray determined structures. Methods Selection of the Reduced Models in the Representation of the Drug-DNA Interactions. The size of the netropsin molecule prevents relevant ab initio analyses from being performed. We have therefore prepared smaller model complexes with equivalent interactions. Three complexes have been investigated in the present work: thymine‚‚‚formamide (A), thymine‚‚‚N-methylacetamide (B), and thymine‚‚‚methyleneiminium cation (C). Hydrogen bonds formed between the amide groups of the drug and DNA are mimicked by A and B, whereas C mimics the hydrogen bonds formed by guanidinium and propionoamidinium groups.

The selection of the model molecules has been performed according to their usefulness in representing the charge distributions of the netropsin groups involved in the interactions of study. Accordingly, electrostatic (potential-derived) charges were determined following Momany’s procedure,29 where the quantum mechanical molecular electrostatic potentials (MEPs) are fitted to a set of monopoles centered at the nuclei. The MEPs were computed from ab initio 6-31G*30 wavefunctions, on points located at four Connolly layers31 placed at 1.4, 1.6, 1.8, and 2.0 times the van der Waals radii. The charges were computed for (a) formamide, (b) the central fragment of netropsin, which contains an amide group attached to two coplanar pyrrole groups (fragment 1), (c) the methyleneiminium cation, (d) the fragment of netropsin which contains the guanidinium group (fragment 2), and (e) the fragment of netropsin which contains the propionoamidinium group (fragment 3). Results are showed in Table 1. The atomic point charges on atoms which are commom to the netropsin fragments and small model molecules point out that both formamide and methyleneiminium cation effectively mimics the behaviour of the drug with respect the parameters of interest. On the other hand, N-methylacetamide was selected as a model molecule in order to mimic the progression of the drug molecule. Accord-

11482 J. Phys. Chem., Vol. 100, No. 27, 1996

Alema´n et al. is the difference of temperature dependent terms which are due to the Bolzmann distribution of the vibrational modes. The distortion energy (EDIS) was computed according to eq 3:

EDIS ) ECOMP - EOPT

Figure 3. Hydration free energy profiles for the complexes thymine‚‚‚formamide (solid line) and and thymine‚‚‚methylenimidium cation (broken line).

TABLE 1: Atomic Point and Global (∑Qi) Charges on Amide and Iminium Groups for Small Model Molecules (Formamide and Methyleneiminium Cation) and Selected Fragments of Netropsina

qN qC qH(N) qH(N) qO(C) ∑Qi a

fragment 1

methyleneiminium

fragment 2

fragment 3

-0.931 0.664 0.426

-0.751 0.655 0.383

-0.580 -0.421

-0.612 -0.325

-1.299 1.143 0.626 0.626

-1.441 1.539 0.575 0.536

-1.115 1.183 0.499 0.524

1.096

1.209

1.091

For more details see text.

ingly, complex B provides a more complete model from a structural point of view. Energy Calculations. The ab initio LCAO-MO method was used for the study of all the complexes. Geometry optimizations were performed at the Hartree-Fock (HF) 3-21G32 and 6-31G* 30 (d-type polarization functions on the heavy atoms) levels. The optimized structures were checked by analysis of harmonic vibrational frequencies obtained from diagonalization of force constant matrices to verify the minimum nature of the stationary points. Single point energy calculations were performed on the HF/6-31G* geometries using the 6-31++G** basis set (two sets of diffuse functions, d-type polarization functions on the heavy atoms and p-type polarization functions on the hydrogen atoms). The counterpoise (CP) correction method33 was applied to the HF/3-21G and HF/6-31G* geometries in order to obtain estimates of the basis set superposition error (BSSE) at these levels. The interaction energy was calculated as the difference of total energy of the complex (EA‚‚‚B) and the sum of the subsystem energies (EA, EB).

EINT ) EA‚‚‚B - (EA + EB)

(1)

All four compounds, i.e., thymine, formamide, N-methylacetamide, and methyleneiminium cation, were also fully optimized at the HF/3-21G and HF/6-31G* levels. The difference in vibrational energy of the complex and the dissociated compound has been taken into account by using vibrational data obtained from force constant analyses. Thus, the thermodynamics in vacuo was computed by correcting the energies of interaction to enthalpies of association at 298 K by using eq 2:

∆H298 ) EINT + ∆ZPE + ∆TC - RT

where ECOMP is the energy of the molecules in the complex and EOPT is the energy obtained from the optimization of the isolated molecule. The net binding energy is the total enthalpy of association computed from eq 2 minus the combined thymine and model molecule distortion energies. Solvation Calculations. The effect of water on the interaction energies was examined following the self-consistent reaction field (SCRF) developed by Miertus, Scrocco, and Tomasi (MST),34,35 which gives good estimations of the free energies of hydratation.36,37 This strategy computes the total free energy as the addition of three contributions: cavitation, van der Waals, and electrostatic:

∆Gsolv ) ∆Gcav + ∆Gvw + ∆Gele

formamide

(2)

where ∆ZPE is the difference of zero-point energies and ∆TC

(3)

(4)

The cavitation was computed by using the Pierotti’s theory,38 and the van der Waals contribution was evaluated from the molecular surface area by using a linear relation found by Tomasi and co-wokers.39 The electrostatic contribution was computed by using the SCRF method at both ab initio 6-31G* (MST/6-31G*) and semiempirical AM1 (MST/AM1) levels. MST/AM1 calculations were performed by using the algorithm developed by Luque and Orozco.36,40 The SCRF method represents the effect of the solvent by a set of virtual charges {Qi} located on the molecular cavity surface, which reacts against the solute charge distribution generating a reaction field. This effect is accounted for by adding a new term into the gasphase solute Hamiltonian (Ho), which lead to the definition of a new Hamiltonian:

H ) Ho + ∑iQi/Ri

(5)

The solute/solvent interface was built by using a moleculeshaped algorithm.41 According to previous studies36,42,43 the cavity was located at 1.25 (MST/6-31G*)/1.20 (MST/AM1) times the van der Waals radii of the different atoms [C,N ) 1.5 Å, O ) 1.4 Å, H ) 1.2 Å, H(bound to heteroatoms) ) 1.0 Å]. All the macroscopic parameters defining the water in the continuum calculations were taken at 298 K. Analysis of Structural Data. Nineteen structures of oligonucleotide-minor-groove binder complexes were chosen from the Nucleic Acid Database.44 The crystallographic R factors and the maximum atomic resolutions are in the range 14-21% and 1.90-2.50 Å, respectively.1-5,7-13 The analysis was performed first, by selecting those thymines or cytosines that had any interaction with donor groups of the ligand. Thereafter, subsets of atomic coordinates were generated from each occurrence that included the nucleotide base and the nitrogen donor atoms with its immediate neighboring bound atoms within the drug molecule. The subsets were transformed to the coordinate system of the calculated model, taken as a fixed reference. All the structures were analyzed by using the graphics program CHAIN45 and subsequently classified into two major groups: type I, in which the donor nitrogen was part of an amide moiety and type II, in which the donor nitrogen was part of an amidinium or guanidinium moiety. Initially, a restrictive criteria for the hydrogen bond distances of the selected interactions was applied. However, for the benefit of a meaningful statistics

Drug-DNA Recognition Mechanism

J. Phys. Chem., Vol. 100, No. 27, 1996 11483

TABLE 2: Intermolecular Geometriesa of the Three Complexesb Investigated by Using ab Initio Quantum Mechanical Calculations and of the Crystal Interactionsc,d Used as Starting Points in the Geometry Optimizations (Tabernero et al., 1993) Ab Bb Cb N1‚‚‚O2(T9)c N6‚‚‚O2(T19)d

dO‚‚‚H

dO‚‚‚N

∠O‚‚‚H-N

2.14 2.04 1.72 1.63 1.78

3.04 2.95 2.68 2.52 2.70

149.0 148.7 154.5 145.4 154.6

a Bond lengths in Å; bond angles in deg. b A, HF/6-31G* optimized geometry of the thymine‚‚‚formamide complex; B, HF/3-21G optimized geometry of the thymine‚‚‚N-methylacetamide; C, HF/6-31G* optimized geometry of the thymine‚‚‚methyleneiminium cation. c Electrostatic interaction in the d(CGCAAATTTGCG)‚‚‚netropsin crystal structure (Tabernero et al., 1993). d Hydrogen bonding interaction in the d(CGCAAATTTGCG)‚‚‚netropsin crystal structure (Tabernero et al., 1993).

analysis, a cutoff criteria of < 3.5 Å was generously adopted, and only those potential hydrogen bonds that lack an appropiate geometry were discarded. It is worth commenting that most of these structures have been solved at moderate atomic resolution, and therefore the accuracy for the nonbonded interatomic distances is somehow poor. Ab initio calculations were performed with Gaussian-9246 and Monstergauss47 computer programs. Semiempirical calculations were carried out with a modified version of the MOPAC48,49 computer package. MEPs and electrostatic charges were determined with the MOPETE/MOPFIT set of programs.50 Hydrogen atom positions for each donor group were generated in the analysis of structural data by using the standard parameters in X-PLOR.51 All the calculations were performed on the CRAY-YMP and IBM/3090 at the Centre de Supercomputacio´ de CAtalunya (CESCA). Results Geometries and Energies. Starting structures for the geometry optimizations of the three complexes were taken from the complex of the dodecamer d(CGCAAATTTGCG) with netropsin.12 The sugar attached to thymine was replaced by a hydrogen atom, whereas the drug was replaced by formamide (A) and N-methylacetamide (B) in the N6(netropsin)‚‚‚O2(thymine 19) position and by methyleneiminium cation (C) in the N1(netropsin)‚‚‚O2(thymine 9) position. Figure 2 shows the HF/6-31G* optimized structures of the A and C, as well as the HF/3-21G structure of B. The three complexes were characterized as minima by analysis of harmonic vibrational frequencies. Direct comparison of the theoretical geometries and the X-ray determined structures should be carried out with some caution, since the latter are for thymine attached to a furanose ring and the amide and charged groups are forming part of the drug. Furthermore, experimental geometries have been determined from molecules embedded in crystals, and under such circumstances molecular geometries are subject to external packing forces. Considering these uncertainties, good overall agreement between both HF/3-21G and HF/6-31G* ab initio data and crystallographical bond lengths and angles were found for the three complexes. Table 2 shows the geometric parameters defining the molecular interactions of: i) the three complexes and (ii) N1‚‚‚O2(T9) and N6‚‚‚O2(T19) in the d(CGCAAATTTGCG)‚‚‚netropsin complex.12 The agreement between the theoretical and X-ray values is quite good. Indeed, the small differences found may be attributed to structure refinement and limited

resolution of X-ray data. Indeed, hydrogen bonded geometries recently reported in a study of lexitropins analogues suggest that crystal packing may be distortive for certain complexes.52 It is interesting to note that the intermolecular distances for complex C are much smaller than for complexes A and B, indicating that interaction involving a charged group is much stronger than interactions involving amide groups. Comparison between A and B reveals that both formamide and N-methylacetamide give a good representation of the amide group of the drug, indicating that the end-methyl groups of the latter do not play a crucial role in the interaction. The close resemblance between A and B suggests that the HF/6-31G* geometry optimization of the latter is not necessary to improve the description of the hydrogen bond given by A. In Table 3 we report the interaction energies and enthalpies of association for the three complexes investigated. We denote the calculations by specifying the method used for geometry optimization after a double slash and the method used at this single point for determination of the energy before the double slash; e.g., HF/6-31G*//HF/3-21G denotes a single point calculation with the 6-31G* basis set at a geometry obtained from a 3-21G basis set. As can be noted the 3-21G basis greatly exaggerates the interaction energies, especially in A and B. The values calculated with the extended basis sets 6-31G* and 6-31++G** are more reasonable. The application of the CP correction with the 6-31G* basis permits the recovery of 1.43, 2.12, and 1.66 kcal/mol for complexes A, B, and C, respectively. However, with 3-21G the CP method is able to recover 6.04, 7.51, and 4.02 kcal/mol for A, B, and C, respectively. These values compare reasonably with those of Hrouda et al.53 for the Watson-Crick AT base pair, who found at the HF/MINI-1// HF/MINI-1 and HF/MIDI-1//HF/MINI-1 levels BSSEs of 3.9 and 6.4 kcal/mol, respectively. More recently, Gould and Kollman54 using the 6-31G* basis set found superposition errors of 2.42 kcal/mol for the Watson-Crick and Hoogsteen AT’s and 2.86 kcal/mol for the Watson-Crick GC. In general, the effect of enlarging the basis set may be estimated to high accuracy (say small BSSE) at the HF level.55 Thus, the addition of diffuse functions will not lead to an important increase of the BSSE.56 For the 6-31++G** interaction energies an estimate of the BSSE was taken as that calculated with the 6-31G* basis set. The enthalpies of association at 298 K for the three complexes are reported in Table 3. The best estimates for A and C are denoted as HF/6-31++G**//HF/6-31G* and lead to values of -3.78 and -23.55 kcal/mol, respectively. For B the best estimate of ∆H298 is denoted HF/6-31++G**//HF/3-21G and lead to a value of -2.36 kcal/mol. The energy difference found between the two hydrogen bond complexes must be attributed basically to the different computational levels used in the geometry optimizations. Thus, A and B give similar ∆H298 (5.65 and 5.14 kcal/mol, respectively) when compared to geometries optimized at the same level (HF/3-21G//HF/3-21G). The inclusion of d-type polarization functions on geometry optimizations results in a better estimate of the association enthalpies.54 The hydrogen bond strengths of complexes A and B compare well with those reported for the cytosine‚‚‚formamide complex at the HF level using the standard 6-31G basis set augmented by d polarization functions on the amino group nitrogens.57 The authors obtained interaction energies ranging from -2.19 to -4.18 kcal/mol depending on the planarity of the cytosine amino group nitrogens. However, it must be stressed that in such cases both cytosine and formamide molecules were kept in the same plane, mimicking a base pair

11484 J. Phys. Chem., Vol. 100, No. 27, 1996

Alema´n et al.

TABLE 3: Interaction Energies and Enthalpies of Association (in kcal/mol) of the Three Complexes Investigated ∆ZPEb

∆TC - RTb

∆H298b

1.56 1.70

-0.57 -0.56

-5.65 -4.96 -3.78d

-0.53

-5.14 -3.94e -2.36e

-1.37 -1.39

-31.24 -24.32 -23.55d

energy//geometrya

EINTb

HF/3-21G//HF/3-21G HF/6-31G*//HF/6-31G* HF/6-31++G**//HF/6-31G*

-12.68 -7.53 -6.34

HF/3-21G//HF/3-21G HF/6-31G*//HF/3-21G HF/6-31++G**//HF/3-21G

-13.46 -6.87 -5.29

Thymine‚‚‚N-Methylacetamide (B) -5.95 1.34 -4.75 -3.17

HF/3-21G//HF/3-21G HF/6-31G*//HF/6-31G* HF/6-31++G**//HF/6-31G*

-34.70 -25.58 -24.81

Thymine‚‚‚Methyleneiminium (C) -30.68 0.81 -23.92 0.99 -23.15c

EINT,CPb Thymine‚‚‚Formamide (A) -6.64 -6.10 -4.92c

a Level of energy calculation//level of geometry optimization. b E INT denotes the interaction energy, EINT,CP the interaction energy including a counterpoise correction, ZPE the zero point energy, TC the temperature dependent part of the vibrational partition function, and ∆H298 the enthalpy of association at 298 K. c BSSE calculated with the 6-31G* basis set. d Includes correction for zero point energy and temperature-dependent terms determined at the HF/6-31G*//HF/6-31G* level. e Includes correction for zero point energy and temperature-dependent terms determined at the HF/3-21G//HF/3-21G level.

TABLE 4: Energies of Distortion and Net Enthalpies of Association at 298 K (in kcal/mol) for the Three Complexes Investigated HF/3-21G// HF/3-21Ga

HF/6-31G*// HF/6-31G*a

Thymine‚‚‚Formamide (A) EDIS(Thymine) 0.28 EDIS(Formamide) 0.62 ∆H298b -5.65 ∆H298,netb -4.75 ∆H298,net(6-31++G**)c

0.14 0.28 -4.96 -4.54 -3.36

Thymine‚‚‚N-Methylacetamide (B) EDIS(thymine) 0.33 EDIS(N-methylacetamide) 1.26 b ∆H298 -5.14 ∆H298,netb -3.55 ∆H298,net(6-31++G**)c

-3.94d -2.35e -0.84

Thymine‚‚‚Methyleneiminium (C) EDIS(thymine) 1.71 EDIS(methyleneiminium) 1.33 ∆H298b -31.24 ∆H298,netb -28.20 ∆H298,net(6-31++G**)c

1.44 0.51 -24.32 -22.37 -21.60

a Level of energy calculation//level of energy optimization. b ∆H 298 denotes the enthalpy of association at 298 K (see Table 2), ∆H298,net c the net enthalpy of association. ∆H298,net(6-31++G**) denotes the best estimates of the net enthalpy of association at 298 K (see text). d Computed at the HF/6-31G*//H3-21G level (see Table 2). e Computed by using the distortion energies at the HF/3-21G//HF/3-21G level.

interaction, whereas in the present work thymine and formamide are twisted, mimicking a base-drug interaction. In order to obtain a more reliable estimation of the bonding interaction of the various complexes, the net enthalpies of association were used. They were calculated by substrating the thymine distortion and model molecule distortion energies. These distortion energies reflect induced fits that permit stronger intermolecular interactions. Table 4 shows thymine and model molecule distortion energies obtained by substracting the energies of optimized isolated molecules from the energies that occur in the complex. As can be noted complex C presents the larger distortion energies, which give a measure of the strength of this interaction. Thus, energies of 1.71 and 1.44 kcal/mol were obtained for thymine at the HF/3-21G and HF/ 6-31G* levels respectively, in comparison with 0.28 and 0.14 kcal/mol for complex A. The same behavior was observed for the model molecules. Such large energies are attributed to structural changes caused by the interactions that link thymine with model molecules. Comparison between A and B points

out the large flexibility of the end-methyl groups, since the distortion energy for N-methylacetamide is 0.64 kcal/mol greater than that of formamide. The best estimates of ∆H298 minus the combined thymine and model molecule distortion energies at the HF/6-31G*//HF/6-31G* level (HF/3-21G//HF/3-21G for complex B) give the best estimates of the net enthalpies of association. These are denoted by ∆H298, net (6-31++G*) in Table 4 and lead to -3.36, -0.87, and -21.60 kcal/mol for A, B, and C, respectively. Energy Calculations in Water. A detailed knowledge of the free energy associated with hydrogen bonds involving charged and noncharged groups in aqueous environment is essential to our understanding of the energetics and recognition mechanisms in drug-DNA association. In this section we compare our in vacuo results with electrostatic interaction and hydrogen bond formation in water. The free energies of solvation were computed for complexes A, B, and C by using the MST/6-31G* and MST/AM1 methods. Since the change of molecular geometry upon solvation has a negligible effect on the thermodynamic parameters,58,59 geometries obtained from gas-phase optimizations at the HF/6-31G* level were kept frozen during SCRF calculations. The total free energy was divided into steric and electrostatic contributions in order to facilitate the analysis. Within the first term we included the expansion term, the cavitation, and the van der Waals contributions to the free energy. The second term accounts for electronic effects. The ∆G of solvation for the complexes and their fragments are displayed in Table 5. The ∆G of solvations recalls the great solubility of the three complexes and fragments, which is particularly remarkable for complex C and methyleneiminium cation. Thus, the introduction of a positive charge in the system yields large polarization effects, with the consequent increase of the electronic contribution to the ∆G of solvation. The ∆∆Gs of solvation for complexes were determined from ∆G of solvation for fragments and are also included in Table 4. The results indicate that both types of interactions are greatly modulated by the solvent. Thus, fragments are better hydrated than complexes, suggesting a poor tendecy to form molecular associations. It is worth noting that complexes A and B give very similar ∆∆G values. This fact together with the gas-phase results (see Table 3) indicate that A is a good model for the description of the hydrogen bonding interaction between noncharged groups. It is encouraging to note the remarkable agreement between the MST/6-31G* and MST/AM1 results. This points to the suitability of the MST/AM1 to capture the essential features of

Drug-DNA Recognition Mechanism

J. Phys. Chem., Vol. 100, No. 27, 1996 11485

TABLE 5: Steric [van der Waals + cavitation + (∆PV)] and Electrostatic Contributions to the Total Free Energy of Hydration (in kcal/mol) of Complexes A and Ca Obtained from MST/6-31G* and MST/AM1 (in Parentheses) Calculations A

B

C

complexb

∆G(ster)

∆G(elec)

∆G(tot)

complex thymine formamide ∆∆Gc complex thymine N-methylacetamide DDGc complex thymine methyleneiminium cation ∆∆Gc

4.46 (4.09) 3.40 (3.21) 2.77 (2.31) -1.71 (-1.43) 4.58 (4.26) 3.40 (3.21) 3.10 (2.66) -1.92 (-1.61) 3.39 (3.19) 3.40 (3.22) 1.60 (1.48) -1.61 (-1.51)

-21.89 (-18.82) -17.28 (-13.64) -11.65 (-11.16) 7.04 (5.98) -20.21 (-17.54) -17.56 (-13.84) -13.01 (-9.82) 10.36 (6.12) -67.87 (-66.59) -17.83 (-14.03) -74.37 (-74.07) 24.33 (21.51)

-17.43 (-14.73) -13.88 (-10.43) -8.88 (-8.85) 5.33 (4.55) -15.63 (-13.28) -14.16 (-10.63) -9.91 (-7.16) 8.88 (4.51) -64.48 (-63.40) -14.43 (-10.81) -72.77 (-72.59) 22.72 (20.00)

a Geometries optimized in gas-phase at the HF/6-31G* level. b A, thymine‚‚‚formamide complex; B, thymine‚‚‚N-methylacetamide complex; C, thymine‚‚‚methyleneiminium cation complex. c ∆∆G(complex) ) ∆G(complex) -[∆G(thymine) + ∆G(ligand)].

TABLE 6: Statistical Results of the Comparison between the Hydrogen Bonding Geometry of the Calculated Models with That Determined Experimentally in Crystal Structures d(H‚‚‚O), Å

∠CdO‚‚‚H, deg

∠CdO‚‚‚N, deg

∠O‚‚‚H-N, deg

Type I 115.9 157.6 41.7 139.3 14.1

124.5 171.2 46.7 152.1 15.1

88.6 141.4 52.9 118.6 15.8

149.0

114.7

105.2

91.58 171.9 80.3 139.4 25.8

83.56 145.0 61.4 117.2 23.8

149.8

158.7

d(N‚‚‚O), Å

a

minimum maximum range mean std devb

1.79 3.08 1.29 2.42 0.45

2.58 3.37 0.79 2.99 0.32

model Ac

2.14

3.04

minimum maximum range mean std devb

1.64 3.38 1.73 2.49 0.54

2.49 3.48 0.99 3.03 0.33

model Cc

1.72

2.68

Type IIa 83.76 170.2 86.4 141.3 24.0 154.5

a

Type I group refers to hydrogen bonds in which the donor nitrogen was part of an amide moiety; Type II group refers to hydrogen bonds in which the donor nitrogen was part of an amidinium or guanidinium moiety. b std dev refers to standard deviation. There are a total of 10 and 14 interactions for Type I and Type II respectively. c Model A refers to the HF/6-31G* optimized geometry of the thymine‚‚‚formamide complex; and model B refers to the HF/6-31G* optimized geometry of the thymine‚‚‚methyleneiminium cation.

the solvation process at a reasonably computational cost. Figure 3 shows the hydration energy profiles for the approach of the two fragments in complexes A and C relative to the isolated fragments computed at the MST/AM1 level. The reaction coordinate (d) used to describe the approach of the two fragments was chosen as the intermolecular distance O(thymine)‚‚‚H(model molecule). All the remaining parameters were kept frozen at the ab initio HF/6-31G* minimum energy structures. The results clearly shows that water destablizes the formation of the complexes. Thus, the approach of the two molecular fragments leads to a systematic increase in the hydrophobicity of the system, which produces a positive ∆∆Gsolv. Two important features can be observed from comparison of hydration energy profiles for complexes A and C. First, ∆∆Gsolv of complex C is around 15 kcal/mol larger than ∆∆Gsolv of complex A. This must be attributed to the polarization effects induced by the positive charge of methyleneiminium cation. Second, ∆∆Gsolv of complex C displays stronger changes with the reaction coordinate than that of complex A. For instance, a small increment of 0.2 Å from the gas-phase minimum energy distance produces a reduction of 1.80 and 0.23 kcal/mol in ∆∆Gsolv of complexes C and A, respectively. The comparison of both gas-phase and solution results reveals an important feature. The hydrogen bond formation between uncharged groups is unfavored since the ∆∆Gsolv is larger than the gas-phase hydrogen bond energy. A similar situation occurs for the hydrogen bond involving a charged group. However,

the positive charge of the complex involving a charged group induces a strong polarization effect, which makes that a small increment in the equilibrium distance with respect to the gasphase value can easily give a more favorable interaction. These results led to important conclusions regarding the role of the electrostatic interactions in the drug-DNA recognition mechanisms. Structural Analysis. In order to compare the hydrogen bonding geometry of the calculated models with that determined experimentally in the crystal structure, we have analyzed the CdO(thymine/cytosine)‚‚‚HsN(ligand) hydrogen bonding geometry for two cases: amide groups (type I) and charged groups (type II). Table 6 shows the statistics for I and II types. Although in both cases the heavy-atom hydrogen bonding geometries have reasonabe values, when hydrogens are placed in their ideal positions the hydrogen bonding parameters become much worse. In particular, the average CdO‚‚‚H distances and angles (2.4-2.5 Å and 139-141°) are far from their expected values in ideal hydrogen bonding geometry. In view of these results we have re-examined all the structures looking for reasons for such a deviation. Two major features of the drugDNA hydrogen bonding interactions become apparent: (1) most of the interactions implicate bifurcated hydrogen bonds and therefore the bonding parameters deviate significatively from the linear or close-to-linear geometry expected in ideal hydrogen bonds; (2) in those cases in which the drug needs to comprise its position in order to fulfill several competing hydrogen bonds, there seems to be a preference to obtain a better geometry for

11486 J. Phys. Chem., Vol. 100, No. 27, 1996 those involving charged groups, hence sacrifying the geometry of amide groups nearby. Even if only distance criteria are applied (cutoff 3.5 Å), we have found that 18 out of 21 possible hydrogen bonds (85.7%) are formed when the involved group in the drug is a charged one, whereas only 13 out of 19 (63%) are made when drug amide groups are considered. Those observations reinforce the accepted idea of a leading role for the charged groups in the binding mechanism of the drugs to the DNA minor groove. Discussion The results obtained in this work indicate that the interaction thymine‚‚‚iminium is significantly stronger than the interaction thymine‚‚‚amide. It is instructive to compare the best estimates of ∆H298 found for both types of interactions. The net enthalpies of association computed at the HF/6-31++G**//HF/6-31G* level are -21.60 and -3.36 for the interactions involving charged and noncharged groups, respectively. These results indicate that the interaction between DNA and the charged ends of a drug are at least six times as strong as the hydrogen bond between DNA and the amide groups of the drug. Results are quite consistent with experimental data, since it is known that SN1871 bonds DNA efficiently, notwithstanding the absence of hydrogen bonds with bases.27,28 The effect of the solvent on the electrostatic and hydrogen bonding interactions has been studied by using two SCRF methods, MST/AM1 and MST/6-31G*. The accurate representation of the solvation effect is problematic because the usual force-field techniques, like Monte Carlo and molecular dynamics, are not able to describe the effect of polar solvent on the intrinsic reactivity of the solute. The SCRF procedures, implemented within quantum mechanical algortihms, have proved to be very powerful in representing the solute-solvent polarization. Thus, the average errors in the free energies of hydration obtained from MST/6-31G* were less than 1 kcal/ mol for neutral molecules37 and 5.7 kcal/mol for monovalent cations.59 The absolute average error for monovalent cations was larger than that for neutral molecules due basically to the magnitude of the absolute free energies of hydration of cations, which make it difficult to obtain small absolute errors. Indeed, the average relative error in monovalent cations is less that 7%. Estimates of the free energies of hydration indicate that amide hydrogen bonds are unfavored in water. However, the interaction involving a charged group would be stable if a small increment of the equilibrium distance with respect to the gasphase value is considered. These results suggest that interactions involving charged groups have a crucial role at the first steps of the drug-DNA recognition mechanisms. It is gerenally believed that the primary mode of recognition and binding between basic minor-groove drugs like netropsin and AT-rich regions of DNA is via hydrogen bonding rather than simple electrostatic interactions. This hypothesis is supported by the stabilization of the DNA secondary structure and the amide‚‚‚DNA hydrogen bonding interactions identified in crystallographic structures. However, our calculated interaction energies are not in total agreement with such a statement. Given the fact that hydrogen bonds between amide groups are, due to the presence of water, the less favored interactions in the first steps of the recognition mechanism between netropsin and DNA, it seems appropiate to think that the approach of the drug to DNA occurs through interactions involving charged groups. These would be between the charged centers of the drug and both the acceptor groups of the DNA and the phospate negative charges. It is also likely that the strength of the ionic interactions with phospate groups are damped by the placement of counter-

Alema´n et al. ions. Contrary to common belief, in the step of approach the hydrogen bonds would be rather passive in the sense that the hydrogen bond is a hydrophobic interaction. This hyphotesis is experimentally supported by the results of SN 18071 which is unable to form such interactions. Conclusions Quantum mechanical calculations on three model systems taken from an oligonucleotide‚‚‚netropsin complex show that interactions involving a charged group are energetically more favored than interactions between non-charged amide groups. Furthermore, calculations considering an aqueous solvent indicate the hydrophobicity of the hydrogen bond between noncharged groups, whereas interactions involving a charged group are energetically possible if a small change in the intermolecular geometry is considered. The intermolecular geometric parameters found for the charged complex from these calculations are in very good agreement with the experimental distributions, suggesting the leading role of the electrostatic interactions on the drug‚‚‚DNA recognition and binding mechanism in its earlier steps. Acknowledgment. This work was supported by the Centre de Supercomputacio´ de Cataluya (CESCA). We are indebted to Drs. Modesto Orozco and F. Javier Luque for making available to us their computer programs. References and Notes (1) Brown, D. G.; Sanderson, M. R.; Skelly, J. V.; Jenkins, T. C.; Brown, T.; Garman, E.; Stuart, D. I.; Neidle, S. EMBO J. 1990, 9, 1329. (2) Brown, D. G.; Sanderson, M. R.; Garman, E.; Neidle, S. J. Mol. Biol. 1992, 226, 481. (3) Coll, M.; Frederick, C. A.; Wang, A. H.-J.; Rich, A. Proc. Natl. Acad. Sci. U.S.A. 1987, 84, 8385. (4) Coll, M.; Aymami, J.; van der Marel, G. A.; van Boom, J. H.; Rich, A.; Wang, A. H.-J. Biochemistry 1989, 28, 310. (5) Larsen, T. A.; Goodsell, D. S.; Cascio, D.; Grzeskowiak, K.; Dickerson, R. E. J. Biomol. Struct. Dyn. 1989, 7, 477. (6) Kennard, O.; Hunter, W. N. Quart. ReV. Biophys. 1989, 22, 327. (7) Kopka, M. L.; Yoon, C.; Goodsell, D.; Pjura, P.; Dickerson, R. E. Proc. Natl. Acad. Sci. U.S.A. 1985, 82, 1376. (8) Nunn, C. M.; Jenkins, T. C.; Neidle, S. Biochemistry 1993, 32, 13838. (9) Pjura, P. E.; Grzeskoviak, K.; Dickerson, R. E. J. Mol. Biol. 1987, 197, 257. (10) Quintana, J. R.; Lipanov, A. A.; Dickerson, R. E. Biochemistry 1991, 30, 10294. (11) Sriram, M.; van der Marel, G. A.; Roelen, H. L. P. F.; van Boom, J. H.; Wang, A. H.-J. EMBO J. 1992, 11, 225. (12) Tabernero, L.; Verdaguer, N.; Coll, M.; Fita, I.; van der Marel, G. A.; van Boom, J. H.; Rich, A.; Aymami, J. Biochemistry 1993, 32, 8403. (13) Teng, M.; Usmann, N.; Frederick, C. A.; Wang, A. H.-J. Nucleics Acids Res. 1988, 16, 2671. (14) Vega, M. C.; Garcı´a-Sa´ez, I.; Aymamı´, J.; Eritja, R.; van der Marel, G. A.; van Boom, J. H.; Rich, A.; Coll, M. Eur. J. Biochem. 1994, 222, 721. (15) Jenkins, T. C.; Lane, A. N.; Neidle, S.; Brown, D. G. Eur. J. Biochem. 1993, 213, 1175. (16) Lane, A. N.; Jenkins, T. C.; Brown, T.; Neidle, S. Biochemistry 1991, 30, 1372. (17) Klevit, R. E.; Wemmer, D. E.; Reid, B. R. Biochemistry 1986, 25, 3296. (18) Liquier, J.; Mchami, A.; Taillandier, E. J. Biomol. Struct. Dyn. 1989, 7, 119. (19) Patel, D. J. Proc. Natl. Acad. Sci. U.S.A. 1982, 79, 6424. (20) Portugal, J.; Waring, M. J. Eur. J. Biochem. 1987, 167, 281. (21) Taylor, J. S.; Schultz, P. G.; Dervan, P. B. Tetrahedrom 1984, 40, 457. (22) Berman, H. M.; Neidle, S.; Zimmer, C.; Thurm, H. Biochem. Biophys. Acta 1979, 561, 124. (23) Caldwell, J.; Kollman, P. A. Biopolymers 1986, 25, 249. (24) Wartell, R. M.; Larson, J. E.; wells, R. D. J. Biol. Chem. 1974, 249, 6719. (25) Zakrzewska, K.; Lavery, R.; Pullman, B. Nucleic Acids. Res. 1984, 11, 8825.

Drug-DNA Recognition Mechanism (26) Zimmer, Ch. In Progress in nucleic acids research and molecular biology; Cohn, W. E., Ed.; Academic Press: New York, 1975; p 285. (27) Baguley, B. C. Mol. Cel. Biochem. 1982, 43, 167. (28) Braithwaite, A. W.; Baguley, B. C. Biochemistry 1980, 19, 1011. (29) Momany, F. A. J. Phys.Chem. 1978, 82, 592. (30) Hariharam, P. C.; Pople, J. A. Theor. Chim. Acta 1973, 28, 203. (31) Connolly, M. QCPE Bull. 1981, 1, 75. (32) Binkely, J. S.; Pople, J. A.; Hehre, W. J. J. Am. Chem. Soc. 1980, 102, 939. (33) Boys, S. F.; Bernardi, F. Mol. Phys. 1970, 19, 553. (34) Miertus, S.; Scrocco, E.; Tomasi, J. Chem. Phys. 1981, 55, 117. (35) Miertus, S.; Tomasi, J. Chem. Phys. 1982, 65, 239. (36) Luque, F. J.; Negre, M.; Orozco, M. J. Phys. Chem. 1993, 97, 4386. (37) Orozco, M.; Jorgensen, W. L.; Luque, F. J. J. Comp. Chem. 1993, 14, 1498. (38) Pierotti, R. A. Chem. ReV. 1976, 76, 717. (39) Floris, F. M.; Tomasi, J.; Pascual-Ahuir, J. L. J. Comp. Chem. 1991, 12, 784. (40) Negre, M.; Orozco, M.; Luque, F. J. Chem. Phys. Lett. 1992, 196, 27. (41) Pascual-Ahuir, J. L.; Silla, E.; Tomasi, J.; Bonaccorsi, R. J. Comp. Chem. 1987, 8, 778. (42) Gao, J.; Luque, F. J.; Orozco, M. J. Chem. Phys 1993, 98, 2975. (43) Bachs, M.; Luque, F. J.; Orozco, M. J. Comp. Chem. 1994, 15, 446. (44) Berman, H. M.; Olson, W. K.; Beveridge, D. L.; Westbrook, J.; Gelbin, A.; Demeny, T.; Hsieh, S.-H.; Srinivasan, A. R.; Scheneider, B. Biophys. J. 1992, 63, 751. (45) Sack, J. S. J. Mol. Graphics 1988, 6, 244.

J. Phys. Chem., Vol. 100, No. 27, 1996 11487 (46) Frisch, M.; Trucks, G. W.; Head-Gordon, M.; Gill, P. M. W.; Wong, M. W.; Foresman, J. B.; Johnson, B. G.; Schlegel, H. B.; Robb, M. A.; Replogle, E. S.; Gomperts, R.; Andres, J. L.; Raghavachari, K.; Binkley, J. S.; Gonza´lez, C.; Martin, R. L.; Fox, D. J.; Defrees, D. J.; Baker, J.; Stewart, J. J. P.; Pople, J. A. Gaussian-92, Revision F.4, Gaussian Inc., Pittsburgh, PA, 1992. (47) Cammi, R.; Bonaccorsi, R.; Tomasi, J. 1987. Modified version of MonsterGauss. Peterson, M.; Poirier, R. University of Toronto, Ontario, Canada. (48) Stewart, J. J. P. QCPE Bull. 1983, 4, 109. (49) Luque, F. J.; Orozco, M. MOPAC/CMS computer program. University of Barcelona, 1992. (50) Orozco, M.; Luque, F. J. MOPETE/MOPFIT computer programs. University of Barcelona, 1992. (51) Bru¨nger, A. T.; Kuriyan, J.; Karplus, M. Science 1987, 235, 458. (52) Jenkins, T. C.; Parrick, J.; Porssa, M. Anti-Cancer Drug Des. 1994, 9, 477. (53) Hrouda, V.; Floriani, J.; Hobza, P. J. Phys. Chem. 1993, 97, 1542. (54) Gould, I. R.; Kollman, P. A. J. Am. Chem. Soc. 1994, 116, 2493. (55) Neusheuser, T.; Hess, B. A.; Reutel, C.; Weber, E. J. Phys. Chem. 1994, 98, 6459. (56) Novoa, J. J.; Planas, M.; Whangbo, M.-H. Chem. Phys. Lett. 1994, 225, 240. (57) Sponer, J.; Hobza, P. J. Am. Chem. Soc. 1994, 116, 709. (58) Orozco, M.; Luque, F. J. J. Am. Chem. Soc. 1995, 117, 1378. (59) Alema´n, C.; Navarro, E.; Puiggalı´, J. J. Org. Chem. 1995, 60, 6135.

JP960126B