Influence of conventional hydrogen bonds in the intercalation of phenanthroline derivatives with DNA: The important role of the sugar and phosphate backbone

Abstract The influence of hydrogen bonds in model intercalated systems between guanine‐cytosine and adenine‐thymine DNA base pairs (bps) was analyzed with the popular intercalator 1,10‐phenanthroline (phen) and derivatives obtained by substitution with —OH and —NH2 groups in positions 4 and 7. Semiempirical and Density Functional Theory (DFT) methods were used both including dispersion effects: PM6‐DH2, M06‐2X and B3LYP‐D3 along with the recently developed near linear‐scaling coupled cluster method DLPNO‐CCSD(T) for benchmark calculations. Our results given by QTAIM and non‐covalent interaction analysis confirmed the existence of hydrogen bonds created by —OH and —NH2. The trends in the energy decomposition analysis for the interaction energy, ΔE int, showed that the ΔE elstat contributions are equal or even a little bit higher than the values for ΔE disp. Such important ΔE elstat attractive contribution comes mainly from the conventional hydrogen bonds formed by —OH and —NH2 functional groups with DNA not only with bps but specially with the sugar and phosphate backbone. This behavior is very different from that of phen and other classical intercalators that cannot form conventional hydrogen bonds, where the ΔE disp is the most important attractive contribution to the ΔE int. The inclusion of explicit water molecules in molecular dynamics simulations showed, as a general trend, that the hydrogen bonds with the bps disappear during the simulations but those with the sugar and phosphate backbone remain in time, which highlights the important role of the sugar and phosphate backbone in the stabilization of these systems.


| INTRODUCTION
The antitumoral effects of some metal complexes have arose a great interest for treatment of cancer disease. [1][2][3] The cytotoxic properties of such complexes have allowed to develop several effective treatments against certain types of cancer as breast, testicular, or ovarian cancers. [1][2][3] However, some of these treatments remain aggressive and are systematically at the origin of side effects. Therefore, more research in this field to reduce such aggressive side effects is a scientific multidisciplinary challenge. Indeed, the origin of the cytotoxicity of these complexes is globally known, they interact with the nuclear duplex DNA and produce the inhibition of its replication/ transcription and thus, the inhibition of the multiplication of the cancer cells. 4,5 However the knowledge on how the molecule interacts with DNA, the induced structural changes that rise up when this interaction is produced, or the explanation of the perturbation of biological processes by these complexes are still dark aspects that need to be enlighten. Thus, one of the challenges of the scientific community related to this field is to discover how the ligands of metal complexes interact with DNA and how such interaction yields cytotoxic effects to better understand and improve the efficacy of chemotherapy treatments against cancer.
Simulation and experimental techniques 6,7 have allowed to categorize two families of molecules that interact with DNA (see Figure 1). That is, the groove binders and the intercalators of DNA.
The groove binders interact with the major groove, from now MG, and minor groove, from now mg, of the DNA by forming covalent bonds through cross-links or by means of weak interactions. On the other hand, the intercalators of the DNA are included between adjacent base pairs (bps). Finally, there is a last mode of interaction called insertion that occurs when some small molecule is included between 2 bps of DNA that are not matched. 6,8 In such case, any flat molecule interacts from the mg side and ejects the mismatched DNA bases and such flat molecule acts as its π-stacking replacement (see Figure 1).

The interactions of the groove binders and intercalators result in
changes in parameters of the structure of the duplex DNA, which are the origin of the therapeutic effects. 9,10 Many studies focus on the incorporation of the 1,10-phenanthroline (phen) ligand and derivatives (see Scheme 1) in metal complexes having such systems significant antitumoral activity. [11][12][13][14][15] Among the different interactions of the metal complexes with DNA, intercalation is an important binding mode 16 in the use of these metal complexes for chemotherapy treatments and this mode of interaction is affected by the planarity of the ligand, type of donor atom and metal coordination geometry. 17 Recent experimental works on the interaction of metal complexes including phen ligands with DNA have been consistent with the intercalation of this flat ligand between DNA bps. 12,18 Linear dichroism, viscosity, and NOESY and NMR spectra also showed that methylated derivatives of phen and dipyridol[3,2-a:2 0 ,3 0 -c]phenazine (dppz) intercalated between DNA bps. 14,19 Finally, latter determination of crystal structures have given very important information about the intercalation of dppz ligand of [Ru(tap) 2 (dppz)] 2+ , [Ru(bpy) 2 (dppz)] 2+ , and [Ru(phen) 2 (dppz)] 2+ between bps of duplex DNA. 8,[20][21][22] On the other hand, it has been found in the bibliography 6,[23][24][25][26] that the process of intercalation needs prior formation of some non-covalent mg binding system. Such process is governed by kinetics. That is, while mg is a very fast process (tenths of milliseconds), the intercalation process is slower (few milliseconds). Thus, the interaction between any small molecule and DNA is best described as some reversible kinetic equilibrium defining the residence time of such interaction. 27 In general, the favored direction of the equilibrium is related to the strength of the formed hydrogen bonds and/or other weak interactions. 28 The entropic and steric factors and desolvation energies have been also found to play an important role in the kinetics of the process. 6,[23][24][25][26] In this sense, the driving force for the intercalation may be related to the energy derived from the removal of the small molecule from the aqueous medium together with the intrinsic contributions to the interaction between the small molecule and DNA (Pauli, van der Waals, electrostatic and polarization). [29][30][31][32][33][34][35][36] Thus, since the cytotoxic effect of the intercalators depends on the time of residence of the drug between bps, 27 whereas mg binding has not usually cytotoxic effects, 28,37 the design of efficient drugs should aim at favoring weak interactions to stabilize the intercalated state and destabilize the mg binding. Smart changes in the intercalators by means of the addition of different functional groups in number and position should increase the residence time of the drug between bps while preventing mg binding. Therefore, because the antitumoral activity of phen complexes has been already analyzed, 11 and different studies have appeared in the bibliography dealing with the antitumoral activity of phen derivatives [38][39][40] it is worth to study the effect of substitution of phen by means of groups capable to form conventional hydrogen bonds with DNA.  41 Mukherjee et al., 23 Robertazzi et al., 42 Vargiu et al., 43 Galindo-Murillo et al., 44 Sasikala et al., 25 52 Langner et al., 53 Hill et al., 54 Ambrosek et al., 55 Biancardi et al., 56 Deepa et al., 58 or even our previous works on this topic. [30][31][32][33]35,36 Thus, taking into account the interest of using intercalators for a range of medical applications, a structural comprehension of the interactions of such drugs with DNA results important. 59 On the other hand, it is also desirable to gain insight on how to tune the activity of such drugs by means of the introduction of substituents, which are able to form conventional hydrogen bonds with DNA when they are intercalated between bps.
Thus, the aim of this work is to carry out computational calculations to analyze the effects that the inclusion of OH and NH 2 groups in phen, which are able to from conventional hydrogen bonds and strength the interaction, produce in the intercalation process by considering the so-called ring models 56 with not only two bps but also the sugar and phosphate backbone. This model is one of the most advanced models for QM calculations on this intercalating systems since the seminal studies by Bondarev et al., 45 and Rĕha et al. 46 took into account the simplest three-body models consisting of the two bases of the DNA base pair and the intercalator, whereas the socalled sandwich model 56 considers the two bps and the intercalator but does not take into account the sugar and phosphate backbone which is crucial to explain the interaction for intercalating systems that are able to form conventional hydrogen bonds as we will see in our present study. Thus, we will show in our work that consideration of the seminal three-body models and even sandwich models to perform QM studies to analyze the electronic structures and the behavior of intercalators that may produce conventional hydrogen bonds is not enough and would lead to erroneous conclusions. Geometric parameters and energetics were analyzed and interaction properties were described in terms of the topology of electronic density, non-covalent interaction (NCI) analysis and charge transfer, whereas the stabilization of the conventional hydrogen bonds in time were monitored by means of MD simulations. Previous computational studies have focused on the study of the interaction energy and the study of frontier orbitals of different DNA models by using intercalators in which the most important interaction to focus on was the π-π stacking between bps and the intercalator. [45][46][47][48][49][50][51][52][53][54][55][56][57][58] However, the present study aims at trying to understand how derivatization of phen, a typical studied flat intercalator, by means of functional groups that are able to form conventional hydrogen bonds with DNA ( NH 2 and OH), affects in the intercalation process and gain insight on how the incorporation of such groups will modify the cytotoxicity. We shall use 4,7-(OH) 2 phen and 4,7-(NH 2 ) 2 phen as intercalators for our study (see Scheme 1) and we will see how important is to take into account the sugar and phosphate backbone in the study because of the produced conventional hydrogen bonds with the intercalator. It means that even if our study is based on QM methods, which generally requires more reduced models than MM approaches, to analyze properties that are not possible to analyze with MM methods like molecular orbitals, QTAIM topologies and electron densities, Pauli repulsions in the interaction energy, and so forth, at least, we have to consider the ring model for a correct non-dynamic QM study of the intercalated system, being the threebody models (still used in 2011) and sandwich models, still used currently, not appropriate because they could lead to misleading conclusions. We expect that this work will help to understand how substitution of phen with OH and NH 2 can tune the stability of the intercalation and thus affects the cytotoxicity of the ligand and how important is, even for QM calculations in which it is more difficult to work with large systems, to take into consideration the S C H E M E 1 Phen and phen derivatives studied in this work sugar and phosphate backbone of the DNA models used for the study of the intercalation of small molecules that may form conventional hydrogen bonds with DNA.

| RESULTS AND DISCUSSION
We carried out semi-empirical and DFT calculations including dispersion effects (see Computational Details) for different structures of adenine-thymine/thymine-adenine (AT/TA) and guanine-cytosine/ cytosine-guanine (GC/CG) bps including the 4,7-(OH) 2 phen, 4,7-(NH 2 ) 2 phen and phen intercalators by means of ring models. 56 Moreover, we have also studied two orientations for the intercalation: through the mg and via MG (see Figure 2). Thus, in order to identify each system, the following nomenclature has been used along the work: (bps/intercalator/bps)groove (e.g., (GC/4,7-(OH) 2 phen/CG)mg would be the system corresponding to the intercalation of 4,7-(OH) 2 phen between GC/CG bps through the mg). In the following subsections, we shall discuss geometries and topological analysis of the studied structures, energetics, charge distribution and the evolution of the hydrogen bonded systems with time.

| Geometries, topological QTAIM analysis of the electron density and NCI analysis
The reader can find the cartesian coordinates of the optimized studied systems in the Supporting information. The identification of hydrogen bond interactions between the intercalators and bps was based on the QTAIM topological analysis of the electron density (ρ) 60 bonding scheme and the obtained distances for such H-bonds are shown in Figure 3 for the ring models of AT/X/TA and GC/X/CG systems via mg and MG. It must be said that a more complete QTAIM study of all the studied systems can be found in the Figures  It is observed that the distortion of geometries leads to a fewer number of π-π interactions showing the most coplanar arrangements higher numbers of π-π interactions. In both systems, GC/X/CG and AT/X/TA, the intercalation via mg yields more distorted geometries than the intercalation via MG, (the distortion is especially marked for the intercalation between AT/TA bps). It is observed that this distortion is due to the presence of interactions of NH 2 and OH groups not only with atoms contained in the bps but also, and more impor-

| Energies
The interaction energy ΔE int between fragments can be split by the EDA into several contributions: the repulsive ΔE Pauli , which takes into account the destabilizing interactions between occupied orbitals, the electrostatic contribution, ΔE elstat , which is the classical electrostatic interaction between the unperturbed charge distributions of the rigid fragments and the orbital interaction contribution, ΔE orb , which accounts for the charge transfer and polarization contributions. Moreover, if an explicit correction term for dispersion interactions is employed, such extra term, ΔE disp , appears in the scheme.
Such EDA was carried out with the ADF software by means of  between the intercalators and DNA appeared. However, in these CH/π, CH/n, and HÁÁÁH interactions the dispersive nature of such weak interactions rules the interaction, whereas the electrostatic compound plays minor role. [68][69][70][71][72][73][74][75][76][77][78][79][80][81] Now, having OH and NH 2 functional groups in the phen intercalator, conventional hydrogen bonds (between hard acids and hard bases) appear in which the electrostatic forces play a major role in the interaction 82,83 and it is reflected in the EDA with an increase of the ΔE elstat contribution. On the other hand, we observe another characteristic behavior for these systems with OH and NH 2 in the ΔE orb contribution, which becomes more important than for the unsubstituted phen. Because ΔE orb is associated to polarization and charge transfer processes, these more negative values of ΔE orb for the systems with OH and NH 2 can be attributed again to the conventional hydrogen bonds formed between the 4,7-(OH) 2 phen and 4,7-(NH 2 ) 2 phen ligands and the DNA, in which a more important charge transfer is produced between the intercalator and the DNA by means of these conventional hydrogen bonds as we will see below.
Taking into account all the contributions, in general, the total interaction energy (ΔE int ) gets reinforced after the substitution with To gain insight about the effect of the solvent in the intercalation process, the contribution of the desolvation penalty, ΔE Solv , was considered and added to the total interaction energy, ΔE int , from EDA to obtain the final energy in solution, ΔE aq (see Tables 1 and 2).
We observe that for the AT/X/TA systems (Table 1)  considering desolvation penalty is À15.5 kcal mol À1 for the MG. In this case the trends when taking into account solvent effects are not the same than without consideration of solvent effects (see Figure 5).
It means that even though the ΔE elstat achieves an important role, which is reflected in the ΔE int , for these intercalators forming conventional hydrogen bonds, the solvent effects by means of the ΔE Solv penalty could change the final trends. Therefore, the main conclusion

| Analysis of the charge distribution
We can go one step beyond on the analysis of the electronic structure with the information extracted from charge distribution of the studied systems. Table 3 shows, for the studied systems, the charge distribu-   The main conclusion that arises from this MD simulation is that the π-π stacking of the successive bps avoid the deformation observed above with the smaller ring systems where the solvent molecules interacted with the hydrogen bonds of the confronted bps.
Attending to the evolution of the hydrogen bonds for the intercalated4,7-(OH) 2 phen in the tetramer DNA structure via mg, we observe that the hydrogen bond formed with the phosphate group remains unchanged along all the simulation (see Figure 11 and movie in the Supporting information). On the other hand, the other OH group that is located between the two bps presents a different behav- process with the eventual improvement of their cytotoxicity. We considered two possibilities for the orientation of the intercalation: (1) through the mg and (2)

| COMPUTATIONAL DETAILS
In order to build our models, first of all, we considered the 4e1u structure from the Protein Data Bank (PDB). Because we are interested only in the intercalation, for this 4e1u structure we cleaned and removed all the atoms with the exception of the dppz intercalator and the surrounding atoms forming the intercalation pocket. That is, the GC/CG bps, sugars and phosphates. Sodium atoms were also added in order to balance the negative charges of the two phosphate groups and the dppz ligand was transformed manually into 4,7-(OH) 2 phen and the ring model 56 (7). This process was repeated until complete the eight conformations for the OH (4) and it gave us a total of 64 conformations associated to the C C O H rotations of the two OH functional groups (see Figure S1). Because a similar situation is found for the C C N H torsion angle of the amino group, the same process of systematic conformational search was applied for the 4,7-(NH 2 ) 2 phen ligand. Subsequently, all the structures obtained from the systematic conformational search were optimized at PM6-DH2 level. 97 We chose this semiempirical method because it takes into account corrections to dispersion forces and it was observed in the bibliography and in our previous works that it describes well this kind of systems in which dispersion forces are important. 30,33,35,98 Once we optimized all the conformations at PM6-DH2 level we just took into account the most stable structure of each studied system forming several hydrogen bonds with DNA for our subsequent analyses. All semiempirical calculations were carried out with the MOPAC software. 99 Net atomic charges were calculated by taking into account the Hirshfeld scheme, 100 which was observed to perform correctly for this kind of systems. 32 Localization of BCPs 60 to perform QTAIM topologies and the NCI 65 analyses were carried out by using the AIMALL software package. 101 Hirshfeld charges and the wave function to perform the topological calculations were generated at M06-2X/6-31 + G(d,p) level from the optimized PM6-DH2 geometries with Gaussian software. 102 The so-called EDA 103 for the intercalated systems was also carried out with the ADF software [104][105][106] in order to calculate the contribution of the different kind of energies to the interaction energy between the considered fragments: intercalator and intercalation pocket (bps, sugars, phosphates and sodium atoms). The EDA splits the interaction energy into contributions related to orbital, Pauli, and electrostatic contributions following a Morokuma-type energy decomposition method. 103,107 Because it is more useful for the discussion to have an explicit ΔE disp term for dispersion the results of the EDA are reported with the B3LYP-D3 functional with the explicit Grimme's D3 correction to dispersion contribution [108][109][110][111] and the TZP basis set. Solvent effects were described by the continuous solvent model COSMO 112 as implemented in ADF software. [104][105][106] The interaction energies obtained from the EDA were compared with benchmark DLPNO-CCSD(T) calculations in which the cc-pVTZ basis set was used and the RIJCOSX approximation for coulomb integrals and HF exchange along with the cc-pVTZ/C and def2/J auxiliary basis sets.
These benchmark DLPNO-CCSD(T) calculations were performed with ORCA software. [113][114][115] We also determined the stability of the conventional hydrogen the frame of ARAID researchers. Adrià Gil is also grateful to Prof.
Emilio Artacho for fruitful discussions on this work and to Prof. Maria José Calhorda also for fruitful discussions on this work and for the years he spent in Prof. Clahorda's lab.

DATA AVAILABILITY STATEMENT
The data that supports the finding of this study are available in the Supporting information of this article.