Capture of a Transition State Using Molecular Dynamics: Creation of an Intercalation Site in dsDNA with Ethidium Cation

The mechanism of intercalation and the ability of double stranded DNA (dsDNA) to accommodate a variety of ligands in this manner has been well studied. Proposed mechanistic steps along this pathway for the classical intercalator ethidium have been discussed in the literature. Some previous studies indicate that the creation of an intercalation site may occur spontaneously, with the energy for this interaction arising either from solvent collisions or soliton propagation along the helical axis. A subsequent 1D diffusional search by the ligand along the helical axis of the DNA will allow the ligand entry to this intercalation site from its external, electrostatically stabilized position. Other mechanistic studies show that ethidium cation participates in the creation of the site, as a ligand interacting closely with the external surface of the DNA can cause unfavorable steric interactions depending on the ligands' orientation, which are relaxed during the creation of an intercalation site. Briefly, such a site is created by the lengthening of the DNA molecule via bond rotation between the sugars and phosphates along the DNA backbone, causing an unwinding of the dsDNA itself and separation between the adjacent base pairs local to the position of the ligand, which becomes the intercalation site. Previous experimental measurements of this interaction measure the enthalpic cost of this part of the mechanism to be about −8 kcal/mol. This paper reports the observation, during a computational study, of the spontaneous opening of an intercalation site in response to the presence of a single ethidium cation molecule in an externally bound configuration. The concerted motions between this ligand and the host, a dsDNA decamer, are clear. The dsDNA decamer AGGATGCCTG was studied; the central …GATG… site was the intercalation site.


Introduction
DNA is a highly flexible and dynamically responsive molecule. Numerous cellular metabolic pathways, which include transcription and/or translation of DNA, necessarily must respond to the local and temporal environment of the cell, and this necessitates a tightly controlled interaction between DNA and these pathways and networks that act to regulate it. To respond quickly and precisely, DNA must react to participate in molecular interactions in ways that are precise, diverse, and controlledly reversible.
The flexibility of this biomolecule is central to its responsiveness. The DNA backbone, a polymer made up of repeating subunits of alternating deoxyribose sugar rings and phosphate moieties, is quite flexible about its single bonds. The base pairs, tethered at one end by a rotatable glycosidic bond to a sugar, may each be imagined as a somewhat flat "paddle", tethered across the helix to its partner by either two or three hydrogen bonds, each of which has a small but nonnegligible energetic cost to rupture and reform. The base pairs constrain the phosphate backbone, which by itself would be a linear (and informationless) chain, not found in the familiar twisted helical formation. This is not only due to the rigidity of the bases but also due to their hydrophobic nature, which makes it energetically favorable to reduce the nucleobase exposure to solvent waters by the stacking of the bases between the outer backbone chains.
As a tool to delineate and visualize the interactions and motions available to this molecule, molecular dynamics has been used to study thermally equilibrated DNA (e.g., see [1][2][3]), as well as to understand interactions between this molecule and intercalating ligands that bind to it. Delineating the step-by-step mechanism of intercalation has been of interest since Lerman's discovery of the acridines 2 Journal of Nucleic Acids [4] and his subsequent proposal of their mechanism of interaction with DNA via intercalation. This mechanism has been well described in papers and review articles, and three major mechanistic approaches have been proposed, which differ in their kinetic microstates (references [5,6] describe single relaxation mechanisms, [7][8][9] two relaxations, [10,11] three relaxations). These mechanistic pathways differ in their kinetic microsteps in describing the mechanism of intercalation. In one of these proposed mechanisms, the "two-step" mechanism, the intercalator electrostatically binds externally to the DNA in a fast process, followed by a slow insertion step to the intercalation site (D = dsDNA, E + = ethidium cation): In the three-step mechanism, the intercalated complex forms a secondary intercalated complex, either via a geometrical change (such as a rotation when the molecule has two moieties by which it can intercalate) within the intercalator itself, or by . the intercalator moving from intercalation site to another intercalation site without re-entering the solvent: One concern with using molecular dynamics as a tool to fully describe the complex potential energy surfaces of molecular interactions is the time limitation imposed by current computational techniques. In general, an MD simulation of around 0.5 to 1 microsecond of a modestly sized system of protein or DNA (less than 5000 atoms) with solvent (either implicit or explicit) is the length that researchers with access to computational capabilities such as a supercomputer or reasonably sized CPU cluster can reasonably carry out. Yet many interesting dynamical features of such trajectories may only be observed on longer timescales, as long as seconds or in certain cases minutes. A notable feature of such dynamics simulations is that most of the time during the simulation the system is undergoing nonproductive random motions which only result in a reaction if certain energetic barriers are overcome. Transition state theory shows that reactants must both acquire energy from their surroundings as well as be in geometrically favored configurations in order to react. Only this situation results in a transition state that may proceed to product (or back to reactants).
Most of the time a given molecule in any system is either in its reactant or product state, as the transition state itself is usually fleeting and occurring on the order of nanoseconds or even femtoseconds. Yet the transition state is precisely the state with the most information about the mechanism of how the reaction proceeds. In long MD simulations, in order to learn anything about the mechanism it is essential to reach a transition state.
A computational simulation due to hardware or software constraints, or other real-world limitations, may not be able to run long enough to reach such a transition state, but there are methods by which a system may be driven to a transition state, for example if a particular mechanism is hypothesized and is being tested for viability.
It must be appreciated that undriven dynamical simulations of interactions between molecules, such as between a ligand and drug, typically yield trajectories which include very long stretches of nonproductive motions from which no mechanistic information may directly be gleaned. Of course, flexibility, movements, and positions of side-chains, excursions of particular molecular moieties and/or loops, and other conformational and geometric features may be observed from such trajectories [12], and this information is useful indeed in characterizing and/or understanding the molecule's inherent flexibility and features as well as in proposing or ruling out possible mechanisms of interaction with known ligands. Therefore, molecular dynamics alone is not usually the tool of choice for direct observation of a molecular mechanism.
To improve the ability of MD simulations to capture interesting transition state structures or simulate processes on a longer timescale, a number of methods, such as adiabatic melting, cross-correlation analysis [13], replica exchange molecular dynamics (REMD) [14], and forced motions along a known trajectory [15], have been developed. Some of these methods require prior knowledge of the mechanism being studied, and although guesses may be made to test hypothetical transitions, the results would then be suggestive but not conclusive. Inspection of the movie shows a remarkable concerted interaction between the ethidium and the target intercalation site on the dsDNA decamer 5 AGGATGCCTG3 , which was the central . . .GATG. . . sequence. This interaction was not driven by external simulation parameters such as constraints, applied force constants, or any other method. As this trajectory was not replicable, it remains an interesting phenomenon.

Results and Discussion
This data could be an artifactual simulation run, representing a nonphysical system. In fact, the b-DNA structure has become very slightly distorted and the helix is widened by almost 1 angstrom. However, inspection of the potential, kinetic, and total energy during the equilibration stage of the MD run shows these energies to be stable. No undue force or steric strain is found in the complex, and the final structure is stable.
Viewing this trajectory one can see a remarkable concerted interaction between the ethidium and the DNA. At 220 femtoseconds (fs), the AT base pair which forms the intercalation site begins to disrupt their interbase pair hydrogen bonds by altering their positions with respect to one another, accompanied by the ethidium approaching the Journal of Nucleic Acids helix from its initial position. At 1300 fs, the phenyl group on the ethidium can be seen to have twisted, followed by the rest of the ethidium molecule twisting to be perpendicular to the helical z-axis of the DNA. At 3700 fs, this twist is completed and the ethidium has rotated to be in the same plane as the plane of the space between the base pairs. At this point, the ethidium begins sliding towards this nascent intercalation site, which is not yet open. At 4100 fs, the base pairs are disrupted but not open into the intercalation configuration. The ethidium cation can be seen moving around the site. At 5600 fs, the base pairs can be seen to open in an asymmetric fashion. At 7200 fs, the complex has fully intercalated. The balance of the dynamics run is seen to simply be the thermal fluctuations available to this complex. This astonishing simulation shows much about a possible mechanism for intercalation. The experimental conditions represented in this simulation are low concentration of ligand (ethidium), physiological pH, low salt, STP, and aqueous solvent. Previous experimental studies done under these particular conditions [5,6] favor the two-step mechanism. This simulation likewise shows a two-step mechanism, where the ethidium does not move away from its externally bound position but twists in place and slips into an intercalation site, that is, created in response to steric and electrostatic forces resulting from the twisting motion of the charged ethidium. Therefore, this simulation provides support for the two-step mechanism of ethidium intercation over the threestep mechanism.

Experimental Section
The original purpose of these calculations was to look at the stability of a proposed major-groove outer-bound complex between ethidium and DNA, and this particular calculation was the initial run to test this hypothesis. Ten other runs were done with similar randomly placed starting configurations for the externally bound ethidium, and in all such runs the externally hound ethidium equilibrated and remained in an externally hydrogen-bound state as described in [16]. In addition, the exact starting configuration that is the topic of this paper was used in three subsequent dynamics runs, and 10 additional runs with the externally bound ethidium placed within 2 angstroms RMSD of the initial position in the configuration which is being reported in this paper were also completed. None of those runs replicated the intercalative trajectory found in this particular run.
The commercial software Sybyl [17] was used to prepare this structure. Ethidium charges were imported from MOPAC. These calculated charges were compared against a detailed experimental study that determined the charges on the heavy atoms of ethidium [18]. Because charges on the hydrogens were absent from this paper, MOPAC was used and then comparisons were made to assure that the charges assigned by this method were physically reasonable. The bform dsDNA sequence 5 AGGATGCCTG3 was constructed in Sybyl and minimized using the Amber99 forcefield. The DNA was immersed in a 30-angstrom waterbox, using 2253 TIP3P waters. Sodium counterions were initially placed at a distance of 5 angstroms from and bisecting the P-O-P bond of the phosphate backbone. These were allowed to relax individually during nested minimizations to assist physiochemically reasonable placement.
Initially, the waters were allowed to relax while the solute molecules (DNA, ethidium cation, and counterions) were held by harmonic constraints. After this converged to an RMSD gradient of 0.5 angstroms, the solvent waters were constrained while the rest of the system was allowed to minimize. This method of nested minimizations was carried out for 3 rounds. A final round of minimization with the entire system unconstrained was run until the criterion of an RMSD of less than 0.1 angstrom was attained. The density of the waterbox did not change over the course of the minimization.
A 500 ps dynamics run was then started from this fully minimized complex, using an NVT ensemble, SHAKE off, to a final temperature of 300 K.

Conclusions
In this paper an anomalous molecular dynamics simulation which seems to have captured the transition state of an intercalative interaction between ethidium and a dsDNA decamer is presented. The initial placement of the ligand in this particular simulation, along with a random initial Boltzmann seed and the particular consequent distribution of energies in this run, has allowed the observation of a previously proposed but not yet observed series of kinetic microsteps along the mechanistic pathway of ethidium intercalation-the concerted opening of an intercalation site as the direct effect of the presence of this ligand.