## Abstract

The formation of bubbles in nucleic acids (NAs) is fundamental in many biological processes such as DNA replication, recombination, telomere formation and nucleotide excision repair, as well as RNA transcription and splicing. These processes are carried out by assembled complexes with enzymes that separate selected regions of NAs. Within the frame of a nonlinear dynamics approach, we model the structure of the DNA duplex by a nonlinear network of coupled oscillators. We show that, in fact, from certain local structural distortions, there originate oscillating localized patterns, that is, radial and torsional breathers, which are associated with localized H-bond deformations, reminiscent of the replication bubble. We further study the temperature dependence of these oscillating bubbles. To this aim, the underlying nonlinear oscillator network of the DNA duplex is brought into contact with a heat bath using the Nosé–Hoover method. Special attention is paid to the stability of the oscillating bubbles under the imposed thermal perturbations. It is demonstrated that the radial and torsional breathers sustain the impact of thermal perturbations even at temperatures as high as room temperature. Generally, for non-zero temperature, the H-bond breathers move coherently along the double chain, whereas at *T*=0 standing radial and torsional breathers result.

## 1. Introduction

The vital life processes of nucleic acids (NAs) have increasingly been addressed and become intensively studied over recent years by biologists, chemists, mathematicians and physicists. The interest of the latter is mainly focused on the actual dynamical processes involved in information retrieval from the NA duplexes associated with structural transitions to single-stranded sequences. Consequently, the two complementary strands of double‐helix DNA must be separated (or at least locally melted) and rearranged in order to yield single-stranded sequences to be used by other molecules, requiring the breaking up of the hydrogen bonds in the corresponding region of the double helix. Therefore, the opening of a select region of the twisted (Watson–Crick) DNA double helix begins with a partial unwinding at an area called the replication fork, and is observed as a bubble. In general, the unwinding of NAs is achieved by a type of enzyme belonging to the large family of NA helicases (see Stryer *et al*. 2002; more than 60 different types).

The latter are ubiquitous and versatile enzymes that, in conjunction with other components of the macromolecular machine, carry out important biological processes, such as DNA replication; DNA recombination; nucleotide excision repair (i.e. UV damage); RNA editing and splicing; the transfer of a single-stranded nucleic acid (ssNA) to another ssNA or a protein or its release into solution (von Hippel & Delagoutte 2001); and the formation of T-loop telomeres (Wu & Hickson 2001). The NA helicases have several structural varieties (monomeric, dimeric, trimeric, tetrameric and closed hexameric), but all of them use the hydrolysis of nucleoside triphosphate to nucleoside diphosphate as the preferred source of energy. Our present study deals with a nonlinear dynamics model of the formation of oscillating bubbles in regions of duplex DNA acted upon by an enzyme. Utilizing models based on nonlinear lattice dynamics to study the formation of DNA opening has been of great interest recently (Englander *et al*. 1980; Agarwal & Hennig 2003). The existence of localized modes, such as solitons and breathers (describing energy localization and coherent transport), makes the nonlinear lattice approaches appealing. With a view to vibrational excitation in DNA, the Peyrard–Bishop model (Peyrard & Bishop 1989) and its successors (Dauxois *et al*. 1993; Barbi *et al*. 1999*a,b*; Cocco & Monasson 1999) have been successfully applied to describe moving localized excitations (breathers) which reproduce typical features of the DNA opening dynamics, such as the magnitude of the amplitudes and the time-scale of the oscillating ‘bubble’ preceding full strand separation.

We are interested in the bubble formation process initiated by structural deformations of selected regions of the parental DNA duplex that serves as a template for the replication. The process begins when the replication apparatus identifies the starting point, and then is combined to it. The replication apparatus is an assembled complex that forms the replication fork and opens it directionally. During this process, the helical DNA is unwound and replicated. One of the components of that complex is helicase, an enzyme which is unable to recognize the origin of the replication *per se*, and which requires the participation of specific proteins to lead it to the initiation site (Delagoutte & von Hippel 2003).

It is assumed that this enzyme operates in that way: a local unwinding in combination with (rather small) stretching of the H bonds in this region. We aim to demonstrate that there are indeed initial deformations which give rise to localized vibrations, constrained to a region of DNA, matching the properties of oscillating bubbles observed experimentally in DNA. These oscillating bubbles, with their temporarily extended but yet unbroken H bonds, serve as the precursors to the replication bubble. Furthermore, we are also interested in whether the stable radial and torsional breathers persist under imposed thermal perturbations.

## 2. Methods

The nonlinear oscillator network model for the DNA double helix used in this paper is explained in detail in several studies (Barbi *et al*. 1999*a*; Cocco & Monasson 1999; Agarwal & Hennig 2003; Hennig & Archilla 2004). The equilibrium position of each base within the duplex configuration is described in a Cartesian coordinate system by , and . The index pair (*n*, *i*) labels the *n*th base on the *i*th strand with *i*=1, 2 and 1≦*n*≦*N*, where *N* is the number of base pairs considered. Displacements of the bases from their equilibrium positions are denoted by *x*_{n,i}, *y*_{n,i} and *z*_{n,i}. The potential energy, taking into account the interactions between the bases, consists of four parts. The potential energy of the hydrogen bond within a base pair is modelled typically by a Morse potential,(2.1)where the variable *d*_{n} describes dynamical deviations of the H bonds from their equilibrium lengths *d*_{0} (see Hennig & Archilla 2004). The site-dependent depth of the Morse potential, *D*_{n}, depends on the number of involved H bonds for the two different pairings in DNA, namely the G–C and the A–T pairs. The former pair includes three H bonds, while the latter includes only two. *α*^{−1} is a measure of the potential‐well width.

The energies of the rather strong and rigid covalent bonds between the nucleotides *n* and *n*−1 on the *i*th strand are modelled using harmonic potential terms,(2.2)*l*_{n,i} describes the deviations from the equilibrium distance between two adjacent bases on the same strand. *K* is the elasticity coefficient.

Effects of stacking, which, owing to the backbone rigidity, impede the sliding of one base over another (Stryer *et al*. 2002), are incorporated in the following potential terms:(2.3)

The supposedly small deformations in the longitudinal direction can be modelled by harmonic elasticity potential terms given by(2.4)

The kinetic energy of a nucleotide is determined by(2.5)where *m* is the mass and denotes the (*x*, *y*, *z*)-component of the momentum.

The model Hamiltonian then becomes(2.6)with(2.7)and the summation in equation (2.6) is performed over all nucleotides and the two strands.

For a typical equilibrium configuration, the rotation angle, by which each base is rotated around the central axis, is given by *θ*_{0}=36°, the distance between base pair planes is *h*=3.4 Å and the inter-base distance (the diameter of the helix) is *d*_{0}=20 Å. For the average mass of one nucleotide, we use *M*=300 amu=4.982×10^{−25} kg. The arbitrary base pair sequence is reflected in randomly distributed bivalued H-bond coupling strengths *D*_{0} (for A–T) and *D*_{1}=2*D*_{0} (for G–C). Note that a ratio of 1.5 is often used, owing to the fact that the G–C base pair has three hydrogen bridges and the A–T two, but the energy depends on the angles and distances of the H-bonds, and on the polarity of the bases, among other factors. Quantum chemical calculations (Šponer *et al*. 2001) lead to a ratio of 2, as used here. However, a ratio of 1.5 leads to qualitatively similar results. The aperiodic arrangement of the two different base pairs, coding the genetic information, renders the DNA chain to be of *A–B-disorder* type. Following Barbi *et al*. (1999*a*) we set *α*=4.45 Å^{−1}, *D*_{0}=0.04 eV, and *K*=1.0 eV Å^{−2}. The value of the parameter *S* is reasonably taken as *S*=2*K* (Cocco & Monasson 1999), and a credible value for *C* is given by *C*=*S*/10 (Agarwal & Hennig 2003).

With the time-scale , one passes to a dimensionless formulation with quantities(2.8)

(2.9)

(2.10)

(2.11)

Subsequently, the tildes are dropped.

For the temperature dependence, we assume that the system is in contact with a heat bath which, based on molecular dynamics techniques, is simulated with the Nosé–Hoover method (Nosé 1984*a*,*b*; Hoover 1985). According to this method, the DNA lattice system is coupled to two additional variables, *p*_{s} and *s*, such that the equations read(2.12)

(2.13)

(2.14)

(2.15)

(2.16)

(2.17)

(2.18)

(2.19)

Note that the heat bath variable *s* does not explicitly enter the equations of motion regarding the actual dynamics of the bases. On the other hand, is an augmented Hamiltonian, the conservation of which can be used to assure accuracy of the numerical computations. The parameter *Q* determines the time-scale of a thermostat with temperature *T*, and *k*_{B} is the Boltzmann constant.

The values of the scaled parameters are *K*=0.683, *r*_{0}=44.50, *h*=15.13 and *l*_{0}=31.39. One time unit of the scaled time corresponds to 0.198 ps of physical time.

We turn now to the study of the creation of bubbles in DNA. The starting point is a DNA molecule for which a certain segment experiences initially angular and radial deformations owing to the action of some enzyme to which a region of the DNA is bound. In order to simulate the deforming action of enzymes, we assume that initially a number of consecutive sites in the centre of the DNA lattice (hereafter referred to as the *central region*) are subjected to forces acting in angular and radial direction, such that in this region the molecule experiences twist reduction together with radial stretching. These structural deformations can be extended over a region encompassing up to 30 bp and, as is going to be demonstrated, give rise to the formation of H-bridge breather solutions (extending over 15–20 bp), reproducing the oscillating ‘bubbles’ observed for the DNA-opening process (Barbi 1998). In figures 1 and 2, the localized initial distortions are shown. Both the angular and radial deformation patterns are bell-shaped, owing to the fact that the enzymatic force is exerted locally, i.e. in the extreme case to a single base pair for which only the H bond is deformed.1 The first one, with non-positive amplitudes, is linked with reduced twist, while the latter one, with non-negative amplitudes, is associated with radial stretchings. The radial and angular deformation patterns are centred at the central lattice site (base pair) at which the H-bond stretching and twist-angle reduction is at a maximum. On either side of the central site, the amplitudes progressively approach zero. The deformation energy amounts to 0.0362 eV. The set of coupled equations (2.12)–(2.19) is integrated with the help of a fourth-order Runge–Kutta method. For the simulation, the DNA lattice consists of 799 sites, and open boundary conditions were imposed. The same initial conditions were used both at zero and non-zero temperatures.

## 3. Results

First, we consider the zero-temperature regime. In figure 3, we depict the spatio-temporal evolution of the distance, *d*_{n}(*t*), between two bases of a base pair, which measures the variation of the length of the corresponding H bond, expressed in ångströms. The dynamics start from a non-equilibrium configuration so that energy redistribution takes place. Most noticeably, this is established in the immediate emission of (small-amplitude) phonon waves from either side of the localized radial pattern situated around the central lattice site. These primary phonon waves travel uniformly towards the ends of the lattice with uniform velocity. At the approach of an equilibrium regime, a small amount of excitation energy is dispersed in the form of secondary phonon waves in the rest of the DNA lattice. However, the vast majority of the excitation energy remains contained in the central region.

We observe that the localized bell-shaped radial pattern is preserved and performs periodic oscillations in time, that is, a radial breather is formed. This periodically oscillating pattern of the radial variable, *d*_{n}(*t*), is attributed to successive stretchings and compressions of the H bonds.

Note that the stretching of a base-pair distance is larger than the compression typical for the evolution in a Morse potential (see also Barbi et al. 1999*a*).

With regard to the associated pattern of the deviations of the twist angles from their equilibrium values, *θ*_{n}(*t*)−*nθ*_{0}, expressed in radians, we find that in the initial phase the amplitudes increase. With the emanation of the primary radial phonon waves from the localized radial pattern (cf. figure 3), there is a linked split up of the standing angular bell-shaped pattern, so that two primary torsional waves of negative localized angular deformations, symmetrical to the central site, are produced, which travel in unison with the radial phonon waves in the direction of the lattice ends. Since the radial phonon waves possess positive amplitudes corresponding to H-bond stretching, the region of the double helix traversed by them experiences a reduction of the twist angle, i.e. *θ*_{n}(*t*)−*nθ*_{0} is negative. Similarly, the second set of radial phonon waves are connected with two small negative-amplitude angular waves following the primary ones away from the central region. Around the central lattice site, the (initially negative) amplitudes of the twist angle grow steadily. Finally, after nearly 150 time-units, the twist angles of the helix in the breather region become positive, this being equivalent to increased twist in a region comprising 20 lattice sites to either side of the central site. Note that the process of amplitude breathing of the torsional part of the lattice proceeds with a longer period than the breathing of the radial localized pattern. This behaviour of *θ*_{n}(*t*)−*nθ*_{0} for zero temperature is shown in figure 4.

Remarkably, in the non-zero-temperature regime, when thermal energy is injected into the DNA lattice, the radial breather essentially remains in a localized shape. Moreover, its amplitudes and energy grow. Therefore, the breather extracts energy from the thermal bath. Simultaneously, we observe that excess energy, not to be carried by the radial breather, is ejected from the central region, where it disperses into wider parts of the lattice as phonons. Nevertheless, a radial breather of fairly large amplitude prevails over the small-amplitude noisy background caused by the heat bath on the DNA lattice. Interestingly, under the impact of thermal perturbations, the radial breather starts to move towards one end of the DNA lattice in a coherent fashion, with the maximum amplitude being larger than the *T*=0 case. This feature is illustrated for the case where *T*=100 in figure 5. Correspondingly, the propagating radial breather is accompanied by a torsional breather leaving the double helix consecutively in over-twisted and under-twisted shape. The direction of the movement of the breather is determined by the actual realization of the random sequence of the DNA base pairs. This is reflected in our model with the bivalued dissociation levels *D*_{1}=2*D*_{0} that are assigned to the sites in a random fashion along the backbone, so that the translational invariance of the lattice is broken by disorder. For the simulations presented here, which were performed with the same realization of disorder for all considered temperatures, the breather moves towards the right end of the DNA chain. However, we found that there also exist realizations of the random *D*_{0}–*D*_{1} sequence for which the motion is towards the left end. The fact that the breather movement is not diffusive, i.e. has random changes of directions, means that the moving breather is a stable new entity produced by the perturbation of the standing one. In mathematical terms, the thermal bath brings about the crossing of a separatrix between the basins of attraction of the two breathers. Varying the temperature yields qualitatively equal results, i.e. the radial and torsional breathers become eventually mobile, travelling towards one end of the DNA lattice. Generally, it holds that the higher the temperature, the more the amplitudes of the breather grow and the more strongly the helix unwinds. For ambient temperature, i.e. *T*=300 K, the maximal amplitude of the radial breather is ten times larger than that for the zero-temperature case, being related to the stretching of the associated H bond by 0.6 Å away from its equilibrium value. Moreover, at the sites (base pairs) where the helix is radically stretched to the greatest extent, it is unwound by an angle of 12°.

Quantitatively, the results regarding the temperature dependence of breather evolution are suitably summarized by the time evolution of the first momentum of the energy distribution, defined as(3.1)where the energy *E*_{n,i} is defined in equation (2.7) and *n*_{c} is the site index corresponding to the centre of the DNA lattice.

This quantity describes the temporal behaviour of the position of the centre of a breather. Thus, it represents a measure for the mobility of the breathers. Generally, the higher the temperature, the larger the amplitude of the radial breather and the faster the radial (and torsional) breather travels along the DNA lattice. This is illustrated in figure 6. As far as the degree of energy localization in the breathers is concerned, the energetic partition number serves as an appropriate mean. It is defined as(3.2)*P* quantifies how many sites are excited to contribute to the breather pattern. In the *T*=0 case, the partition number changes only slightly and oscillates in an interval rather close to its initial value at all times. The relatively small growth of *P*(*t*) accounts for the phonons being emitted from the central region.

For *T*>0 we find that *P*(*t*) rises with time. While up to times *t*≲80 the mean growth of *P*(*t*) proceeds linearly, afterwards, when the DNA lattice system comes closer and closer to equilibrium, the participation number increases only in a logarithmic fashion. Eventually, the value of *P*(*t*) reaches a plateau around the time value of 150. Compared with the initial width of the localized distortion pattern, i.e. *P*(0)=29, the energy spread has risen by ≲44%. This energy distribution over the lattice, however, is mainly due to the thermalization of the lattice by the thermal perturbation, rather than by an actual spread of the radial breather. The latter maintains its degree of localization at all times. The time evolution of *P* is shown in figure 7. Remarkably, in the range of 50–300 K, we observe almost equal temporal evolution of *P* regardless of the value of *T*. Thus, the thermal perturbations do not strongly influence the width of the breather anymore—evidence of the stability of the breather and the pronounced storage capability of the DNA lattice.

## 4. Discussion

In summary, we observed that, out of an initial non-equilibrium situation, for which the hydrogen bonds of an under-twisted segment of the DNA lattice have been stretched, breathers develop in the radial and angular displacement variables. For zero temperature, the breathers remain standing in the initially excited region. Interestingly, for *T*>0, the breathers start to travel coherently towards one end of the DNA duplex. The observed breathers realistically represent the oscillating bubbles found prior to the complete unzipping of DNA, i.e. they oscillate with periods in the range 0.3–0.8 ps, are on average extended over 10–20 bp, and possess maximal amplitudes of the order of ≲0.6 Å. Our results demonstrate that the action of some enzyme, mimicked by localized radial and torsional distortions of the DNA equilibrium configuration, in fact initiates the production of oscillating bubbles in DNA. Moreover, these oscillating bubbles sustain the impact of thermal perturbations. Whether the bubbles, when meeting each other, can merge in such a way that a larger-amplitude radial breather results which resembles the replication bubble with broken H bonds is still being investigated.

## Acknowledgments

The authors acknowledge support by the Deutsche Forschungsgemeinschaft via a Heisenberg fellowship (He 3049/1-1). J.M.R. acknowledges a Study Licence from the Junta de Andalucia, Spain.

## Footnotes

↵It turns out that this shape is also more stable in a thermal bath than a rectangular shape, as in Hennig & Archilla (2004).

- Received June 16, 2004.
- Accepted November 16, 2004.

- © 2005 The Royal Society