LIGAND-INDUCED DNA CONFORMATIONAL CHANGES IN PROFLAVINE MINOR GROOVE-BOUND COMPLEXES STUDIED BY MOLECULAR DYNAMICS SIMULATION

Background: Minor groove binding is a rate-limiting step in proflavine-DNA intercalation reaction. This step is believed also to be responsible for the sequence-dependent kinetics of proflavine binding to DNA. At the same time, most studies are focused on the final stage of the reaction – the intercalation complex, and there is a lack of data concerning the structure and stability of proflavine-DNA minor groove-bound complexes. Objectives: The objective of this study was to investigate the stability of proflavine minor groove-bound complexes with DNA oligonucleotides of different sequence by molecular dynamics simulation and to analyze the DNA conformational changes caused by the proflavine binding. Materials and methods: The molecular dynamics simulations of proflavine minor groove-bound complexes with poly(dA)·poly(dT) and poly(dCG)·poly(dCG) oligonucleotides of 30 bp length were done in program package AMBER12 with explicit water (SPC/E) and ions (NaCl 0.15 M) using force fields FF14SB for DNA and GAFF for ligand. The starting configurations of complexes were obtained by docking method in AutoDock 3.05. After multi-stage equilibration protocol, each system was simulated at T=300 K and p=1 bar for a 50 ns production phase. Then trajectories were post-processed in AMBERTools17 and VMD-1.9.3 packages. Results: Our simulations confirm that proflavine-DNA minor groove-bound complexes are stable in the 50 ns time range but there are some structural rearrangements in them with respect to the initial structures. The narrowing of the DNA minor groove is observed in the proflavine binding site. In proflavine-poly(dCG)·poly(dCG) complex it is more pronounced and is accompanied by the BI/BII transitions in DNA and the reorientation of ligand. In proflavine-poly(dA)·poly(dT) complex the specific intermolecular hydrogen bonds are formed, which are optimized by the changes in opening and propeller twisting of involved AT-base pairs. Complexes are stabilized by the van der Waals and hydrophobic interactions, which are more favorable in the proflavine-poly(dA)·poly(dT) complex. Conclusions: Our results show that the binding of proflavine to a minor groove of DNA induces the conformational changes in the DNA that are important for the resulting complex stability.

Proflavine (PF) (Fig. 1) is a small aromatic molecule that is known to exert antibacterial and mutagenic properties due to its interaction with DNA [1]. Depending on drug-DNA concentration, two binding modes are observed for PF: the strong binding of individual PF molecules to DNA at low drug-DNA ratios and weak binding of PF aggregates to DNA at high drug-DNA ratios [2]. To explain the strong binding of PF to DNA, an intercalation model was proposed by Lerman [3], according to which the planar drug is inserted between the base pairs into the DNA double helix with drug plane being perpendicular to the helix axis and long axis of the drug chromophore being almost parallel to the long axis of the base pair. This model became generally accepted due to the large number of experimental facts supporting it [3][4][5][6]. Later the X-ray structures of PF complexes with dinucleoside monophosphates [7][8][9] and short DNA duplex [10] were resolved also confirming the Lerman intercalation model. Despite the overall agreement on the structure of PF-DNA intercalation complex, the mechanism of the intercalation reaction has been controversial. Early kinetic studies on strong PF binding to calf-thymus DNA, performed by Li [11] using joule-heating temperature-jump method with absorption detection, revealed two relaxation times in sub-millisecond and millisecond range with different dependence on concentration. The two-step scheme of the reaction was proposed: the strong outside-bound complex is formed between the PF and DNA on the first stage and the intercalation complex is formed on the second stage (Scheme 1). The absorption spectrum of the outside-bound complex was similar to that of the intercalation complex, implying that there was at least a partial overlap of bases with PF in the outsidebound complex.
( ) ( ) 23 The subsequent investigations showed that the kinetics of the PF intercalation reaction depends strongly on the DNA base composition [12,13]. Using joule-heating temperaturejump method with fluorescence detection, Ramstein found that there was only one relaxation time in the relaxation signal of PF-poly[d(GC)] system, and the inverse of the relaxation time varied linearly with concentration [13]. He suggested a single-step reaction for binding of PF to poly[d(GC)] and the resulting complex was interpreted as a strong outside-bound complex (Scheme 2). For the PF-poly[d(AT)] system also only one relaxation time was found, though its inverse leveled off at high phosphate concentration, and the two-step mechanism of the reaction proposed by Li was found to be appropriate (Scheme 3). An additional slower 8 K.V. Miroshnychenko, A.V. Shestopalova relaxation time independent on concentration was seen for the interaction of PF with natural DNAs. Ramstein supposed that it was due to the PF exchange between GC-and AT-sites, having different kinetics of dye binding.
Later, the kinetics of PF binding to DNA was reinvestigated by Marcandalli using iodinelaser temperature-jump technique with absorption [14] and fluorescence detection [15]. One relaxation time was detected regardless of the choice of the method used. The dependence of reciprocal relaxation time on DNA concentration was linear in the case of GC-rich DNA and tended to a plateau at high DNA concentrations for AT-rich DNA. To explain the results Marcandalli suggested that there were two types of sites in DNA: B 1 sites(connected to GCpairs) and B 2 sites (connected to AT-pairs). But unlike Ramstein, he proposed two-step mechanism of the reaction for both sites, and the complex formed at the final stage of the reaction was believed to be of intercalation type in both cases (Scheme 4, 5). The complexes formed on the first stage, E and E', were rapidly formed outside-bound complexes, but in contrast to the outside-bound complex from the Li scheme [11], these complexes were spectrophotometrically similar to free dye molecules, implying that they were water-solvated and there was no direct overlap of PF with DNA bases in them. Marcandalli showed that the magnitude of K e B i (K e -is the equilibrium constant of the first step of the reaction (Scheme 4, 5), B i -number of corresponding sites) determined the dependence of reciprocal relaxation time on DNA concentration. In the case of GC-rich DNA the outside-bound complex of PF was weak resulting in linear concentration dependence. For the AT-rich DNA strong outsidebound complex was formed between PF and DNA, which led to the leveling off of reciprocal relaxation time at high DNA concentrations and slower intercalation compared to the GCsites. Marcandalli assumed that the greater stability of the outside-bound complex of PF with AT-rich DNA could be due to the hydrogen bonding of PF with purine N3 and pyrimidine O2 atoms in the minor groove of DNA. The N2 amino group of guanine protruding into the minor groove would interfere such interactions causing the weaker outside binding of PF with GC-rich DNA.
Nowadays the detailed picture of drug-DNA interactions at the atomic level can be obtained using computer simulation methods. Most of the computational studies of PF-DNA interactions considered the final stage of the reaction -the intercalation complex [16][17][18][19][20][21]. The outside binding of PF was modeled in the pioneering works as a hydrogen bonding to phosphate group [22], based on the 3:2 PF-CpG structure of Berman [8], or as an incomplete intercalation from the minor groove side [23]. To the best of our knowledge, the first modeling study that investigated the outside-bound complex of PF with DNA oligomer was done by Sasikala [24]. Using metadynamics method he studied the dissociation of PF from the minor groove-bound state and with appropriate configurational restraint managed to achieve the intercalation of PF from the minor groove-bound state. In the subsequent work of Sasikala, the intercalation of PF from the major groove-bound state was achieved and the processes of intercalation, de-intercalation and dissociation of PF from both grooves of DNA were compared [25]. It was found that the intercalation and de-intercalation of PF proceeded through the major groove side but the timescale of the intercalation reaction was determined by the stable pre-intercalative minor groove-bound state. According to Sasikala, the binding of PF to DNA can be described by the following scheme: first, drug forms quickly (in a nanosecond range) groove-bound states with the population of the minor groove-bound state being higher due to its greater stability (compared to the major groove-bound state), but the direct intercalation from the minor groove-bound state is extremely slow due to the high energy barrier, therefore drug dissociates from the minor groove-bound state in a millisecond range and forms the intercalation complex from the major groove side (Scheme 6).
Recently the outside-binding of multiple PF molecules to DNA was simulated by Sasikala [26]. The formation of PF aggregates around DNA in nanosecond range was observed. Aggregates were seen mostly in the major groove of DNA, though a few stacks of PF molecules in the minor groove of DNA were also found. Due to the low phosphate to dye ratio (P/D=1) used in that study, the simulated process can be attributed to the weak binding of PF to DNA. It would be interesting though to investigate the strong monomeric minor groove binding of PF in more detail. There is a lack of simulation data on minor groove binding of PF to different DNA sequences. At the same time, there is multiple experimental evidence that the kinetics of PF binding to DNA depends on the DNA composition [12,13,15], and the analysis of literature suggests that it is the minor groove-bound state of PF that is responsible for the timescale of the reaction and the sequence-dependent kinetics [15,25].
Previously we used molecular docking method to obtain the minor groove-bound (MGB) complexes of PF with DNA oligonucleotides of different sequence built according to the Xray fiber data [27]. Depending on the geometry of the minor groove of DNA, two types of PF MGB complexes were observed. For DNA-targets with narrow minor groove (mostly poly(dA)·poly(dT) sequences), the strong MGB complex of PF was formed: its energy was comparable or even more favorable than that of the intercalation complex of PF. In this type of MGB complex, the PF chromophore was laying along the groove and its amino groups were participating in hydrogen bonding with DNA bases. For DNA-targets with wide minor groove (DNA of GC-composition), the second type of PF MGB complex was obtained: the PF chromophore was located across the minor groove and its amino groups formed hydrogen bonds with one of the sugar-phosphate chains. This type of complex was less energetically favorable than the intercalation complex of PF and the first type of PF MGB complex. The aim of this study was to investigate the stability of PF MGB complexes formed with DNA oligonucleotides of different sequence using molecular dynamics (MD) simulation with explicit solvent and physiological salt concentration and to analyze the DNA conformational changes caused by the PF binding.

MATERIALS AND METHODS
The DNA oligonucleotides of 30 base pairs long were constructed in the module NAB of AMBERTools17 [28] or program FIBER of 3DNA package [29]. The protonated form of PF molecule was optimized at the B3LYP/6-31G* level in the Gaussian 03 software [30], and the RESP atomic charges were calculated at the HF/6-31G* level of theory using RED III package [31]. The structures of PF MGB complexes were obtained by molecular docking method in the AutoDock 3.05 package [32] using Lamarckian genetic algorithm with the following parameters: the size of population = 150 individuals, the maximum number of energy evaluations = 10 7 , the maximum number of generations = 27000. The Kollman partial atomic charges were used for DNA-targets and the RESP atomic charges -for PF molecule during docking procedure. All hydrogen atoms in the structures of target and ligand were considered explicitly and treated as non-polar. Under these conditions the NDB structures of PF complexes with dinucleoside monophosphates (NDB codes: ddb009, ddb033, ddb034, drb008, drd004 [33]) were well reproduced in test re-docking calculations. The docking search space was confined by the grid box of 70x70x70 points centered at the fifteenth step of the DNA oligonucleotide with distance between adjacent grid points equal to 0.375Å. 100 dockings were performed for each DNA-target, and the obtained complexes were ranked by the docking energy and clustered using a tolerance of 1 Å. The best cluster containing more than 10 conformations was used as a starting structure for MD simulations.
The MD simulations were performed for MGB complexes of PF with poly(dA)·poly(dT) and poly(dCG)·poly(dCG) oligonucleotides, and separately for DNA oligonucleotides and the PF. The program package AMBER12 [34] was used with the force fields FF14SB [35] for DNA and GAFF [36] for PF. Each system was solvated with SPC/E water molecules [37] in a truncated octahedral box, and periodic boundary conditions were applied. The minimal distance from solute atoms to the edge of the box was greater than 15 Å. The sodium counterions were used to neutralize the DNA charge, and an additional amount of Na + and Clions was added to model the 0.15M concentration of NaCl salt. The Joung parameters for ions were used [38]. The initial positions of ions were randomized so that there were no ions within 6 Å of solute, and the minimal distance between the ions was greater than 4 Å. The long-range electrostatic interactions were treated using particle-mesh Ewald method [39] with a direct space cutoff of 9 Å. The bonds involving hydrogen were fixed using SHAKE algorithm [40], and a 2 fs integration step was applied. To equilibrate the system, the multistage protocol was used. First, the positions of water and ions were minimized keeping the solute fixed with weak restraints (25-20 kcal/(mol·Å 2 ). Then the system was heated during 100 ps from 100 to 300 K under constant volume condition and 15 kcal/(mol·Å 2 ) restraints on solute. After that, the restraints on solute atoms were gradually relaxed from 15 to 1 kcal/(mol·Å 2 ) during a 5 series of alternating stages of 5000 steps minimization and 50 ps NPT MD at T = 300 K and p = 1 bar. To control the temperature and pressure, the Berendsen thermostat and barostat [41] were used with a coupling constant of 0.2 ps for both parameters. The next step was the 50 ps NPT MD with 0.5 kcal/(mol·Å 2 ) restraints on solute followed by the 100 ps unrestrained NPT MD. The last stage of the equilibration protocol was the unrestrained 1 ns NPT MD at T = 300K and p = 1bar with the coupling constants increased to 5 ps. Then the simulations were extended under the same conditions for a 50 ns production phase.
The MD trajectories were analyzed using CPPTRAJ module [42] of AmberTools17 [28] and visualized in VMD-1.9.3 [43]. The hydrogen bonds (HB) were determined using the geometric criteria: the distance between the HB donor and the HB acceptor should be less then 3.2 Å, and the angle formed by the HB donor, the hydrogen and the HB acceptor should be greater than 120°. The ion was considered to be located in the minor groove of a DNA helical step when it was within 3 Å of N3 purine or O2 pyrimidine atoms of bases of this helical step.
The binding free energy estimates were obtained according to the molecular mechanics Poisson-Boltzmann surface area methodology (MMPBSA) using MMPBSA.py script [44] from AmberTools17. Under this approach, the binding free energy is calculated as a difference between the free energy of complex and the sum of the free energies of the unbound receptor and ligand: Each term in this equation can be represented as a sum of three components: the gas energy, calculated from the force field, the solvation energy and the entropy contribution: Ligand-induced DNA conformational changes in proflavine minor groove-bound… 11 The solvation free energy can be divided into the polar and non-polar part: (9) The polar part of the solvation free energy was calculated by solving nonlinear Poisson-Boltzmann equation. The ionic strength was set to 0.15 M. The dielectric constants of solute and solvent were equal to 4 and 80, correspondingly. The ratio between the longest dimension of the rectangular finite-difference grid and that of the solute was equal to 4.0. The grid spacing of 0.5 Å was used.
The non-polar part of the solvation free energy was estimated as a sum of the repulsive term related to the formation of the solute-sized cavity in the solvent and the attractive dispersion term describing the attractive solute-solvent interactions. The repulsive term was calculated as where SAV -is the molecular volume enclosed by the solvent-accessible surface, p=0.0378 kcal/(mol·Å 3 ) -is the solvent pressure parameter, and c=-0.5692 kcal/mol. The dispersion term was computed with a surface-based integration method [45].
Each energy term was calculated as an average over the ensemble of representative structures taken from the trajectory. Both single and multiple trajectory approaches were used and compared. We used for the calculation every 10 th frame of the trajectory production phase (at every 20 ps, 2500 frames total).
The entropy contribution was calculated as a sum of translational, rotational and vibrational terms using standard statistical mechanics formulas for the rigid rotor harmonic oscillator ideal gas approximation [46]. Due to the computational cost, the entropy calculation was done only for a subset of 50 trajectory frames taken at every 1 ns. These frames were minimized using implicit solvent and then subjected to the normal mode analysis. It took more that 40 hours on 8-core node (2x Intel® Xeon® CPU 5420) to complete the minimization and normal mode calculation for 50 frames.

Average MD structures and root-mean-square deviation of PF from the initial position
The starting structures of PF MGB complexes with poly(dA)·poly(dT) and poly(dCG)·poly(dCG) oligonucleotides obtained by docking method are shown in Fig. 2 (a) and (b). As one can see, the position of PF in the DNA minor groove is different in these two complexes. The poly(dA)·poly(dT) oligonucleotide has narrow minor groove, and in complex with it, the PF chromophore is located along the groove with its protonated ring nitrogen and amino groups being directed into the groove and making hydrogen bonds with DNA. The poly(dCG)·poly(dCG) oligonucleotide has wider minor groove compared to the poly(dA)·poly(dT), and in this case, the PF chromophore is positioned across the groove with its protonated nitrogen interacting with one of the sugar-phosphate chains.
To estimate the stability of PF MGB complexes in MD simulation, the root-mean-square deviation (RMSD) of PF heavy atoms with respect to the starting structures was calculated (the structures were aligned by the 4 base pairs closest to the bound PF) (Fig. 3).From the RMSD values it follows that the PF MGB complexes are stable in the 50 ns time interval. The average RMSD value of PF in complex with poly(dA)·poly(dT) is quite small (1.38±0.28 Å) indicating that the structure of the docking complex is well conserved during the MD trajectory. Indeed, the average MD structure of PF MGB complex with poly(dA)·poly(dT) is very much alike to the starting structure ( Fig. 2 (a) and (c)). On the contrary, the average RMSD value of PF in complex with poly(dCG)·poly(dCG) is relatively high and equals to 4.09±0.93 Å. This means that though the ligand is in close contact with poly(dCG)·poly(dCG) during all the trajectory time there are some structural rearrangements in this complex with respect to the starting structure. Visual inspection of the trajectory reveals that during the MD simulation the reorientation of ligand takes place that causes the narrowing of the poly(dCG)·poly(dCG) minor groove in the PF binding site. The chromophore of ligand becomes located along the groove, like in the complex of PF with poly(dA)·poly(dT), except that amino groups of PF are directed outside the groove (Fig. 2 (d)). As it is evident from the time dependence of the RMSD, the reorientation of ligand occurs already at the equilibration phase (Fig. S1).The different orientation of PF amino groups in the average MD structures of PF MGB complexes can be explained by the fact that poly(dCG)·poly(dCG) has more shallow minor groove compared to the poly(dA)·poly(dT) due to the NH 2 -groups of guanine protruding into the minor groove.

Hydrogen bonds
The orientation of PF amino groups affected the formation of intermolecular HB in PF MGB complexes. There were no intermolecular HB in the PF-poly(dCG)·poly(dCG) complex during the most part of the simulation time. And in the MD trajectory of PFpoly(dA)·poly(dT) complex, the average number of HB between the PF and DNA was 1.77. In Fig. 4, the snapshots from the MD trajectory of the PF-poly(dA)·poly(dT) complex with depicted intermolecular HB are given. As shown, there are two PF groups, which make HB with DNA (thymine O2 or adenine N3 atoms) at the same time: the protonated ring nitrogen and one of the PF amino groups. This observation is in agreement with the prediction of Marcandalli that strong outer complex with AT-rich DNA is formed only for those acridine derivatives that have HB donor groups simultaneously in the positions 6 and 10 of the acridine chromophore [15].

Interactions with solvent and ions
In both PF MGB complexes, water-mediated intermolecular HB (water bridges) were seen during the MD simulation. The average number of water bridges between the PF and DNA was equal to 1.2 for the PF-poly(dA)·poly(dT) complex and 1.6 for the PFpoly(dCG)·poly(dCG) complex. The snapshots of the PF MGB complexes with the maximum observed number of water bridges are shown in Fig. S2. In the PF-poly(dCG)·poly(dCG) complex, all the PF nitrogens were solvated and could interact via water bridges with DNA. And in the PF-poly(dA)·poly(dT) complex, only amino-groups of PF were solvated whereas the protonated ring nitrogen was unavailable to water because it was deeply inserted into the DNA minor groove and most part of the time participated in direct hydrogen bonding with DNA. The lifetime of the HB between PF and water molecules was in range 10-20 ps. Due to the PF-DNA interactions, the total number of water molecules bound via HB to each MGB complex was ~5 molecules less than the sum of the waters bound via HB to its free components (Table S1). The intrusion of the PF into the DNA minor groove caused not only the change of the hydration pattern but also the redistribution of the sodium ions in the DNA minor groove (Fig. 5). More noticeable changes were seen for the poly(dA)·poly(dT), which had more Na + ions in its minor groove during the simulation compared to the poly(dCG)·poly(dCG) in line with other simulation studies [47,48]. The residence times of Na + ions within the minor groove of the PF MGB complexes were about one-and-half to two times shorter than that for the free DNA oligonucleotides.

Helical parameters
To study the structural deformations of DNA caused by the PF minor groove binding, the helical parameters of DNA in complexes were calculated for every 40 ps of the MD trajectory (using CPPTRAJ) and compared to the parameters obtained from the MD trajectories of free DNA oligonucleotides. The most obvious change was the narrowing of the DNA minor groove in the PF binding site and the neighboring helical steps with respect to the free DNA oligonucleotides (Fig. 6). The effect was more pronounced for the PF-poly(dCG)·poly(dCG) MGB complex: the groove narrowing in the PF binding site was up to 2 Å compared to the free DNA oligonucleotide. In complex with poly(dA)·poly(dT), the narrowing of the minor groove caused by PF was circa 1 Å.
A partial widening of the minor groove observed at the initial steps for the poly(dA)·poly(dT) oligonucleotide both in free state and in complex with PF ( Fig. 6) was due to the fraying of the terminal AT base pair and the subsequent formation of the stacked noncanonical structure between the terminal thymine and adenine during the MD simulation. The fraying of the terminal base pairs in MD simulations is a known force field artifact [49]. However, we expect that it slightly affects the final geometry of the PF binding site because PF was bound in the central part of the 30bp DNA oligonucleotide.
The narrowing of the minor groove was accompanied by the variations of other DNA helical parameters. The increased twist (overwinding) was found in the PF binding site, which was partially compensated by unwinding in the adjacent steps (Fig. S3). The changes of other helical parameters involved decrease of roll with change of sign from positive to negative, more positive x-displacement and more positive slide in complex with poly(dCG)·poly(dCG) (Fig. S4-S6). Those changes are characteristic for BI/BII transitions of DNA backbone that are caused by coupled changes of dihedral angles ε and ς from (t,g-) in BI state to (g-,t) in BII state [50,51]. The difference ε-ς is around -90º for BI state and around 90º for BII state. The analysis of DNA backbone angles in the PF binding site showed that, in the MD trajectory of PF-poly(dCG)·poly(dCG) complex, two nucleotides (C15 and G46) had significant population of BII state (Fig. 7). The transitions to the BII state caused a shift in the distribution of sugar puckers of cytosines (C13, C15, C45, C47) in the PF binding site to the south region (Fig.  S7). The adaptation of the double helix structure to the backbone transitions led to the enhanced fluctuations of propeller twist near the bound PF in the PF-poly(dCG)·poly(dCG) complex ( Fig. S8 (a)). The BI/BII transitions are often seen upon binding of proteins to DNA and are believed to play an important role in the indirect DNA recognition [52][53][54]. Unlike the PF-poly(dCG)·poly(dCG) complex, the nucleotides in the PF binding site in its complex with poly(dA)·poly(dT) were mostly in the BI state during the MD simulation. This difference can be rationalized by two reasons. First, the poly(dA)·poly(dT) already has a narrow minor groove, and small conformational changes are needed to provide tighter proflavine binding. Second, the BI/BII conformational transitions are sequence specific: the CG and GC steps are known to undergo BI/BII transition more easily than AA or TT steps [55,56]. The base pairs of the PF binding site in the PF-poly(dA)·poly(dT) complex were characterized by more negative propeller twist and an increased opening (Fig. S8 (b) and S9  (b)). These changes most likely were due to the hydrogen bonding of PF with N3 atoms of adenine (A17, A16) and/or O2 atoms of thymine (T44, T45) (Fig. 4). The sugar pucker distributions of A18 and A17 nucleotides were shifted to the north region (Fig. S10) possibly as a result of the interaction of their O4' atoms with the aromatic system of the PF chromophore (Fig. S11).

MMPBSA free energy estimates
To estimate quantitatively the stability of PF MGB complexes, we calculated the PF binding free energies using MMPBSA approach. The results are given in Table 1. As one can see, the MGB complex of PF with poly(dA)·poly(dT) is enthalpically more favorable than that with the poly(dCG)·poly(dCG). The stabilization of complexes was due to the van der Waals interactions and hydrophobic interactions (there is no need to form an extra cavity in solvent because PF fits well into the DNA minor groove) (Table S2). At the same time, the intrusion of PF disturbed the DNA-solvent interactions in the minor groove that caused the unfavorable ΔG disp term (Table S2).
The entropy penalty was less in the case of PF-poly(dCG)·poly(dCG) complex. Analysis of the entropy components showed that the major unfavorable contribution (~22 kcal/mol) was due to the loss of translational and rotational degrees of freedom, which was partially compensated by the favorable vibrational term (Table S3). It should be noted that the entropy gain due to the release of tightly bound water molecules upon PF binding from the spine of hydration of poly(dA)·poly(dT) or ribbon of hydration in the case of poly(dCG)·poly(dCG) was not taken into account. Therefore we suppose that the obtained entropy penalties are overestimated.
Values of the energy terms calculated using multiple trajectory approach showed high fluctuations. But this approach allowed to access the cost of the DNA conformational changes upon the PF binding. For the PF-poly(dCG)·poly(dCG) complex, the unfavorable conformational term (ΔG conf ) was obtained that is consistent with the observed in the MD trajectory notable conformational changes of the poly(dCG)·poly(dCG) minor groove in the PF binding site (Table S2). All energy terms are given in kcal/mol. ST and MT are the single and multiple trajectory approaches, correspondingly. ΔH is the sum of all energy terms except entropy contribution; -TΔS is the entropy penalty due to the loss of translational, rotational and vibrational degrees of freedom upon ligand binding; ΔG Total = ΔH-TΔS.

CONCLUSIONS
We performed the MD simulation of the PF MGB complexes with DNA oligonucleotides poly(dA)·poly(dT) and poly(dCG)·poly(dCG). The stability of those complexes is responsible for the timescale of the PF intercalation reaction [25]. The results of our simulation confirm, that PF MGB complexes are stable in the 50 ns time interval but there are some structural rearrangements with respect to the initial docking structures, especially in the case of PFpoly(dCG)·poly(dCG) complex.
The binding of PF to the DNA minor groove is sequence specific. The PFpoly(dA)·poly(dT) MGB complex is enthalpically more favorable than that with the poly(dCG)·poly(dCG). It is stabilized by the van der Waals and hydrophobic interactions. Besides, there are specific hydrogen bonds in PF-poly(dA)·poly(dT) complex formed by the protonated ring nitrogen and one of the amino-groups from the PF side and O2 thymine and N3 adenine atoms of poly(dA)·poly(dT), which were predicted earlier by Marcandalli [15]. The PF-poly(dCG)·poly(dCG) MGB complex is weaker due to the interfering amino-groups of guanine in the minor groove of poly(dCG)·poly(dCG), which prevent deeper binding of ligand.
The binding of PF launches the DNA conformational changes through the induced fit mechanism. The narrowing of the DNA minor groove in the PF binding site is observed. In the PF-poly(dCG)·poly(dCG) complex, it is more pronounced and is accompanied by the BI/BII backbone conformational transitions. During the induced fit process, the reorientation of ligand in the PF-poly(dCG)·poly(dCG) complex takes place: the PF chromophore becomes located along the minor groove, like in PF-poly(dA)·poly(dT) complex, except that its aminogroups are directed outside the groove.
In the case of PF-poly(dA)·poly(dT) complex, the conformational changes are not so large but they are still important for tighter PF-DNA binding: the change of the propeller twist and the increased opening in the PF binding site optimize specific hydrogen bonds between the PF and poly(dA)·poly(dT). Recently Ramakers and coworkers reported that the binding of Hoechst33258 to the A-tract DNA induces subtle conformational changes in the DNA structure, which maximize favorable ligand-DNA interactions [57]. The results of our study and that of Ramakers indicate that the model of rigid DNA-target is not a good choice even for small minor groove binders. An induced fit model gives better description of the binding process. One should keep in mind the DNA conformational changes induced by ligand binding upon design of new minor groove binders.   All energy terms are given in kcal/mol. ST and MT are the single and multiple trajectory approaches, correspondingly. ΔG conf is the sum of the bonded energy terms (bond, angle, dihedral) and 1-4 non bonded interactions; ΔG vdW is the van der Waals energy; ΔG Coul +ΔG PB -is the sum of the Coulombic energy and the polar solvation energy (Poisson-Boltzmann energy); ΔG rep -is the repulsive cavity term; ΔG disp -is the dispersion energy term describing the attractive solute-solvent interactions; ΔH is the sum of all energy terms except entropy contribution; -TΔS is the entropy penalty due to the loss of translational, rotational and vibrational degrees of freedom upon ligand binding; ΔG Total = ΔH-TΔS.