A theoretical study of methylation and CH/π interactions in DNA intercalation: methylated 1,10-phenanthroline in adenine–thymine base pairs

The study of CH/p interactions in DNA model intercalated systems has been carried out with the popular intercalator 1,10-phenanthroline (phen) and several methyl derivatives, changing number and position, and the adenine–thymine tetramers (AT/TA) where thymine also contains a methyl group. Density Functional Theory (DFT) was used for the calculations, by means of improved functionals including dispersion effects. Our results given by the AIM analysis confirm the existence of these CH/p interactions and the energy decomposition analysis shows a perfect direct correlation between the number of CH/p interactions found and their DEint. Moreover, despite the important role of dispersion energy in the systems with more methyl groups, it is not yet enough to compensate the Pauli repulsion term, DEPauli, and the orbital contribution, DEorb, and the electrostatic contribution, DEelstat, become crucial for the stabilization of the structures in the intercalation process.


Introduction
Despite being among the weakest interactions in the family of hydrogen bonds and in a grey area considered by some authors to be between the hydrogen bond and the London van der Waals interaction, 1 the CH/p interaction 2,3 has been shown to be important in several chemical and biochemical processes. For the last three decades, several experimental 3 and theoretical studies, 4-10 which support the existence of such an attractive molecular force, have been reported. Being a "donor-acceptor" interaction between a so acid (CH) and a so base (psystem), 2,11 it can be considered a hydrogen bond, as demonstrated by experiments on the electronic effect of a substituent on stereoselectivity, [12][13][14][15] conformational equilibrium, 16 crystal packing, 17 enantioselectivity, 18,19 and coordination chemistry. 12 Moreover, statistical analyses of crystal structures show that the C-H bond tries to point toward the p system. 20 Ab initio calculations of the benzene clusters with small hydrocarbon molecules 8,[21][22][23] conrmed this preference and the role of weak electrostatic interactions in stabilizing this orientation. Nevertheless, while electrostatic forces play a major role in conventional hydrogen bonds, 24,25 the CH/p interaction has been mainly attributed to the dispersion forces, 8,20,[26][27][28][29][30][31][32][33][34][35][36][37] namely when aliphatic or aromatic CH groups are involved. Electrostatic contributions may become more relevant for species involving strong CH donors such as chloroform or acetylenic CH. Thus, the interaction energy depends on the nature of the molecular fragments, CH as well as p-groups: the stronger the proton donating ability of the CH group, the larger the stabilizing effect. The dual nature of the interaction, with electrostatic and dispersion terms, is the basis for the ubiquitous existence of this molecular force in chemistry, responsive to surrounding circumstances, especially when it operates cooperatively.
Directionality is a requisite for the CH/p hydrogen bond, distinguishing it from the classic London dispersion force. Nevertheless, in general, the CH/p interaction is a very weak interaction and its limited directionality comes from dipole/ quadrupole and charge-transfer interactions, 38 while the conventional hydrogen bonds are stronger interactions with strong directionality. 39 Dependence of an interacting system oen follows the order of the strength: the stronger the bond, the stronger the trend for the linearity. 39,40 Thus, the different nature shows that the roles of CH/p interactions in controlling structures of molecular assemblies should not be discussed in complete analogy with the hydrogen bond.
As said above, these CH/p interactions become important in several chemical and biochemical systems and several experimental and theoretical studies showed their key role for controlling crystal packing, 20,41-46 organic reactions, 47,48 conformational analysis, 49,50 and molecular recognition processes [51][52][53][54][55][56][57][58][59][60][61] mainly owing to their cooperative effect, and their importance in chemistry has been well recognized. 38,62 Moreover, CH/p interactions are also important in the stability of biological structures 63,64 because they include dispersion and polarization contributions and persist in both highly polar and apolar environments such as those found in the interior of proteins. 65 This is of paramount importance in the consideration of the effect in biochemistry. Thus, some authors highlighted the role of CH/p interactions in the stability of proteins, 66-68 the binding of carbohydrates to proteins, 69,70 and the packing of the adenine ring 71 or guanine nucleotide 72 in protein structures. In addition, CH/p contacts of the methyl group of thymine or the backbone sugars with DNA or proteins have been also described. [73][74][75] Moreover, several theoretical calculations 8,22,23,31,[33][34][35]37 and gas phase spectroscopic measurements of the CH/p interactions have been reported 22,23,33,35,[76][77][78] and it is expected that such CH/p interactions are able to play an important role in drug design.
The aim of this study is the analysis of the CH/p interactions between different methylated derivatives, in number and position, of the 1,10-phenanthroline (phen) intercalator drug (see Scheme 1) and the adenine-thymine (AT) base pairs, taking into account that thymine also contains a methyl group, which is able to interact with the intercalators by means of CH/p interactions.
We focused our study on phen methylated derivatives because it was shown that both phen and its organometallic complexes have signicant antitumoral activity. [79][80][81][82][83] In particular, its Mo(II) complexes were effective in vitro against different human tumor cell lines. 79 Viscosity, linear dichroism, NMR and NOESY spectra experiments showed that cytotoxic phen methylated derivatives intercalate between DNA base pairs and the efficiency depends on the number and position of methyl groups. 84,85 Several reviews on drugs based on intercalation have appeared recently, [86][87][88][89][90] showing how the growing interest in the use of these systems for a range of medical applications prompted a renewed need for a structural understanding of their interactions with bases of DNA. 91 The development of a strategy to quantify this process and calculate the relevant energy terms would be very useful. In this work we carry out the study of the intercalation of methylated phen systems between AT/TA stacked base pairs of DNA by DFT methods with selected functionals that take into account dispersion forces, which are crucial for the study of the intercalation and CH/p interactions, to calculate the contributions to the energy and the nature of the interaction between the active cytotoxic intercalator and base pairs.
As far as we know, this is the rst study addressing the important role of CH/p interactions in intercalated systems by means of DFT methods including dispersion contribution. The knowledge of intrinsic molecular interactions constitutes an important prerequisite to understand the role of intercalators between nucleobases in DNA structure. Therefore, we expect that this work will help to understand how methylated phen ligands behave in biological processes of intercalation, and to interpret the results of several experimental techniques.

Computational details
Several approaches reported in the literature may be used to model the intercalation of ligands in DNA.  Large systems (DNA chains from decamers to pentadecamers including the intercalators) have been modelled at low level of theory by using force elds in molecular mechanics (MM) or classical molecular dynamics (MD). 92,93,[95][96][97][98]113 On the other hand, small systems (models consisting in just one base pair and the intercalatorthree-body model) have been treated at high level of theory by means of HF, MP2 or DFT level of calculation. [99][100][101][102][103][104][105] Finally, other approximations have been found in the bibliography between the treatment of a DNA sequence of various steps and the intercalator at low level of theory and the treatment of the three body systems with the base pair and the intercalator at high level of theory. That is, the use of a model consisting of the intercalator and two base pairs (either includingring modelor not consideringsandwich modelthe phosphate backbone and sugars) at HF, MP2 and dispersion corrected DFT levels of theory, at DF-SAPT0 level or at QM/MM level. 94,[106][107][108][109][110][111][112] Here we shall use these last models with two base pairs and the intercalator and DFT methods including dispersion corrections.
To build the models, we used the optimized geometries of AT/TA (q ¼ 36 ) described in our previous work, 114 which correspond to Watson-Crick base pairing described in DNA. For such structures, we increased the distance between the two base pairs (z coordinate) to twice the initial value. Then the intercalator was inserted manually to be equidistant from both base pairs and to achieve maximum overlap with them. Only two different orientations were taken into account for the model systems (intercalator + two DNA base pairs). In one, the intercalation takes place from the minor groove. This situation is reproduced when the N 1 and N 10 atoms of the intercalator (Scheme 1) are close to the minor groove. The localization of the side corresponding to the minor groove in the base pairs is shown in Fig. 1. In the other orientation, the intercalation takes place from the major groove side (N 1 and N 10 , see Scheme 1, close to the major groove; see Fig. 1 for the localization of the side corresponding to the major groove in the base pairs). The structures corresponding to these two orientations were called (AT/intercalator/TA)mg (intercalation from minor groove) and (AT/intercalator/TA)MG (intercalation from major groove). The intercalator can be phen, 4-Mephen, 5-Mephen, 4,7-Me 2 phen, 5,6-Me 2 phen or 3,4,7,8-Me 4 phen. Full geometry optimizations were performed using a DFT approach 115 with the M06-2X functional, 116,117 which includes some dispersion effects, and the 6-31+G(d,p) basis set. 118 We already used with success this functional to study the stacking of AT/TA base pairs and the role of CH/p interactions in a previous study. 114 Moreover M06-2X functional and its precursor M05-2X were used by other authors to study intercalated systems 110,112 to avoid the use of the more time consuming MP2 method, which had been extensively used 99,100,[102][103][104][105][106][107]110 before the appearance of DFT methods including dispersion. Because we are not studying excitation/ decay processes as uorescence, phosphorescence, etc. we do not take into account multireference methods. In fact, we tested the stability of the wavefunction for different systems and the wavefunction was stable at our M06-2X/6-31+G(d,p) level of calculation.
The optimized geometries were characterized by the harmonic vibrational frequencies to verify that all frequencies are real and correspond to minima structures in the potential energy surface. The counterpoise (CP) method 119,120 was used to estimate the basis set superposition error (BSSE) in some calculations. These calculations were performed with the GAUSSIAN 09, Revision A.02 soware. 121 Net atomic charges were obtained using the natural population analysis of Reed et al. 122,123 by means of the NBO 5.0 soware. 124 The localization of bond critical points 125 and topological analysis were carried out using AIM2000 program. 126 The Energy Decomposition Analysis (EDA) to partition the interaction energy between the fragments was performed using the ADF program. [127][128][129] In this analysis, the interaction energy is divided into orbital, Pauli and electrostatic terms following a Morokuma-type energy decomposition method. 130,131 These calculations were also performed with the M06-2X functional with an uncontracted triple-z basis set of Slater type orbitals augmented with a polarization function (TZP). Since the calculations with M06-2X functional gave some convergence problems, we also carried out the EDA with the B3LYP-D3 functional with the explicit Grimme's D3 correction to dispersion forces, 132-135 with the same basis set (see ESI †). Previous calculations on the intercalation of phen between base pairs showed that when comparing the EDAs obtained with M06-2X and B3LYP-D3 functionals the conclusions are very similar. 94 Moreover, we think that the discussion becomes simpler with an explicit term for dispersion (DE disp ).

Results and discussion
In the next sections we shall analyze the effect of methylation in phen, namely number and position of methyl groups, when intercalation occurs between AT/TA base pairs. First, we focus on the AIM analysis of the optimized geometries of the systems. We then analyze the energetics emphasizing the changes in the stability and in several energetic contributions. Finally, a global discussion of the results is given.

Geometries and AIM analysis
The structures of all the systems optimized at M06-2X/6-31+G(d,p) level are shown in Fig. 2  Several CH/p interactions are observed aer optimization of the structures and for some of them, as (AT/phen/TA)mg and (AT/4-Mephen/TA)mg, reorganize drastically the geometry. In the case of (AT/phen/TA)mg system the adenine in the lower pair is extruded in order to form some CH/p interaction with the methyl group of the upper thymine and, in the case of (AT/4-Mephen/TA)mg, even folding of the adenines occurs to achieve this CH/p interaction now with the methyl group of the intercalator (see Fig. 2).
The identication of CH/p interactions was based on the AIM 125,136,137 topologies (see some examples in Fig. 4 and all the structures in the ESI † along with the values for the electronic density, r, and its laplacian (V 2 r) at each bond critical point (BCP)). The length of CH/p interactions goes from 2.33 to 3.48 A. The (AT/phen/TA)MG system has only one CH/p, whereas the system with more CH/p interactions is the (AT/3,4,7,8-Me 4phen/TA)MG system with 8 bond paths associated to the CH/p interactions (see Table 2, Fig. 4 and ESI †). As a general trend we nd more CH/p interactions for the structures where intercalation is produced from the minor groove than from the major groove. The only exception is the above mentioned case of (AT/ 3,4,7,8-Me 4 phen/TA)MG with 8 bond paths associated to CH/p interactions, whereas for the (AT/3,4,7,8-Me 4 phen/TA)mg only 5 bond paths are found.
CH/p interactions are characterized by values of r at the BCPs that go from 0.0044 a.u. to 0.0122 a.u. and a positive Laplacian (see from Tables S4-S15 of the ESI †). We found between 4 and 7 bond paths (BPs) associated to these CH/p interactions in all the systems that contain 4,7-Me 2 phen, 5,6-Me 2 phen or 3,4,7,8-Me 4 phen intercalators and 2 BPs associated to CH/p interactions in the monosubstituted systems (see Fig. 4 and ESI †). The methyl group of thymine is also involved in CH/p interactions with the intercalator. The values of r in the associated BCPs go from 0.0062 a.u. to 0.0073 a.u. with a positive  In each case, a side view (above) is presented to appreciate better the rise (R in A) displacement of base pairs, and the intercalator$$$base pair distance (in parentheses) and a top view (below) is also shown to appreciate better the twist (q in ) motion. Hydrogen bond distances between base pairs are also collected. Laplacian. These CH/p interactions will favor in general the intercalation process and they go from 1 in the case of the (AT/ phen/TA)MG to 8 for the (AT/3,4,7,8-Me 4 phen/TA)MG. The importance of stabilizing CH/p interactions was already stressed for the methane/adenine systems. 138 Energetics Table 3 collects in the second column the energies of formation, E f of each structure from the separated fragments (the two optimized base pairs and the optimized intercalator in the case of the systems with intercalator or just the two optimized base pairs when the intercalator is not included in the system). The many-body analysis of the interaction energy is also depicted in Table 3. Thus, the third column corresponds to the interaction energy of the bodies, E int(bod) , which is calculated by subtracting from the total energy the energy of the isolated fragments (4 in the case of stacked base pairs without intercalator and 5 for the systems with intercalator) with the geometry they have in the nal system. The other columns contain the analysis of the interaction energy, obtained by splitting it into different twobody and many-body interactions, all of them calculated with the geometries in the total system. E HB is the two-body hydrogen-bonding interaction in each base pair, calculated as the difference between the energy of the base pairs and the sum of the energies of the two bases. E S has the other four two-body interactions between the four bases: two intrastrand and two interstrand stacking interactions in the case of the systems without intercalator, the interstrand terms being usually repulsive and, in the case of the systems with the intercalator, the four two-body interactions between the four bases and the intercalator, all of them being attractive. Finally, E MB is the multibody interaction energy which was calculated by substracting the two-body terms (E HB and E S ) from the total interaction energy of the bodies, E int(bod) .
As expected E int(bod) are more negative than E f , because the deformation energy is not taken into account in E int(bod) . The behavior of E f differs when comparing intercalation from the minor groove with the intercalation from the major groove. In the latter case E f correlates with the number of methyl groups and the energies go from À12.7 to À20.9 kcal mol À1 . Also, the number of CH/p interactions increases (see Table 2, topologies of Fig. 4 and ESI †) with the number of methyl groups and it can be directly related to the stabilization of E f when going from phen to 3,4,7,8-Me 4 phen. On the other hand, when intercalation takes place from the minor groove, E f energies go from À11.9 to À21.4 kcal mol À1 and whereas they do not correlate with the number of methyl groups, they do with the number of CH/p interactions. This indicates that, more than by the number of methyl groups in the systems, the stabilization is determined by their capacity to form CH/p interactions. This correlation is only broken in the (AT/4-Mephen/TA)mg system, the most stabilized system but with only 4 CH/p interactions. We attribute this anomalous behavior of (AT/4-Mephen/TA)mg Table 1 Distances ( A) for the CH/p interactions found in the optimized (AT/intercalator/TA)mg systems at M06-2X/6-31+G(d,p) level System d(Y-CH 3 /Z(B)) a a Y is the fragment with the CH 3 group that produces the interaction (thymine, T, or the intercalator, I). Z is an atom of the fragment acting as p-system for which a bond path is associated, whereas B is the fragment (A for adenine, T for thymine and I for the intercalator).  a Y is the fragment with the CH 3 group that produces the interaction (thymine, T, or the intercalator, I). Z is an atom of the fragment acting as p-system for which a bond path is associated, whereas B is the fragment (A for adenine, T for thymine and I for the intercalator).
to the formation of strong N 6 -H/N 3 hydrogen bonds (see Fig. 1, 2 and S3 of ESI †) between both extruded adenines, which can neither be formed for (AT/4,7-Me 2 phen/TA)mg nor for (AT/ 3,4,7,8-Me 4 phen/TA)mg structures. These systems have also a methyl group in position 4 but the adenines are not folded enough. E int(bod) shows similar trends as E f . Both values of E f and E int(bod) show that intercalation through the minor groove is favored for disubstituted systems and monosubstituted AT/4-Mephen/TA system, whereas intercalation through the major groove is favored for the most substituted AT/3,4,7,8-Me 4 phen/ TA system and the system without methyl groups AT/phen/TA. For (AT/5-Mephen/TA) the energies are very similar and no preference is observed for any. These trends can be associated to the number of CH/p interactions achieved. The structures that form more CH/p interactions when intercalate via minor groove will have more negative E f and E int(bod) than their counterparts intercalated through the major groove (4-Mephen, 4,7-Me 2 phen, and 5,6-Me 2 phen). The 3,4,7,8-Me 4 phen structure, having the highest number of CH/p interactions for the intercalation via major groove, has the lowest E f and E int(bod) when intercalation is produced through this side. The inuence of the CP correction in E int(bod) was also calculated for these systems (see Table 3). The values go from 7.3 to 8.8 kcal mol À1 , which would be $11% of the E int(bod) . Thus, the effect of the BSSE does not change dramatically E int(bod) at our level of calculation.
Let us now discuss the partition of the E int(bod) into two-body and many-body terms (see Table 3). For all the systems there are two similar values of E HB , which correspond to hydrogen bonds of the upper and lower pair of bases and the sum is given in parentheses. As a general trend, all the values of E HB for methylated systems are very similar to those for the systems with non-methylated phen. Thus, methylation of phen does not affect the forces related to the hydrogen bonds, in agreement with the small changes observed in hydrogen bonds aer methylation (compare values in Fig. 2 and 3).
The values of E S are more negative than E f because they do not contain the energy relaxation term. Such values of the E S increase with the number of CH/p interactions and thus with the number of methyl groups when intercalation is produced via major groove. On the other hand, when intercalation takes place from the minor groove, the number of CH/p interactions is not proportional to the number of the methyl groups (see Table 1 and Fig. S3 and from Tables S4-S9 of the ESI †) but the stabilization of E S is totally proportional to the number of effective CH/p interactions (compare the number of CH/p interactions in Table 1 and the values of E S in Table 3). Thus, it can be concluded that the stabilization of E S depends on the number of CH/p interactions in the structure.
Finally, the E MB contribution is residual compared to E HB and E S with the exception of the (AT/4-Mephen/TA)mg system. For this structure the value of E MB (À10.8 kcal mol À1 ) is roughly 1/3 of the value for E HB or E S and such high value can be attributed to stabilizing N 6 -H/N 3 hydrogen bond interactions between the extruded and folded adenines. Also, E MB is cooperative (E MB < 0) for all the systems.
The energy of the interaction between the rigid fragments can also be decomposed by the EDA into several contributions: DE elstat corresponds to the classical electrostatic interaction between the unperturbed charge distributions of the rigid fragments, DE Pauli comprises the destabilizing interactions between occupied orbitals, and the orbital interaction contribution DE orb , accounts for charge transfer and polarization terms. The EDA provides another interpretation of the energy terms and was performed with ADF 127-129 (see the Computational details), using different functionals. As reported previously, 94,131,139 if an explicit correction term for dispersion interaction is employed, the results of the EDA remain unchanged and the dispersion contribution appears as an extra   Fig. 5. DE disp values are around À50 kcal mol À1 . This DE disp contribution is more important for 3,4,7,8-Me 4 phen intercalator (À55.5 kcal mol À1 and À60.1 kcal mol À1 ), not surprisingly because 3,4,7,8-Me 4 phen, with four methyl groups, has the highest value of polarizability (209.9 a.u. À3 , see Table S3 of ESI †) and thus the strongest dispersion forces. DE disp for (AT/4,7-Me 2 phen/TA)MG has also a quite different value (about À55 kcal mol À1 ) because this ligand with two methyl groups has also a high polarizability (179.0 a.u. À3 , see Table S3 of ESI †). The DE orb component is very similar for all the systems (between À12.4 and À17.0 kcal mol À1 ), while DE elstat ranges from À25.8 to À36.1 kcal mol À1 . This attractive contribution has a determining inuence in the energy balance dening the total DE int , as it accounts for roughly 1/3 of the attractive forces and parallels the behavior of DE int . This was also observed in another work 140 addressing the role of electrostatics in stacked DNA base pairs and in our previous work on the intercalation of phen. 94 This term stabilized the DE int since the DE disp contribution is not sufficient to compensate the repulsive DE Pauli term and DE elstat provides the other considerable stabilizing contribution to DE int . Also, the values of DE elstat are more negative when intercalation takes place from the major groove with the only exception of (AT/4-Mephen/TA)MG (À29.5 kcal mol À1 ), for which its counterpart (AT/4-Mephen/ TA)mg has a more negative DE elstat value (À31.2 kcal mol À1 ).
DE int has the smallest values (in absolute value) for the systems without methyl groups. On the other hand, the most negative values of DE int correspond to the 3,4,7,8-Me 4 phen intercalator for the intercalation via major groove (À38.9 kcal mol À1 ) and to the 5,6-Me 2 phen system when intercalation takes place from the minor groove (À37.7 kcal mol À1 ). Also, the disubstituted 5,6-Me 2 phen and 4,7-Me 2 phen systems have more negative values than the monosubstituted intercalators. We notice again that the stabilization of the DE int reects the number of CH/p interactions and it increases with the number of methyl groups when intercalation is produced through the major groove. On the other hand, when intercalation takes place from the minor groove there are more CH/p interactions for dimethyl systems (6 CH/p interactions, see Table 1, Fig. 4 and ESI † about AIM results) than for the tetramethyl system (5 CH/p interactions) and the DE int is more negative for dimethyl systems. Thus, in these and the other systems, there is a perfect and direct t between DE int and the number of CH/p interactions, highlighting the importance of these CH/p interactions for such biological processes.
Finally we looked at the relation between the stabilization of the DE int not only with the number of CH/p interactions but also with the strength of these CH/p interactions. We can quantify the strength of CH/p interactions by the values of the electronic density (r) at BCPs points associated to CH/p interactions. 62 That is, as the value of r in the BCP increases, the strength of the CH/p increases. Thus, we can dene the cumulative electronic density as the summation of all the values for r in each bond critical point associated to CH/p interactions of any structure and this cumulative electronic density can give Table 3 Formation energy (E f ) and multi-body analysis of the interaction energies, in kcal mol À1 , for the free stacked base pairs AT/TA and including the intercalator at M06-2X/6-31+G(d,p) level of theory us an idea about the total strength achieved by the sum the individual CH/p interactions in such structure. Thus, we plotted the cumulative electronic density r vs. the values of DE int for (AT/intercalator/TA)mg and (AT/intercalator/TA)MG systems in Fig. 6 (see ESI † for the r values). Both plots show a perfect correlation between the stabilization of DE int and the increasing value of r in the bond critical points associated to CH/p interactions.

Discussion
In this exhaustive DFT study we analyzed the geometries, the energy of the bodies, the EDA, and the CH/p interactions of methylated phen systems intercalated between AT/TA stacked base pairs, which also contain methyl groups capable to interact with the p system of the intercalator. We showed the importance of the effective CH/p interactions, triggered by the methyl groups to stabilize the nal intercalated structures. Indeed, CH/p interactions, mainly characterized by the AIM analysis and BCPs and bond paths, were shown to be present for the systems where the -CH 3 groups of the intercalator or of the thymine DNA base interacted with either p-system of the base pairs or the intercalator. We established a perfect t between the number of CH/p interactions and the stabilization of E f , E int(bod) , E S or DE int (see denition in the section corresponding to energetics). These terms reect the stabilization of the nal intercalated system. This result  is very interesting because it suggests that, more than the number of methyl groups, it is the number of the effective CH/p interactions that rules the stability of the nal intercalated system. For this reason, in the cases of (AT/4,7-Me 2 phen/TA)mg and (AT/5,6-Me 2 phen/TA)mg with 2 -CH 3 groups, these energy terms are more negative than for (AT/3 84 and Brodie et al. 85 who reported the evidence of modulation of the cytotoxicity with methylation of phen in different number and position, but offered no explanation for these experimental ndings. Here, we describe the direct correlation we found between the number of effective CH/p interactions and the stabilization of E f , E int(bod) , E S , and DE int . If the cytotoxic effect of the intercalators depends on their time of residence between base pairs, 141 then it can be enhanced by the formation of hydrogen bonds or electrostatic interactions. This was used by Snyder et al. 142 to explain the cytotoxicity of Michler's ketone, which depended signicantly on the position of the positively charged N-dimethyl groups. Thus, groups of the intercalator able of forming hydrogen bonds and/or electrostatic interactions dictating the residence time and biological activities of any drug will become essential for cytotoxicity. 143 Some authors even suggested a minimum value of electrostatic energy, setting the limit to hydrogen bond interactions yielding effective intercalation. 142 This nding is similar to those of McFayden et al. 84 and Brodie et al. 85 now for neutral methyl groups. However, in the case of methylated phen systems, because methyl groups are neither cationic nor capable of forming charge assisted strong hydrogen bonds, the idea of dening a electrostatic threshold above which the species become cytotoxic is not realistic owing to the small energy amounts involved. Here, we emphasize the importance of these CH/p interactions and notice that the higher the number of CH/ p interactions, the more stabilized E f , E int(bod) , E S and DE int will be. Thus, the higher number of CH/p interactions contribute to an increased stability of the intercalation complex and we extrapolate that the systems able to form more CH/p interactions will be the most cytotoxic. In fact, the 5,6-Me 2 phen system, which was the most cytotoxic system in the work of Brodie et al. 85 (in that case intercalated between GC/CG via minor groove) has the highest number (6) of CH/p interactions for the intercalation between AT/TA base pairs via minor groove.
We believe that quantication of intrinsic contributions to interaction energy in the intercalated systems can play an important role in favoring intercalation and thus, in discrimination the cytotoxicity of the intercalators. Of course, many other variables will govern the process of intercalation/deintercalation of ligands in DNA aside from the intrinsic forces studied here, namely, solvation, entropy and steric effects as suggested by Sasikala et al. 92 and Franco et al. 113 Nevertheless, the intrinsic interactions that take place directly in the intercalation pocket must have a primary role and their comprehension and rationalization becomes fundamental.
Because the objective of our work was to give new and detailed insight on the nature of the intrinsic interactions that rule intercalation (DE disp , DE elstat , DE orb , and DE Pauli ), we used a model in which the treatment at QM level with our technical resources was possible (two base pairs of DNA and the intercalator). Such interactions were studied in the seminal works of Bondarev et al. 99 andŘeha et al. 100 with three-body models (intercalator + one base pair) and this simplest of the models continued to be used for more than one decade 99-105 until the recent work of Hazarika et al. 104 Nevertheless, as stated by Hill et al., 108 we need the other base-pair to have a better description of the role of the intrinsic forces that rule the intercalation. We understand that the presence of the phosphate backbone would also improve our results because we would have a more realistic model, but such systems cannot be treated at QM level with our technical resources. Despite this limitation, we believe that the main conclusion of this study, without taking into account the phosphate backbone, would remain. Indeed, more than the number of methyl groups, responsible for the increase in the polarizability of the molecule and thus, the dispersion forces, it is the effective number of CH/p interactions that governs the stability of the intercalation of methylated phen between AT/TA base pairs containing -CH 3 groups in thymine. Moreover, we also believe that the systems able to form more effective CH/p interactions will be also the most cytotoxic.

Conclusions
We carried out DFT optimizations on the stacked systems AT/TA incorporating methylated phen intercalators (4-Mephen, 5-Mephen, 4,7-Me 2 phen, 5,6-Me 2 phen, and 3,4,7,8-Me 4 phen) by using the M06-2X functional, which includes the effect of dispersion, to determine the effect of intercalator methylation in the structure, energetics and bond properties of the DNA base pairs where intercalation is produced and to analyze the importance of the CH/p interactions in such systems.
The multibody energy analysis shows that the E HB contribution does not change signicantly with methylation in agreement with the lack of changes in lengths of the hydrogen bonds. On the other hand, the effect of methylation produces more negative E S contributions, and these become more important than the E HB . The E MB term can be considered residual in all the analyzed structures with the exception of the folded (AT/4-Mephen/TA)mg system for which we obtained a value of À10.8 kcal mol À1 ($14% of the E int(bod) ). This increase is due to the formation of N 6 -H/N 3 hydrogen bonds between the extruded and folded adenines and is unlikely to occur in the real systems.
The EDA supports the same conclusion we obtained in our previous work. 94 Even increasing the values of polarizability by methylating phen and thus the attractive DE disp contribution to DE int , this DE disp term is necessary but not enough to balance the repulsive DE Pauli contribution, so that DE orb and specially, DE elst contributions are crucial to achieve the stability of the intercalated system.
Finally, we found a perfect direct correlation between the number of effective CH/p interactions and the stabilization of E f , E int(bod) , E S , and DE int . It suggests that the number of effective CH/p interactions which stabilize the intercalated nal structure is more important than the number of methyl groups. In an attempt to build a bridge between our results and the cytotoxicity of the systems, we suggest that the structures in which more CH/p interactions are present will also lead to the most cytotoxic systems. We hope that our work will help to shed light on understanding these important intercalation processes, which are of the highest current interest.