metal complexes that interact with dna have gained considerable interest as diagnostic agents and chemotherapeutic drugs ( 1,2 ) .
these include complexes of transition metals from groups 812 with dozens of combinations of ligands and oxidation states ( 36 ) . among the transition metals ,
copper has also been extensively studied and proven as a promising candidate for drug development ( 79 ) .
copper toxicity has been hypothesized to come from its ability to produce reactive oxygen species ( 10 ) , to displace other metal ions , to induce lipid peroxidation ( 11 ) and/or to directly cleave nucleic acids ( 1214 ) .
experimental data suggest that these compounds interact directly with dna , however , the specific molecular interactions and modes of binding are not clearly established ( 15 ) .
the casiopenas family of copper compounds ( ccs ) have been in active research since 1980s ( 1618 ) .
these complexes have shown promising biological activity to a myriad of tumors both in vitro and in vivo ( 1923 ) ; because of this , phase i clinical trials are currently underway for two of the compounds .
the general formula is [ cu(nn)(no)]no3 and [ cu(nn)(oo)]no3 where the nn ligand denotes either 2,2-bipyridine or 1,10-phenanthroline ( the aromatic ligand ) no represent an essential amino acid or peptides and the oo represents a non - aromatic ligand either acetylacetonate or salicylaldehydate ( figure 1 ) . in vitro
, assays have been conducted for multiple cell lines ( 22,24,25 ) having a 10-fold enhanced effect when in comparison with the defacto transition - metal anticancer drug cis - platin ( 19,26,27 ) .
experimental observations using casiopenas ( 28 ) and similar copper complexes showed nuclease activity when in contact with dna ( 2931 ) similar to the activity detected in other metal complexes that include co , ni , ru , zn and rh ( 3235 ) .
the intent of studies of the specific interactions between the casiopenas family of complexes and dna is to facilitate the development of drugs with increased specificity and reduced toxic side effects .
querying the protein data bank ( pdb ) for drug - dna structures that contain copper complexes shows structures where copper ions form complexes with the nucleobases in a z - dna configuration ( 36 ) .
schultz and co - workers report a modified b - dna chain with modified residues that forms a complex with added copper ions ( 37 ) .
neidle and co - workers report a very interesting copper ( ii ) salphen complex stacked between two anti - parallel g - quadruplex chains ( 38 ) .
electron paramagnetic rresonance ( epr ) techniques have been used to study the interactions between dna chains and ( 1,10-phenanthroline)-copper(ii)-(amino acids ) complexes ( 29,39 ) .
species aligning to the epr measured g|| axis are almost parallel to the phenanthroline moiety in proximity to dna , however , it is not clear how deep the phenanthroline ring inserts into the dna double helix .
further computational studies modeling the drug - dna interactions have suggested the atomic mechanism by which complexes can interact with either the grooves of the dna or via intercalation between base pairs , and these studies have suggested information about the thermodynamic and energetic properties ( 2,40 ) .
several groups have applied molecular dynamics ( md ) , quantum mechanics ( qm ) and hybrid methodologies ( qm / mm ) to models of copper complexes binding with dna ( 4145 ) . the complex cu[1,10-phenthroline]with multiple functional groups and a serinol link between the ligands have been studied by magistrato and co - workers using md and qm and the simulations and energetic analyses suggest that these complexes bind to the dna in the minor groove with a partial intercalation between base pairs ( 43 ) with related calculations using qm methodologies yielding similar results ( 42 ) . in 2012 , we reported a md - dft - qtaim study to determine the specific site of recognition between a copper complex ( cc ) and dna ( 46 ) .
the formation of the cc  desoxyrribose - phosphate adduct in the minor groove proves to be a good candidate as initial step toward the cleaving of dna chains .
the copper atom within the complex links to an oxygen atom of a phosphate group , whereas the aromatic ligand interacts with the desoxyrribose by means of \documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } $ { \rm c } - { \rm h } \cdots { \rm \pi } , { \rm o } \cdots { \rm \pi } ( { \rm c})\;{\rm and}\;{\rm o } \cdots { \rm \pi } ( { \rm n})$\end{document } contacts . here
we expand on the previous work to more directly understand the intercalation process of a family of 21 casiopenas and to also determine the preference of the cc to remain within the minor groove versus undergo intercalation by a md study .
thanks to recent advances in biomolecular simulation methods applied to dna and access to large - scale computational resources , it is now possible to reproducibly converge the structure of dna ( 47,48 ) and also , as shown in this work , to spontaneously observe intercalation processes of cc 's in md simulations .
( a ) cu[(2,2-bipyridine)(acetylacetonate ) ] , ( b ) cu[(1,10-phenanthroline)(acetylacetonate ) ] and ( c ) cu[(1,10-phenanthroline)(glycinate ) ] .
the 21 compounds investigated in this work are listed in table 1 with representative structures shown in figure 1 .
the starting structure was based on the x - ray crystallography structures of the free cc when available .
all structures were fully optimized using density functional theory with the m052x ( 4951 ) meta - gga hybrid functional and the 6311+g(2d,2p ) basis set . starting with the optimized geometry , the molecular mechanics parameterization for each copper compound
was facilitated using the antechamber program included in the ambertools 14 software and the general amber force field ( gaff)(52 ) for small organic molecules .
atomic charges were generated using the resp methodology at the hf/631 g * level of theory with the aid of the r.e.d server ( 53 ) .
non - bonded lennard - jones parameters for the copper atom (  = 0.0427 kcal / mol , rmin/2 = 1.0330  ) were taken as reported by babu and co - workers ( 54 ) , and previous work using these sets of parameters have yielded results that match experimental structures ( 46 ) . to validate and assess
the model parameters for each compound , 10 ns of md simulation was performed on each compound ( tip3p ( 55 ) water model , 1000 steps of minimization and 5000 steps of heating , 300 k final temperature ) to ensure proper sampling of representative structures .
comparison between the available x - ray crystals , dft - optimized calculations and the average structures from the md simulations is shown in supplementary table 2s of the supporting material .
acetylacetonate ; dime  dimethyl group ; me  methyl group ; tetrame  tetra methyl group ; gly  glycinate .
the x substituent is located in the aromatic ligand ( either bipyridine or phenanthroline in all the cases ) . as a model for understanding the binding of these 21 copper compounds to dna we used the dickerson - drew dodecamer ( ddd ) with
. the high - resolution structure of the free dna was obtained from the pdb database ( pdb code : 1bna ) and md simulations were performed with the parmbsc0 ( 62 ) modifications of the cornell ff99 force field ( 63 ) .
crystallization waters and ions were removed and hydrogen atoms where added according to canonical templates included in the force field using tleap .
each copper compound was manually oriented in the minor groove of the ddd using the accelrys ds visualizer ( 64 ) .
previous work suggested that the preferred binding site for these type of compounds is in the minor groove of dna ( 46 ) .
the starting position for the copper was the same for all the systems and was placed between the atoms dt20:o1p and dt8:o1p .
the distance between these two atoms is 11.6  so that the copper atom was placed 5.8  from both the oxygen atoms and parallel to the direction of the minor groove and models were built both with the aromatic and non - aromatic ligands inside the minor groove ( figure 2 ) . with these configurations , the compounds starting structures are next to the central a6-t7 step . for compounds 04 , 05 , 09 and 1121 , all of which are non - symmetric , the starting position was in either an  or a  position : for the compounds cu[(4-me - phen)(acac ) ] , cu[(5-me - phen)(acac ) ] , cu[(5-phenyl - phen)(acac ) ] , cu[(5-cl - phen)(acac ) ] and cu[(5-no2-phen)(acac ) ]
, the  position corresponds to the substituent of the aromatic ligand placed within the minor groove and  position with the substituent placed away from the minor groove ( the molecule is flipped 180 in the z axis ) . for the compounds with the glycinate ligand ( 1321 ) ,
the  position is with the carboxylate group within the minor groove and  position with the carboxylate group placed away from the minor groove .
representative starting position for the 21 copper complexes in the minor groove of the dna , oriented between the a6t7 and t18a19 base - pair steps .
the copper atom of the complex is at a distance of 5.8  from each oxygen of the phosphate group .
the dna chain and the copper complex are in licorice representation , thymine 8 and 20 are in the ball and stick representation . joung 
cheatham na cations ( 65 ) were added to neutralize the system and the system was solvated within a truncated octahedral periodic box of 71  per side size with the tip3p water model using a 10  distance from any solute atom .
excess na and cl ( 55 ions for each system ) where added to reach a concentration of 150 mm to simulate the intracellular ionic conditions . to avoid any biasing created by the initial placement of the ions ,
the randomizeions command of cpptraj was used in which a randomly chosen water molecule is swapped for an ion and a new set of coordinates is generated for every simulation .
electrostatic interactions were calculated using the particle mesh ewald method ( 66,67 ) with a cutoff value of 9  and default parameters , a 2 fs integration time step and the shake ( 68 ) algorithm to constrain all bonds with hydrogen atoms ( tolerance set to 0.00001  ) .
initial minimization was done using 500 steps of steepest descent and 500 steps of conjugate gradient minimization with a 25 kcal / mol force restraint applied to the dna and copper - compound atoms .
afterward , heating of the system was carried out in a two - step method using the same restraint force constant on the complex , first step raising the temperature to 300 k over 5000 steps and then a 50 000-step equilibration at constant temperature .
next we gradually decreased the harmonic restraints to the complex with 1000 steps of minimization and 25000 steps of md equilibration using 5 , 4 , 3 , 2 , 1 and 0.5 kcal / mol restraint values .
the equilibration procedure was followed by 10 s of production simulation using constant temperature at 300 k ( berendsen thermostat , coupling value of 5.0)(68 ) .
a total of 35 simulations were carried out , trajectory information was saved every 1 ps of md simulation time .
clustering analysis was performed using the dbscan algorithm implemented in cpptraj on the casiopena ( cas residue ) after a root mean square ( rms ) fit to the first frame of simulation using the dna residues .
the fit created a common reference for the dna so the clustering algorithm , which was performed based on the binding molecule , was able to detect root mean square deviation ( rmsd ) differences , which generated the different binding positions .
additional values used were a value of no more than 2  between clusters and a sieve value of 100 frames .
energy interactions between the dna 12-mer and the copper complex were estimated using the mmpbsa.py script ( 69 ) .
ionic strength was set to 150 mm with a solute dielectric constant of  = 1.0 and a solvent dielectric constant of  = 80 .
five thousand frames from the most populated cluster were used for the mmpbsa binding energy calculations . selected exemplar cpptraj inputs used for the analysis
to study the free energy of the copper compound bound to the intercalation site of the dna , we used the representative structure of the most populated cluster from the clustering analysis as a starting point , which has the ligand bound in the intercalation pocket as the initiation point for umbrella sampling analysis .
the reaction coordinate was the distance between the copper atom and the four atoms around it ( n 
n for the aromatic ligand and o  o for the acetylacetonate ) and the center of mass of the binding pocket formed by the four nucleotides .
we used a distance reaction coordinate from 0.2 to 20  using 0.2 increment for a total of 100 windows with a 15.0 kcal / mol  harmonic restraint constant .
each step of the reaction coordinate was 5 ns long after an equilibration step of 1 ns . each window
the potential of mean force profile was obtained using the wham methodology ( 70,71 ) .
all the quantum calculations were performed using the g09-d.01 revision of gaussian ( 72 ) .
md simulations were run using the amber12(73 ) software and trajectory analysis was made using the development version of cpptraj v14 ( 74 ) .
the rsmd with respect to the experimental structure of the free dna is presented in supplementary figure 1s of the supporting material .
the rmsd values do not include base pairs at the edges of the dna to avoid higher rmsd values caused by terminal base - pair fraying events which has been observed in nucleic acids simulations using the amber force fields ( 47,48 ) . in figure 3 , we show five representations of the interactions found in the 35 simulations studied in this work obtained from the representative structure from the most populated cluster using clustering analysis described in the methodology section . each final configuration ( from a  e , as defined in figure 3 and table 2 ) aids with the analysis of the observed results .
the most populated cluster represents the binding conformation that is populated for the majority of time during the 10 s of md simulation .
the representative binding conformations found for each of the simulations are described in table 2 .
the bottom panel highlights close - up detail of the structure depicted in the top panel .
a : stacking on the terminal base pairs of the dna chain , b : minor - groove binding , c : intercalation with base - pair displacement / opening , d : minor - groove binding with stacking of one of the terminal frayed bases , and e : intercalation near the end of the dna chain .
structures shown are a single - structure snapshot at 10 s simulation time from complexes 05 , 12b , 10 , 09 and 01 .
the percentage denotes the amount of time spent in that particular conformation with respect to the total simulation time .
rmsd values of the highest populated structure from the clustering analysis versus the experimental structure ( pdb entry 1bna ) using all residues .
g values obtained from 5000 frames using the most populated cluster trajectory . from the full data set , 23 of the systems generated an interaction of type a were the copper compound detaches from the minor groove and diffuses along the backbone of the dna until reaching one of the chain edges where the compound stacks on the terminal base pairs .
the dna rmsd values do not increase above 0.2  with respect to the starting structure suggesting that the binding interaction does not produce any significant change in the dna chain for this type of interaction ( table 2 ) .
all of these systems explored multiple conformations and detached completely from the dna , although the prevalent structure as found by the clustering analysis was the stacking interaction at the end of the chain .
as shown in table 2 , the type a binding mode is not as energetically favored as the others .
the type b interactions tend to be more favorable than the terminal base - pair - stacked binding mode .
the type b configuration is characterized by binding of the ligands within the minor groove of the dna duplex with configurations very similar to the starting point for all the md simulations .
the compound 19 remained with both the aromatic and non - aromatic ligand facing the groove while the compound 19 shifted the aromatic ligand toward the groove , exposing the glycinate toward the solvent .
after 125 ns , the complex turned 185 on the z axis until the aromatic moiety and the carboxylic side of the glycinate again faced the inside of minor groove and remained in that configuration for the rest of the simulation ( refer to the structures in the supporting material ) .
the rmsd plot for molecule 03 shows a difference of 1.5  starting after 400 ns and lasted for the rest of the simulation .
the interaction of this molecule with dna resulted in a base - pair eversion / opening toward the major groove followed by intercalation of the copper complex in the resulting pocket ( configuration type c ) .
this interaction , which tends to display the most favorable binding affinity , will be discussed in the next section of the article .
type d configuration consists of the ligand bound to the minor groove but close to one of the ends of the chain where the compound forms a stacking interaction with one of the frayed terminal nucleotides ( residues 1 , 24 , 12 or 13 in the ddd ) . for 09 ,
the compound diffuses inside the minor groove throughout the simulation until it reaches the g10c15 position . in this site ,
the interaction lives for 80 ns and then the g12 base flips back to pair again with c13 .
compound 01 presents low rmsd fluctuations with respect to the reference frame , although the clustering analysis showed an intercalating structure near the terminal base pairs ( configuration type e ) .
the top plot presents the distance measures between the complex and residues 12 and 13 at the edge of the dna chain . at the start of the simulation ,
the compound is located in the minor groove , close to residue 6 , which gives a distance of 2224  . after only 80 ns
, the compound detaches from the minor groove and diffuses close to the dna backbone until it finds the edge of the dna and stacks above residues g12 and c13 ( similar to type a configuration in figure 3 ) .
after 1800 ns of a type a configuration , base opening ( fraying ) of residue g12 toward the minor groove causes an increase in the distance between residue and complex ( marked as b in figure 4 ) .
this causes the complex to reposition toward residue c13 , increasing the stacking interaction and lowering the distance between both molecules .
after 100 ns , the complex slides away from residue c13 and stacks with the base pairs c11 and g14 .
residue 13 flips back to stack again with the complex , but now from the opposite side ( marked as c ) .
residue 13 reforms the watson  crick ( wc ) pairing with residue 12 and the intercalation is complete .
the spike marked as d in figure 4 corresponds to a fraying event of both residues 12 and 13 that lasts for 95 ns before both residues move back into position .
the average structure for the most populated cluster is shown in supplementary figure s2 of si for each ligand molecule . for compound 11 ,
the copper complex detached from the minor groove and shifted toward the end of the dna interacting with the fraying base pairs for 3.1 s until the intercalation formed which then remained for the rest of the simulation . the process of intercalation at the terminal base pairs as observed in md simulation .
top : running average ( 500 frames ) of the distance between g12 and copper complex 01 ( black ) and the distance between c13 and the same compound 01 ( red ) .
the effect of the copper complexes and the five different types of binding to dna observed in the simulations can be summarized comparing the structure from the clustering analysis with the crystal experimental structure of the free dna ( table 2 ) .
as already observed in the plots of rmsd versus time available in the supporting material , the systems that have the highest deviation from the experimental structure are the ones in type c configuration ( 3.34.3  ) , where 2 nucleotides are flipped toward the mayor groove .
the type d configuration also presents high divergence in rmsd due to the distortion of the nucleotides present at the end of the dna chain , flipping toward the minor groove .
the rest of the simulations have values between 1.7 and 2.7  , which are expected values , consistent with previous md simulation work , since we are including the terminal base pairs at each end of the dna that have multiple terminal base - pair fraying events .
this is also expected because the copper complex detaches from the dna and forms a stacking interaction with a base pair at the end of the dna , this actually lowers the frequency of fraying , stabilizing the dna and hence lowering the rmsd difference from the crystal structure .
overall , from the 35 simulations performed , the majority ( 66% ) are type a with similar populations observed for the type b ( 9% ) , type c ( 11% ) , type d ( 6% ) and type e ( 9% ) interactions .
since the simulations were performed without any type of biasing or restraint on the dna or the ligand , the interaction where the copper complex is stacking at the end of the dna is the most common one .
this does not mean that the terminal stacking interaction is energetically favored , in fact , the calculated energy - binding values presented in this work ( table 2 ) show that type a binding has overall lower values compared to the other types of binding . the process for type a of interaction has been previously discussed and observed in similar unbiased , unrestrained drug - dna simulations ( 76,77 ) .
we can attribute this interaction as being predominant in our results mainly to three observations : ( i ) a lack of proper binding of the ligand with the minor - groove environment that causes the ligand to shift through the cavity until there is no more dna ; ( ii ) complete detachment of the ligand from the minor groove toward the solvent ; and ( iii ) steric hindrance caused by the substituent present in the aromatic ligand that causes the whole complex molecule to either shift inside the groove or loose the binding completely .
the majority of the complexes that interact in the type a configuration do not have any substituents in the 4 , 4 position ( in the case of the bipyridine aromatic ligand ) nor in the 4 , 7 position ( for the phenanthroline ligand ) .
none of the complexes with the acetylacetonate ligand ( complexes 01 to 12 ) were found in the minor groove ( type b ) , whereas three of the complexes with the glycinate ligand formed stable minor - groove binding : system 17 that has 5,6-dimethyl substituent and both systems using complex 19 which have 4,7-diphenil substituents in the aromatic ligand .
the 17 system remains bound for the whole simulation , fluctuating between the a5t20 and a6t19 base pairs and forming hydrogen - bonding interactions of the methyl group in position 5 of the phenanthroline ligand with the backbone residues and hydrogen bonding of the aromatic ligand with n3 of adenine .
the presence of the extra phenyl ring in compound 19 stabilizes the binding as the ring goes inside the minor groove , supplying the complex with extra binding sites and hydrogen - bond interactions with the dna backbone .
the system 19 , which started with the carboxylate group of the glycinate toward the minor groove made a 180 turn in the z - axis from 120 to 200 ns time so that the carboxylate group now is facing toward the solvent and the amino part toward the minor groove in the same orientation as the 19 system .
the 19 system fluctuated between t7a18 and t8a17 and only minor readjustments of the ligand within the groove were observed .
insertion and base - pair eversion have been observed experimentally by the barton group ( 7881 ) .
cytosine base pair in the central position and both residues are ejected toward the major groove as the bulky metal complex inserts in the resulting cavity .
, they discovered multiple conformational changes in a 9-mer dna , including the eversion of two guanine residues toward the major groove .
they observed after 300 ns of md simulations the same modes of binding as reported in this work ( end binding , minor - groove binding and intercalation , figure 3 ) .
this is due to the shape of the doxorubicin ligand , which aligned with the top- and bottom - matched base pairs in a perpendicular orientation , allowing one of the flipped bases to flip back into position . in our simulations , compounds 02 , 03
, 05 and 10 started in the minor groove and caused an at pairing to flip or base pair open toward the mayor groove .
the molecules 02 , 03 and 10 pushed open the a5t20 pair and complex 05 the t8a17 pair .
the copper compound moved inside the resulting pocket ( shown in figure 5 ) , which corresponds to the average structure of the highest populated cluster for each simulation using the full 10 s trajectory . in figure 5 , we present the first 500 ns of simulation for each of the simulations that generated the flipped nucleotides .
it is interesting to observe that the intercalation process starts early in the simulation ; the 05 and 10 systems start almost at similar time ( 160 ns ) , followed by 02 ( 308 ns ) and finally the 03 system ( 360 ns ) , and the intercalation process takes no longer than 20 ns except for the 10 system which lasts for more than 100 ns .
this system takes longer time to reach the final intercalated position due to the rearrangement of the two bulky phenyl groups in the minor groove .
top : average structure extracted from the most populated cluster using all the simulation frames for systems 02 , 03 , 05a and 10 .
bottom : running average ( 200 frames window ) of the rmsd values of the systems in the configuration type c. rmsd values calculated using the crystal structure as reference and central residues 2 to 11 and 14 to 23 .
the process of intercalation follows the general sequence : ( 1 ) ligand binding in the minor groove , ( 2 ) ligand positioning , ( 3 ) partial opening of base pairs and partial insertion of the ligand , ( 4 ) flipping out of both the bases and full ligand insertion in the cavity , followed by ( 5 ) reorientation of the ligand inside the cavity .
all four systems that were found to produce this type of intercalation follow this same basic process or model , with variation only on how much time a certain system spent in each of the 15 steps .
the main observed difference is with step 3 where the copper complex partially pushes the base pair until the wc hydrogen bonds break and the compound pushes the nucleobases toward the major groove .
this step happens almost instantaneously in systems 03 and 05 , where the at base pair is pushed with atoms 2 , 3 and 4 of the phenanthroline ring .
regardless of this observation , the process of intercalation is similar in each of the four systems . in previous dna 
ligand intercalation studies , the rise , roll and twist helical parameters provide a good description of the intercalation process ( 83,84 ) . in these simulations where a base pair is structurally separated and pushed toward the minor groove ,
there is a coordinated global change in the flipping nucleotides , affecting both inter- and intra - base - pair geometries .
structural parameters of the flipped nucleotide for each case mentioned above and the distance of the copper atom of the ligand ( used as the center of the molecule ) to the o4 and o2 atoms of the thymine base in the intercalation pocket are presented in supplementary figure 3s of the supporting material .
the first step of the intercalation process was the presence of the bound ligand in the minor groove , which in this set of simulations was also the starting position .
the distance between the ligand and the o4 and o2 atoms used as a set of reference points in the intercalation cavity is > 9  in each case ( supplementary figure s3 ) . in step 2 , the ligand reorients multiple times in the minor groove ( a and d in figure 6 ) . in figure 5 , we can see the amount of time the ligand explored the space inside the minor groove until moving to step 3 . until this step ,
the structural parameters in supplementary figure 3s fluctuate close to the canonical values . during step 3 , the opening angle and the stretch distance increase , in addition to the propeller angle , and
this indicates that one of the bases is being pushed ( b and e in figure 6 ) .
step 4 is the flip of both nucleotides toward the major groove and partial insertion of the ligand inside the resulting cavity . in this step ,
the rise of the previous and next base pairs that form the intercalation pocket goes down , which puts the four nucleotides closer together hence reducing the volume of the intercalation pocket .
step 5 presents the deep insertion of the ligand inside the pocket ( c and f in figure 6 ) .
once inside , the ligand shows rotation on the z - axis , reorienting itself multiple times during the simulation . in supplementary figure s4
the orientation of the ligand inside the cavity formed a 83 angle with respect to the bases shifting 11% of the total simulation time .
system 05 shifted continuously between angles of 200 and 90. finally , system 10 shifted conformation upon entering the cavity , forming a 290 and shifting conformation in a small interval of 20 for the rest of the simulation .
this stability inside the pocket for this compound is due because of the steric hindrance caused by both phenyl groups in the aromatic ligand of the complex .
orange and yellow are the top and bottom base pairs forming the intercalation cavity after base flipping .
if the intercalation  base - flipping events are truly an effect of the binding of the casiopenas complexes , simulations were performed using the ligands with no metal center , thus , only the two aromatic ligands and the two non - aromatic ligands : acetylacetonate , glycinate , bipyridine and phenanthroline .
each molecule was parameterized using the same methodology described and manually placed in the same position inside the minor groove of the dna 12-mer in the same manner as the casiopenas compounds .
each of the four systems was run with three independent copies using the same equilibration and production protocol and randomizing the position of the counter ions for each copy .
supplementary figure 5s in the supporting material shows the rmsd values for each of the copies .
visual inspection of the trajectories showed transient and short - lived ( < 50 ns ) interactions with the terminal base pairs of the dna with the aromatic ligands .
the non - aromatic ligands detach from the minor groove early in the simulation ( 25 ns ) and only short - lived interactions with the dna 's backbone are detected . considering the prolonged sampling time used ,
these results indicate that intercalation between two base pairs and intercalation with the opening of internal base pairs are unlikely to be detected using the ligands by themselves .
additional analysis was done to understand the base - pair eversion event of systems 03 and 10 where the copper complex was removed from the opening of the everted nucleotides to examine if a canonical b - dna structure could be reformed .
the starting structure was obtained using only the trajectory frames were the complex was inserted in the dna and clustered the frames using the same methodology as described previously .
a representative structure of the top populated cluster was then extracted , the copper complex was deleted and three independent simulations were run , using the same equilibration and production protocol , randomizing the counter ion positions for each copy .
supplementary figure 6s in the supporting material shows the rmsd value of the dna and distance from the center of mass for each flipped nucleotide of each copy .
system 03 reforms the wc pairing in residues a5t21 in 1820 ns of simulation time as can be observed in the distance top plot in supplementary figure 6s . in two of the simulations of the 10 system , both nucleotides switch to  syn conformation of the  dihedral .
the flipping back of the nucleotide into the main body of the dna is observed , but in the wrong position which causes mismatching between the bases , which in turn causes major dna structure deviations as can be seen in both the rmsd and distance plots of supplementary figure 6s .
the cc insertion and base - pair eversion mode of interaction yielded the most stable binding energies ( table 2 ) .
a trend is detected in the binding - energy values that follows the type of interaction each molecule produces in the simulation .
molecules that produced the base - pair flipping have the most stable binding energy as already mentioned that ranges from -44.6 to -38.6 kcal / mol .
the molecules that remained in the minor groove or formed a stacking interaction inside the minor groove with fraying nucleotides ( type a and type d configuration ) presented a binding - energy range from -34 to -27 kcal / mol . molecules that intercalated in the terminal base pairs have a range of binding energy from -25 to -19 kcal / mol .
the remaining molecules that stacked in the terminal base pairs are grouped together ranging from -21 to -15 kcal / mol .
the potential of mean force calculation for system 02 and system 10 is shown in figure 7 .
a 20  distance between the ligand and the dna is where both are separated and is considered as the reference value . in both the systems studied , the copper compound entered the intercalation cavity from the minor - groove side .
the system 02 has two minimums , one at a distance of 2.9  ( with an interaction value of -9.4 kcal / mol ) and the other at 4.5  ( -9.2 kcal / mol ) .
the first minimum is presented where the ligands are inside the intercalation cavity with the non - aromatic ligand toward the major groove and the copper atom interacting with the oxygen atoms dt19:o4 and dg4:o6 .
the second minimum has the ligand shifted toward the minor groove , forming a stacking interaction with the copper complex aromatic ligand and the base pair a6t19 . for system 10
, a single minimum of -10.2 kcal / mol was observed that corresponds to the same conformation obtained from the clustering analysis . as the separation between the dna and the copper compound increases ,
the 02 system shows a smooth linear free energy increase , whereas system 10 shows a decreases of the linear profile caused by the rearrangement of the bulky phenyl groups present in the aromatic ligand .
the distance is between the center of mass of the intercalation cavity and the central atoms of the copper compound ( see  materials and methods section ) .
each point is the average from five independent simulations , error bars show standard deviation between the five copies . what started as a research project to measure binding energies of the casiopenas complexes to dna in the hope of a deeper understanding of the effects of the different substituents on each molecule has opened a new mindset of interaction mode between planar aromatic complexes similar to the ones studied and dna .
it is important to remember that even most of the systems yielded a binding mode type a in the terminal base pairs , this was slightly favored by the complexes with the glycinate systems versus the acetylacetonate systems ( 14 versus 9 ) .
none of the compounds with acac was found in the minor groove whereas three systems were found with the glycinate ; none of the glycinate systems produced base - pair flipping whereas four of the acac systems did ; also for acac , three systems generated terminal base - pair intercalation .
base - pair eversion caused by systems 03 and 05 , which initiated from the aromatic ligands , suggests that molecules with similar aromatic structures that are considered classic intercalators can also interact in this manner , for example , the common dye ethidium bromide which contains a phenanthridine moiety which is similar to the phenanthroline ligand used in this study .
the highlights of our findings include : 
 the possibility to study the intercalation pathway of a tertiary coordination complex with dna starting from the minor groove.the intercalation was due to the loss of the wc pairing between a base pair and the flipping event toward the major groove.the complex moves inside the cavity resulting from the flipping of the nucleobases.intercalation was observed at the terminal base pairs . the possibility to study the intercalation pathway of a tertiary coordination complex with dna starting from the minor groove .
the intercalation was due to the loss of the wc pairing between a base pair and the flipping event toward the major groove .
the compounds that have methyl groups in the aromatic rings ( mono , di - methyl or tetra - methyl ) with either acac or gly do not show a lasting interaction with the dna .
they suggest that the presence of methyl groups or substituents in the positions 4 and 7 of the phenanthroline ring makes more difficult the process of intercalation due to steric hindrance ( 85 ) .
compounds 10 and 19 that have the bulky phenyl groups show a stabilizing interaction with the minor groove , although , the non - aromatic ligands also have an important role . compound 10 that has the acac ligand is interacting with the dna ( with a 93% population ) .
when the non - aromatic ligand is switched to glycinate , the carboxylate side has to be exposed to the solvent to achieve a stable interaction ( with a population of 95% ) .
if the carboxylate side of the glycinate ligand is toward the minor groove , the complex is destabilized ( population of 44% for 19 ) .
we also hypothesize that with longer dna chains , as it is actually found in the cell , the type a interaction will be considerably lower .
longer dna chains will allow the ligands to sample through the grooves for more time before reaching the terminal base pairs and form the stacking interaction .
this extended time could allow the ligand to find an energetically favored binding configuration , which in turn can promote groove binding and possibly also generate the base - pair eversion .
it is clear that extended md simulations , as presented in this study , explored more sampling space than earlier work , provided detailed insight into the spontaneous intercalation processes and created a new set of possible interaction modes to consider when binding molecules to dna .
the rsmd with respect to the experimental structure of the free dna is presented in supplementary figure 1s of the supporting material .
the rmsd values do not include base pairs at the edges of the dna to avoid higher rmsd values caused by terminal base - pair fraying events which has been observed in nucleic acids simulations using the amber force fields ( 47,48 ) . in figure 3 , we show five representations of the interactions found in the 35 simulations studied in this work obtained from the representative structure from the most populated cluster using clustering analysis described in the methodology section . each final configuration ( from a  e , as defined in figure 3 and table 2 ) aids with the analysis of the observed results .
the most populated cluster represents the binding conformation that is populated for the majority of time during the 10 s of md simulation .
the representative binding conformations found for each of the simulations are described in table 2 .
the bottom panel highlights close - up detail of the structure depicted in the top panel .
a : stacking on the terminal base pairs of the dna chain , b : minor - groove binding , c : intercalation with base - pair displacement / opening , d : minor - groove binding with stacking of one of the terminal frayed bases , and e : intercalation near the end of the dna chain .
structures shown are a single - structure snapshot at 10 s simulation time from complexes 05 , 12b , 10 , 09 and 01 .
the percentage denotes the amount of time spent in that particular conformation with respect to the total simulation time .
rmsd values of the highest populated structure from the clustering analysis versus the experimental structure ( pdb entry 1bna ) using all residues .
g values obtained from 5000 frames using the most populated cluster trajectory . from the full data set , 23 of the systems generated an interaction of type a were the copper compound detaches from the minor groove and diffuses along the backbone of the dna until reaching one of the chain edges where the compound stacks on the terminal base pairs .
the dna rmsd values do not increase above 0.2  with respect to the starting structure suggesting that the binding interaction does not produce any significant change in the dna chain for this type of interaction ( table 2 ) .
all of these systems explored multiple conformations and detached completely from the dna , although the prevalent structure as found by the clustering analysis was the stacking interaction at the end of the chain .
as shown in table 2 , the type a binding mode is not as energetically favored as the others .
the type b interactions tend to be more favorable than the terminal base - pair - stacked binding mode .
the type b configuration is characterized by binding of the ligands within the minor groove of the dna duplex with configurations very similar to the starting point for all the md simulations .
the compound 19 remained with both the aromatic and non - aromatic ligand facing the groove while the compound 19 shifted the aromatic ligand toward the groove , exposing the glycinate toward the solvent .
after 125 ns , the complex turned 185 on the z axis until the aromatic moiety and the carboxylic side of the glycinate again faced the inside of minor groove and remained in that configuration for the rest of the simulation ( refer to the structures in the supporting material ) .
the rmsd plot for molecule 03 shows a difference of 1.5  starting after 400 ns and lasted for the rest of the simulation .
the interaction of this molecule with dna resulted in a base - pair eversion / opening toward the major groove followed by intercalation of the copper complex in the resulting pocket ( configuration type c ) .
this interaction , which tends to display the most favorable binding affinity , will be discussed in the next section of the article .
type d configuration consists of the ligand bound to the minor groove but close to one of the ends of the chain where the compound forms a stacking interaction with one of the frayed terminal nucleotides ( residues 1 , 24 , 12 or 13 in the ddd ) . for 09 ,
the compound diffuses inside the minor groove throughout the simulation until it reaches the g10c15 position . in this site ,
the interaction lives for 80 ns and then the g12 base flips back to pair again with c13 .
compound 01 presents low rmsd fluctuations with respect to the reference frame , although the clustering analysis showed an intercalating structure near the terminal base pairs ( configuration type e ) .
the top plot presents the distance measures between the complex and residues 12 and 13 at the edge of the dna chain . at the start of the simulation ,
the compound is located in the minor groove , close to residue 6 , which gives a distance of 2224  . after only 80 ns
, the compound detaches from the minor groove and diffuses close to the dna backbone until it finds the edge of the dna and stacks above residues g12 and c13 ( similar to type a configuration in figure 3 ) .
after 1800 ns of a type a configuration , base opening ( fraying ) of residue g12 toward the minor groove causes an increase in the distance between residue and complex ( marked as b in figure 4 ) .
this causes the complex to reposition toward residue c13 , increasing the stacking interaction and lowering the distance between both molecules .
after 100 ns , the complex slides away from residue c13 and stacks with the base pairs c11 and g14 .
residue 13 flips back to stack again with the complex , but now from the opposite side ( marked as c ) .
residue 13 reforms the watson  crick ( wc ) pairing with residue 12 and the intercalation is complete .
the spike marked as d in figure 4 corresponds to a fraying event of both residues 12 and 13 that lasts for 95 ns before both residues move back into position .
the average structure for the most populated cluster is shown in supplementary figure s2 of si for each ligand molecule . for compound 11 ,
the copper complex detached from the minor groove and shifted toward the end of the dna interacting with the fraying base pairs for 3.1 s until the intercalation formed which then remained for the rest of the simulation . the process of intercalation at the terminal base pairs as observed in md simulation .
top : running average ( 500 frames ) of the distance between g12 and copper complex 01 ( black ) and the distance between c13 and the same compound 01 ( red ) .
the effect of the copper complexes and the five different types of binding to dna observed in the simulations can be summarized comparing the structure from the clustering analysis with the crystal experimental structure of the free dna ( table 2 ) .
as already observed in the plots of rmsd versus time available in the supporting material , the systems that have the highest deviation from the experimental structure are the ones in type c configuration ( 3.34.3  ) , where 2 nucleotides are flipped toward the mayor groove .
the type d configuration also presents high divergence in rmsd due to the distortion of the nucleotides present at the end of the dna chain , flipping toward the minor groove .
the rest of the simulations have values between 1.7 and 2.7  , which are expected values , consistent with previous md simulation work , since we are including the terminal base pairs at each end of the dna that have multiple terminal base - pair fraying events .
this is also expected because the copper complex detaches from the dna and forms a stacking interaction with a base pair at the end of the dna , this actually lowers the frequency of fraying , stabilizing the dna and hence lowering the rmsd difference from the crystal structure .
overall , from the 35 simulations performed , the majority ( 66% ) are type a with similar populations observed for the type b ( 9% ) , type c ( 11% ) , type d ( 6% ) and type e ( 9% ) interactions .
since the simulations were performed without any type of biasing or restraint on the dna or the ligand , the interaction where the copper complex is stacking at the end of the dna is the most common one .
this does not mean that the terminal stacking interaction is energetically favored , in fact , the calculated energy - binding values presented in this work ( table 2 ) show that type a binding has overall lower values compared to the other types of binding . the process for type a of interaction has been previously discussed and observed in similar unbiased , unrestrained drug - dna simulations ( 76,77 ) .
we can attribute this interaction as being predominant in our results mainly to three observations : ( i ) a lack of proper binding of the ligand with the minor - groove environment that causes the ligand to shift through the cavity until there is no more dna ; ( ii ) complete detachment of the ligand from the minor groove toward the solvent ; and ( iii ) steric hindrance caused by the substituent present in the aromatic ligand that causes the whole complex molecule to either shift inside the groove or loose the binding completely .
the majority of the complexes that interact in the type a configuration do not have any substituents in the 4 , 4 position ( in the case of the bipyridine aromatic ligand ) nor in the 4 , 7 position ( for the phenanthroline ligand ) .
none of the complexes with the acetylacetonate ligand ( complexes 01 to 12 ) were found in the minor groove ( type b ) , whereas three of the complexes with the glycinate ligand formed stable minor - groove binding : system 17 that has 5,6-dimethyl substituent and both systems using complex 19 which have 4,7-diphenil substituents in the aromatic ligand .
the 17 system remains bound for the whole simulation , fluctuating between the a5t20 and a6t19 base pairs and forming hydrogen - bonding interactions of the methyl group in position 5 of the phenanthroline ligand with the backbone residues and hydrogen bonding of the aromatic ligand with n3 of adenine .
the presence of the extra phenyl ring in compound 19 stabilizes the binding as the ring goes inside the minor groove , supplying the complex with extra binding sites and hydrogen - bond interactions with the dna backbone .
the system 19 , which started with the carboxylate group of the glycinate toward the minor groove made a 180 turn in the z - axis from 120 to 200 ns time so that the carboxylate group now is facing toward the solvent and the amino part toward the minor groove in the same orientation as the 19 system .
the 19 system fluctuated between t7a18 and t8a17 and only minor readjustments of the ligand within the groove were observed .
insertion and base - pair eversion have been observed experimentally by the barton group ( 7881 ) .
cytosine base pair in the central position and both residues are ejected toward the major groove as the bulky metal complex inserts in the resulting cavity .
, they discovered multiple conformational changes in a 9-mer dna , including the eversion of two guanine residues toward the major groove .
they observed after 300 ns of md simulations the same modes of binding as reported in this work ( end binding , minor - groove binding and intercalation , figure 3 ) .
this is due to the shape of the doxorubicin ligand , which aligned with the top- and bottom - matched base pairs in a perpendicular orientation , allowing one of the flipped bases to flip back into position . in our simulations , compounds 02 , 03
, 05 and 10 started in the minor groove and caused an at pairing to flip or base pair open toward the mayor groove .
the molecules 02 , 03 and 10 pushed open the a5t20 pair and complex 05 the t8a17 pair .
the copper compound moved inside the resulting pocket ( shown in figure 5 ) , which corresponds to the average structure of the highest populated cluster for each simulation using the full 10 s trajectory . in figure 5 , we present the first 500 ns of simulation for each of the simulations that generated the flipped nucleotides .
it is interesting to observe that the intercalation process starts early in the simulation ; the 05 and 10 systems start almost at similar time ( 160 ns ) , followed by 02 ( 308 ns ) and finally the 03 system ( 360 ns ) , and the intercalation process takes no longer than 20 ns except for the 10 system which lasts for more than 100 ns .
this system takes longer time to reach the final intercalated position due to the rearrangement of the two bulky phenyl groups in the minor groove .
top : average structure extracted from the most populated cluster using all the simulation frames for systems 02 , 03 , 05a and 10 .
bottom : running average ( 200 frames window ) of the rmsd values of the systems in the configuration type c. rmsd values calculated using the crystal structure as reference and central residues 2 to 11 and 14 to 23 .
the process of intercalation follows the general sequence : ( 1 ) ligand binding in the minor groove , ( 2 ) ligand positioning , ( 3 ) partial opening of base pairs and partial insertion of the ligand , ( 4 ) flipping out of both the bases and full ligand insertion in the cavity , followed by ( 5 ) reorientation of the ligand inside the cavity .
all four systems that were found to produce this type of intercalation follow this same basic process or model , with variation only on how much time a certain system spent in each of the 15 steps .
the main observed difference is with step 3 where the copper complex partially pushes the base pair until the wc hydrogen bonds break and the compound pushes the nucleobases toward the major groove .
this step happens almost instantaneously in systems 03 and 05 , where the at base pair is pushed with atoms 2 , 3 and 4 of the phenanthroline ring .
regardless of this observation , the process of intercalation is similar in each of the four systems . in previous dna  ligand intercalation studies , the rise , roll and twist helical parameters provide a good description of the intercalation process ( 83,84 ) . in these simulations where a base pair is structurally separated and pushed toward the minor groove , there is a coordinated global change in the flipping nucleotides , affecting both inter- and intra - base - pair geometries .
structural parameters of the flipped nucleotide for each case mentioned above and the distance of the copper atom of the ligand ( used as the center of the molecule ) to the o4 and o2 atoms of the thymine base in the intercalation pocket are presented in supplementary figure 3s of the supporting material .
the first step of the intercalation process was the presence of the bound ligand in the minor groove , which in this set of simulations was also the starting position .
the distance between the ligand and the o4 and o2 atoms used as a set of reference points in the intercalation cavity is > 9  in each case ( supplementary figure s3 ) . in step 2 , the ligand reorients multiple times in the minor groove ( a and d in figure 6 ) . in figure 5
, we can see the amount of time the ligand explored the space inside the minor groove until moving to step 3 . until this step ,
the structural parameters in supplementary figure 3s fluctuate close to the canonical values . during step 3 , the opening angle and the stretch distance increase , in addition to the propeller angle , and
this indicates that one of the bases is being pushed ( b and e in figure 6 ) .
step 4 is the flip of both nucleotides toward the major groove and partial insertion of the ligand inside the resulting cavity . in this step ,
the rise of the previous and next base pairs that form the intercalation pocket goes down , which puts the four nucleotides closer together hence reducing the volume of the intercalation pocket .
step 5 presents the deep insertion of the ligand inside the pocket ( c and f in figure 6 ) .
once inside , the ligand shows rotation on the z - axis , reorienting itself multiple times during the simulation . in supplementary figure s4
the orientation of the ligand inside the cavity formed a 83 angle with respect to the bases shifting 11% of the total simulation time .
system 05 shifted continuously between angles of 200 and 90. finally , system 10 shifted conformation upon entering the cavity , forming a 290 and shifting conformation in a small interval of 20 for the rest of the simulation .
this stability inside the pocket for this compound is due because of the steric hindrance caused by both phenyl groups in the aromatic ligand of the complex .
orange and yellow are the top and bottom base pairs forming the intercalation cavity after base flipping .
if the intercalation  base - flipping events are truly an effect of the binding of the casiopenas complexes , simulations were performed using the ligands with no metal center , thus , only the two aromatic ligands and the two non - aromatic ligands : acetylacetonate , glycinate , bipyridine and phenanthroline .
each molecule was parameterized using the same methodology described and manually placed in the same position inside the minor groove of the dna 12-mer in the same manner as the casiopenas compounds .
each of the four systems was run with three independent copies using the same equilibration and production protocol and randomizing the position of the counter ions for each copy .
supplementary figure 5s in the supporting material shows the rmsd values for each of the copies .
visual inspection of the trajectories showed transient and short - lived ( < 50 ns ) interactions with the terminal base pairs of the dna with the aromatic ligands .
the non - aromatic ligands detach from the minor groove early in the simulation ( 25 ns ) and only short - lived interactions with the dna 's backbone are detected . considering the prolonged sampling time used ,
these results indicate that intercalation between two base pairs and intercalation with the opening of internal base pairs are unlikely to be detected using the ligands by themselves .
additional analysis was done to understand the base - pair eversion event of systems 03 and 10 where the copper complex was removed from the opening of the everted nucleotides to examine if a canonical b - dna structure could be reformed .
the starting structure was obtained using only the trajectory frames were the complex was inserted in the dna and clustered the frames using the same methodology as described previously .
a representative structure of the top populated cluster was then extracted , the copper complex was deleted and three independent simulations were run , using the same equilibration and production protocol , randomizing the counter ion positions for each copy .
supplementary figure 6s in the supporting material shows the rmsd value of the dna and distance from the center of mass for each flipped nucleotide of each copy .
system 03 reforms the wc pairing in residues a5t21 in 1820 ns of simulation time as can be observed in the distance top plot in supplementary figure 6s . in two of the simulations of the 10 system , both nucleotides switch to  syn conformation of the  dihedral .
the flipping back of the nucleotide into the main body of the dna is observed , but in the wrong position which causes mismatching between the bases , which in turn causes major dna structure deviations as can be seen in both the rmsd and distance plots of supplementary figure 6s .
the cc insertion and base - pair eversion mode of interaction yielded the most stable binding energies ( table 2 ) .
a trend is detected in the binding - energy values that follows the type of interaction each molecule produces in the simulation .
molecules that produced the base - pair flipping have the most stable binding energy as already mentioned that ranges from -44.6 to -38.6 kcal / mol . the molecules that remained in the minor groove or formed a stacking interaction inside the minor groove with fraying nucleotides ( type a and type d configuration ) presented a binding - energy range from -34 to -27 kcal / mol .
molecules that intercalated in the terminal base pairs have a range of binding energy from -25 to -19 kcal / mol .
the remaining molecules that stacked in the terminal base pairs are grouped together ranging from -21 to -15 kcal / mol .
the potential of mean force calculation for system 02 and system 10 is shown in figure 7 . a 20  distance between the ligand and the dna
is where both are separated and is considered as the reference value . in both the systems studied , the copper compound entered the intercalation cavity from the minor - groove side .
the system 02 has two minimums , one at a distance of 2.9  ( with an interaction value of -9.4 kcal / mol ) and the other at 4.5  ( -9.2 kcal / mol ) .
the first minimum is presented where the ligands are inside the intercalation cavity with the non - aromatic ligand toward the major groove and the copper atom interacting with the oxygen atoms dt19:o4 and dg4:o6 .
the second minimum has the ligand shifted toward the minor groove , forming a stacking interaction with the copper complex aromatic ligand and the base pair a6t19 . for system 10 ,
a single minimum of -10.2 kcal / mol was observed that corresponds to the same conformation obtained from the clustering analysis .
as the separation between the dna and the copper compound increases , the 02 system shows a smooth linear free energy increase , whereas system 10 shows a decreases of the linear profile caused by the rearrangement of the bulky phenyl groups present in the aromatic ligand .
the distance is between the center of mass of the intercalation cavity and the central atoms of the copper compound ( see  materials and methods section ) .
each point is the average from five independent simulations , error bars show standard deviation between the five copies . what started as a research project to measure binding energies of the casiopenas complexes to dna in the hope of a deeper understanding of the effects of the different substituents on each molecule has opened a new mindset of interaction mode between planar aromatic complexes similar to the ones studied and dna .
it is important to remember that even most of the systems yielded a binding mode type a in the terminal base pairs , this was slightly favored by the complexes with the glycinate systems versus the acetylacetonate systems ( 14 versus 9 ) .
none of the compounds with acac was found in the minor groove whereas three systems were found with the glycinate ; none of the glycinate systems produced base - pair flipping whereas four of the acac systems did ; also for acac , three systems generated terminal base - pair intercalation .
base - pair eversion caused by systems 03 and 05 , which initiated from the aromatic ligands , suggests that molecules with similar aromatic structures that are considered classic intercalators can also interact in this manner , for example , the common dye ethidium bromide which contains a phenanthridine moiety which is similar to the phenanthroline ligand used in this study .
the highlights of our findings include : 
 the possibility to study the intercalation pathway of a tertiary coordination complex with dna starting from the minor groove.the intercalation was due to the loss of the wc pairing between a base pair and the flipping event toward the major groove.the complex moves inside the cavity resulting from the flipping of the nucleobases.intercalation was observed at the terminal base pairs . the possibility to study the intercalation pathway of a tertiary coordination complex with dna starting from the minor groove .
the intercalation was due to the loss of the wc pairing between a base pair and the flipping event toward the major groove .
the compounds that have methyl groups in the aromatic rings ( mono , di - methyl or tetra - methyl ) with either acac or gly do not show a lasting interaction with the dna .
they suggest that the presence of methyl groups or substituents in the positions 4 and 7 of the phenanthroline ring makes more difficult the process of intercalation due to steric hindrance ( 85 ) .
compounds 10 and 19 that have the bulky phenyl groups show a stabilizing interaction with the minor groove , although , the non - aromatic ligands also have an important role . compound 10 that has the acac ligand is interacting with the dna ( with a 93% population ) .
when the non - aromatic ligand is switched to glycinate , the carboxylate side has to be exposed to the solvent to achieve a stable interaction ( with a population of 95% ) . if the carboxylate side of the glycinate ligand is toward the minor groove , the complex is destabilized ( population of 44% for 19 ) .
we also hypothesize that with longer dna chains , as it is actually found in the cell , the type a interaction will be considerably lower .
longer dna chains will allow the ligands to sample through the grooves for more time before reaching the terminal base pairs and form the stacking interaction . this extended time could allow the ligand to find an energetically favored binding configuration , which in turn can promote groove binding and possibly also generate the base - pair eversion .
it is clear that extended md simulations , as presented in this study , explored more sampling space than earlier work , provided detailed insight into the spontaneous intercalation processes and created a new set of possible interaction modes to consider when binding molecules to dna .
 the family of anticancer complexes that include the transition metal copper known as casiopenas shows promising results . 
 two of these complexes are currently in clinical trials . 
 the interaction of these compounds with dna has been observed experimentally and several hypotheses regarding the mechanism of action have been developed , and these include the generation of reactive oxygen species , phosphate hydrolysis and/or base - pair intercalation . to advance in the understanding on how these ligands interact with dna 
 , we present a molecular dynamics study of 21 casiopenas with a dna dodecamer using 10 s of simulation time for each compound . 
 all the complexes were manually inserted into the minor groove as the starting point of the simulations . 
 the binding energy of each complex and the observed representative type of interaction between the ligand and the dna is reported . with this extended sampling time 
 , we found that four of the compounds spontaneously flipped open a base pair and moved inside the resulting cavity and four compounds formed stacking interactions with the terminal base pairs . 
 the complexes that formed the intercalation pocket led to more stable interactions . 
INTRODUCTION
MATERIALS AND METHODS
RESULTS AND DISCUSSION
Interaction of the copper complexes with DNA
Base-pair eversion and intercalative interactions
Energy analysis between Casiopenas and DNA
SUPPLEMENTARY DATA