## Abstract

Single molecule DNA experiments provide interesting data that allow a better understanding of the mechanical interactions between the strands and the nucleotides of this molecule. In some sense, these experiments complement the classical ones about DNA thermal denaturation. It is well known that the original Peyrard–Bishop (PB) model by means of a harmonic stacking potential and a nonlinear substrate potential has been able to predict the existence of a critical temperature of full denaturation of the molecule. In the present paper, driven by the findings of single molecule experiments, we substitute the original harmonic intra-strand stacking potential with a Duffing type potential. By elementary and analytical arguments, we show that with this choice it is possible to obtain a sharp transition in the classical domain wall solution of the PB model and the compactification of the classical solitary wave solutions of other models for the dynamics of DNA. We discuss why these solutions may improve our knowledge of the DNA dynamics in several directions.

## 1. Introduction

The celebrated Peyrard–Bishop (PB) model to describe DNA dynamics was created as a first attempt to model the *denaturation* of the double-stranded macromolecule (i.e. the partial or even total strand separation). It is well known that this opening phenomenon, which corresponds to the disruption of the H-bonds joining the bases in pairs, is the first step in a chain of events leading to the expression of the genetic information encoded in the double helix. For this reason, our understanding of denaturation is a fundamental issue on the long way toward the unravelling of the DNA secrets.

In a physiological framework, denaturation occurs, for example, in the *transcription* process when the DNA unwinds into separate sugar-phosphate strands over a short region (15–20 bp). This opening action takes a relevant amount of energy to pull the base pairs apart and expose them to water. To synthesize the nascent RNA chain, it is necessary to consider the role of the RNA polymerase, the enzyme copying DNA to RNA, and to study the corresponding *bubble* dynamics. Englander *et al*. (1980) were the first ones proposing to list, among the various goals of the RNA polymerase, the ability to focus thermal energy over a precise region before the transcription, such that the double helix can be opened. Moreover, Englander and co-workers conjectured that the focusing may be achieved by means of *solitons* propagating on the double helix. In the last 25 years, a large amount of research has been devoted to the nonlinear dynamics of DNA and to the quest for coherent structures in these models, as it has been clearly reported in a recent review by Peyrard (2004) or in Yakushevich's book (2004), and still today the research in this field is very active.

In many of the nonlinear models for DNA dynamics, the denaturation bubble is described either as a breather, a pulse or a kink solitary wave solution of a *semilinear* Klein–Gordon type equation. These equations are obtained by considering what we call a *mesoscopic* model of DNA. Indeed, the possibility to simulate the dynamics of DNA by molecular modelling on the relevant timescales of transcription and folding is computationally expensive and not always feasible. Transcription involves a region that extends about 20 bp, and therefore to study this phenomenon with a molecular model we need to consider at least 100 bp. In 2004, a typical simulation protocol was able to analyse oligomers of 18 bp (in the presence of water and Na^{+} counterions plus 100 mM NaCl) for 20 ns with snapshots every picosecond (Lankas *et al*. 2003). Moreover, molecular computations generate a huge amount of data that would be very difficult to synthesize in mesoscopic quantities of practical interest. This is the reason why micro-scale models seem to be prohibitive, and therefore it seems necessary to address our attention to a scale intermediate between the atomic scale and the full macromolecule scale models where the DNA is modelled as a continuum rod (Cuesta-Lopez *et al*. 2005).

Mesoscopic models of nonlinear DNA dynamics are usually formulated considering a Hamiltonian that takes into account some of the various possible modes of deformation of the macromolecule (e.g. stretching, twisting and bending) plus a substrate potential. This potential takes into account, for example, the interaction between the two strands. In the present paper, we restrict our attention to the simplest models proposed for the DNA dynamics, in order to avoid complexities that may hinder the mathematical feasibility of the model.

The aim of the present paper is to investigate in more detail the role of *nonlinear stacking forces* in the mathematical issues associated with wave propagation in the PB and other basic models of DNA dynamics. We remember that the essential role of nonlinear effects in the statistical mechanics of the PB model has been already pointed out in Dauxois *et al*. (1993), where an improvement, including cooperativity effects through anharmonic nearest-neighbour stacking interactions, has been proposed.

Our main goal is to show by simple analytical arguments that the anharmonic stacking is indeed the main reason for sharp denaturation interfaces and for *compactification* of the various solutions of interest in DNA dynamics. Indeed, since observed patterns in nature are usually of finite extent, the usual modelling of a coherent localized structure via a soliton with infinite extent is a shortcoming of a continuum theory and of an inadequate modelling of the physical phenomena (Rosenau 1994). In this paper, we show that a careful mathematical modelling of the interplay between nonlinearity and dispersion is a necessary step toward a more realistic and hopefully useful generation of mesoscopic models for nonlinear DNA dynamics.

In the past, longitudinal soliton-like equations for DNA dynamics have been treated both with the continuum model of a rod and with the discrete Toda model (e.g. Christiansen *et al*. 1990). These studies predicted that, at the physiological temperature, a large number of longitudinal solitons may be indeed (theoretically) generated, but, since the experimental verification of such results was very hard, the main efforts have been addressed toward the transverse motions of the molecule. In fact, transverse motions are easier to observe because they are directly related to the breaking of the H-bonds, which is easily observable, for example, using UV absorbance. Xiao *et al*. (1987) have considered the influence of longitudinal vibrations in the well-known models of DNA dynamics proposed by Takeno & Homma (1983) and Yomosa (1983). The conclusion by Xiao *et al*. was that the internal DNA stretching motions were second-order phenomena. This result is coherent with the well-known fact that stacking interactions are stiffer than the transverse interactions of the hydrogen bonds. We point out that, also at the macroscopic level, nearly all the interesting results in the study of DNA configuration by using the classical theory of elastic rods have been obtained by using the Euler–Kirchhoff's model of the elastica. This model is based on the postulate of *inextensibility*, and the relaxation of this kinematical assumption in the classical theory of beams and rods is not a trivial issue (Antman 1995).

For all these reasons and for using linear stacking potentials, still today the standard approach neglects the stretching deformation mode also in highly refined models of DNA dynamics. This seems to clash with the results by Dauxois *et al*. (1993), who have shown the importance of anharmonic stacking interactions for a realistic description of the thermodynamics of the denaturation. The nonlinear stacking term introduced by Dauxois *et al*. has played a fundamental role in the context of the *vexata quaestio* about the nature of the phase transition in connection with the denaturation critical temperature (see Kafri *et al*. 2000; Theodorakopoulos *et al*. 2000). Moreover, in Barbi *et al*. (1999*a*,*b*), the basic PB model have been generalized in a significant way to take into account the twist–opening interactions due to the helicoidal molecular geometry. These results have then been modified in Cocco & Monasson (1999), with the introduction of fluctuations in the axial distance between successive base-pair planes. Nevertheless, in all these fundamental papers, the role between nonlinearity and localization of waves is not discussed in detail, and, to the best of our knowledge, only the recent paper by De Luca *et al*. (2004) addresses this important issue for similar models. It is worth to note that the results in De Luca *et al*. (2004) are based on numerical simulations and not on a rigorous mathematical discussion as in the present research.

On the other hand, the advent of single molecule experiments has been a deep revolution in the knowledge of the mechanics of DNA (Lavery *et al*. 2002; Strick *et al*. 2003). Hence, the large amount of experimental data coming from these manipulations allows to rethink again about the role of nonlinear stacking potential and of the intra-strand degrees of freedom. For example, it is now clear that the complex configuration of the biological molecule may be explained only by considering the nonlinear elasticity framework. Therefore, it seems no more possible to use a simple reductionist approach where all the various modes of deformation are uncoupled.

Starting from these considerations, in the present paper, we show by means of simple arguments that the role of the anharmonic elementary stacking potential is a consequence of the nonlinear mechanical properties of the DNA. Moreover, we show that this anharmonicity plays a fundamental role in the understanding of the localization of the energy along the biomolecule and in the prediction of the sharpness of the transition fronts between closed and open states. By using analytical methods, we establish that when a nonlinear stacking potential is considered, some of the classical solutions for the PB model clearly change their character. From the biological point of view, the interesting feature of these findings is that the usual solutions, as domain walls, breathers, kinks or pulses, present a strong compactification. Indeed, our new solutions may be special cases of the *compacton* solitary wave solutions introduced by Rosenau & Hyman (1993). These solutions are permanent nonlinear structures coming out from a more faithful description of the original problem in the discrete-to-continuum transition, and for this reason we think they are crucial for a proper description of DNA dynamics.

From the mathematical point of view, the important feature of the new solutions is that they are no more analytical, since they develop some discontinuity in the higher-order derivatives. It is not the first time that in the study of the DNA dynamics this sort of *weak* solutions are considered. For example, this is the case of the *peakon* solution of the rod-like model by Ichiwaka *et al*. (1981) and some of the numerical simulations provided by Takeno (2005), where the expression ‘compact-like modes’ is used explicitly. As far as we know, in these studies, the biological and mechanical role of the various weak solutions considered is not discussed in detail and it is not recognized that their appearance is strictly connected with the nonlinearity of the stacking interactions.

The plan of the paper is as follows. In §2, the basic equations are considered and the problem is discussed in more detail. In §3, we consider the static domain wall solution of the PB model. The nonlinear Schrödinger (NLS) equation obtained from the PB model and its modification when the anharmonic stacking potential is considered are given in §4. In §4, we also discuss when compact-like solitary waves with drop boundary conditions are possible in the PB model. In §5, we consider a different kind of model to show that anharmonic stacking interactions are responsible for the *compactification* of kink-like travelling waves. Section 6 is devoted to concluding remarks.

## 2. Basic equations

The PB model starts considering the transversal displacements of the nucleotides from their equilibrium positions along the direction of the hydrogen bonds. By introducing a transformation to the centre-of-mass coordinates representing the in-phase and out-of-phase transversal motions, the model reduces to a set of uncoupled dynamical equations. In this way, we have to consider only one degree of freedom per base pair, i.e. the stretching of the bond linking the *n*th base pair. The evolution of is determined by the Hamiltonian(2.1)where *m* is the reduced mass of the bases. When , we have a closed base pair; a positive value of indicates that the base pairs are stretched away and for a large value the denaturation may take place. In (2.1), the potential , describing the interaction between the two base pairs, is a Morse potential, i.e.(2.2)where *D* is the dissociation energy of the pair and *a* is a parameter which sets the spatial scale of the potential. The potential describes the interactions between adjacent bases along the DNA. In the original PB model, the following simplest harmonic approximation is considered:(2.3)Peyrard (2004, p. R12) points out that: ‘the harmonic approximation allows easier calculations and it is sufficient to get some interesting results which agree with some experimental observations’.

To compare the PB model with experimental data, Peyrard (2004) takes into consideration the fraction of base pairs that stay bounded at a given temperature. In so doing, it is possible to check that the model is able to predict the existence of the critical temperature showing the complete denaturation of the two strands. This success has been the main reason for a huge activity of research on mathematical models of the DNA macromolecule. Many generalizations of the PB model aiming to introduce more realistic geometries, a higher number of degrees of freedom, non-uniformity of the base pairs, dissipative and stochastic effects, non-local effects and so on, have been proposed in the last 15 years. A partial review of this intensive activity may be found in the Yakushevich's book (2004).

We have to say that one of the problems encountered from the beginning with the PB model was that the predicted curve fraction of closed base pairs versus the temperature was very smooth in contrast with the experimental data showing for the same curve a very sharp transition. The latter is a sort of catastrophic phenomenon, reminiscent of the *snap buckling* bifurcation in structural mechanics (see Antman 1995), indicating that experimentally the number of the denaturated base pairs increases rapidly in proximity of the critical temperature. To overcome the unphysical feature of the model, Peyrard has suggested to improve the description of the stacking interactions by considering, instead of the harmonic longitudinal coupling, the following improved anharmonic potential:(2.4)With this potential, the sharp transition to denaturation can be obtained and many important thermodynamical issues can find an adequate answer (see Theodorakopoulos *et al*. 2000). The choice of (2.4) was motivated by the observation that the stacking energy is not a property of single bases, but of the neighbouring couples of base pairs. Therefore, the expression (2.4) provides a phenomenological description of the cooperative effect: a base pair in the vicinity of an open state has a lower vibrational frequency, which reduces its contribution to free energy, but simultaneously a weaker coupling along the strands gives to the bases more freedom to move independently of each other, causing an entropy increase that drives a sharp transition.

How the base pairs stack on one other is a central point to determine the possible configurations of DNA and we know that subtle motions of the base pairs can be accumulated over series of such base pairs to make different kinds of double helices. This is an important point in the denaturation and then in the transcription process (Calladine & Drew 1997). The advent of techniques to micro-manipulate single molecules has enabled detailed studies of DNA elasticity and of the kinetics of motor enzymes. Recent experiments have measured the rate of replication of DNA catalysed by a single enzyme moving along a stretched template strand. One of the most significant results of such experiments is that it is possible to modulate the rate of DNA replication with the longitudinal mechanical tension (Goel *et al*. 2001). These experimental results clearly show the importance and the need of more refined models for the stacking potential. This need is more evident in a meso-scale model, where, what we call, the stacking potential synthesizes a huge amount of mechanical information.

Hence, of fundamental interest in our paper are the classical single molecule experiments aimed to compute the force-extension curves of double-stranded DNA. The mathematical interpretation of such experiments is usually given in terms of the worm-like-chain (WLC) interpolation formula and its generalizations (see Bouchiat *et al*. 1999; Bouchiat 2006; Ogden *et al*. 2006). The basic WLC interpolation formula is given analytically by(2.5)where *F* is the force applied, is the relative elongation, *L* the total contour length of the molecule, *A* its persistence length, the Boltzmann constant and *T* the absolute temperature. Obviously, in (2.5), . Formula (2.5) is in good qualitative and quantitative agreement with the experimental data in the weak-force regime (less than about 5 pN), which is a substantially linear regime and is in good qualitative agreement in the intermediate regime (from about 5 to about 60 pN), which is a highly nonlinear regime. The very large forces regime beyond the *overstretching transition* is not captured by (2.5).

The experiments connected with (2.5) show that in the elastic regime DNA is a nonlinear material. This means, that in pulling apart the base pairs, the isolation of a single mode of deformation is a strong approximation. Indeed, in nonlinear elastic materials, all possible modes of deformation (stretching, twisting, shearing, bending, etc.) are strongly coupled together. For example, any shear deformation (as torsion) cannot be obtained without normal stress differences (this is the so-called *Poynting effect*, see Ogden 1984). This means that a twist of the strands is always coupled with an extensional deformation of the molecule. These coupling effects are the reason why filamentous structures may show highly packed conformations and amazing solenoidal phases (Ghatak & Mahadevan 2005). Indeed, there is experimental evidence that the enzymes involved in the DNA transcription and replication, when they pull apart the two strands to originate the replication bubble, perform a very complex deformation of the macromolecule and therefore of the sugar-phosphate strand (Bates & Maxwell 2005).

In this framework, it is very hard not only to provide a precise quantitative description of the stacking interactions, but also to clearly understand the qualitative features of these interactions. In a mesoscopic model, and also in the simplest planar schema of the DNA molecule, the stacking potential has to synthesize a huge number of mechanical (and moreover electrostatic) phenomena. Even if (2.4) is more realistic than a harmonic potential, there is no clear evidence that this is the *real* stacking potential. At the moment, we ignore the form of the potential which is in qualitative and quantitative agreement with experiments.

A phenomenological approach first of all has to be driven by mathematical feasibility, in order to highlight the main differences between linear and nonlinear stacking interactions in a simple way. For this reason, in the following sections, we prefer to use the following *Duffing* polynomial form:(2.6)which may be considered the simplest anharmonic and *universal* interaction. Indeed, (2.6) may be used to approximate locally (by a Taylor expansion) any other regular nonlinear potential. The advantage of using (2.6) relies not only on its simple mathematical form. In fact, in appendix A, we show that the interaction (2.6) may be deduced from a phenomenological strain-energy compatible with the theory of nonlinear elasticity.

Therefore, by using (2.6), our main goal in §3 is to investigate the impact of the term in the nonlinear dynamics of DNA looking for exact closed-form results.

## 3. Compact domain wall between open and closed regions

By introducing the dimensionless stretching of the base pair *Y*=*ay* and the dimensionless time , in the continuum limit, the equation resulting from the PB Hamiltonian is given by(3.1)where . This is a classical example of *semilinear* Klein–Gordon equation corresponding to the one-dimensional motion for a linear elastic continuum on a nonlinear elastic substrate. It is well know that equation (3.1) has the exact static solution(3.2)where is an integration constant that determines the position of the wall. Solution (3.2) describes a DNA configuration, where for the molecule is closed, whereas for the base-pair separation grows linearly with space and the molecule is fully *denaturated*. The simple unique static solution (3.2) is a fundamental implication of (3.1), because as shown in Dauxois *et al*. (2002) (see also Peyrard 2004) this solution drives the transition by providing an appropriate balance between the energy cost and the entropy gain in the thermodynamics of the PB model.

Now we modify equation (3.1) by considering the anharmonic coupling between neighbouring sites within the same strand. By introducing the dimensionless variables, the nonlinear elastic potential is given, up to a constant, by(3.3)where and are two non-negative constitutive constants. Representative values of these parameters may be obtained, for example, by a local approximation of the experimental curves of the force-extension data of double-stranded DNA (see appendix A). In so doing, we have that in the weak-force regime and in the intermediate force regime. When the applied force is around 50 pN, we have that , and in this situation it is interesting to investigate what we call the extreme nonlinear case by considering . In such a way, the effects of nonlinearity are clearly identified.

### 3.1 The extreme nonlinear case

We first consider the *extreme* case and , obtaining the full nonlinear Klein–Gordon equation(3.4)In the static case where , from (3.4) we have the ordinary differential equation(3.5)and the consequent first integral(3.6)where *C* is an integration constant.

By considering the asymptotic boundary conditions , *Y*=0 for , we show that in (3.6) it must be . The positive real branch implicitly defined in (3.6) is given by(3.7)By performing the quadrature ( is a new integration constant), from (3.7) we have(3.8)and we obtain for *Y*>0(3.9)By denoting(3.10)in this case by equation (3.9), we have explicitly(3.11)In the case *Y*<0, from (3.8) we have that(3.12)or(3.13)but clearly this case is of no interest in our framework.

The interesting feature of the solutions in (3.11) is that holds in the special *finite* points and . This means that the values(3.14)are zeros of and by (3.7) it is easy to see that they are *double* zeros, i.e. .

In this situation, the domain wall solutions of equation (3.5), in analogy with the domain wall solution (3.2) of (3.1), are given by(3.15)(3.16)(The solutions here considered may be read also in *x* and the critical values are the ones defined in (3.14).)

Since for , solutions (3.15) and (3.16) have an oblique asymptote defined, respectively, by(3.17)in the denaturated zone they are similar to (3.2). Solutions (3.15) and (3.16) together with their asymptotes are reported in figure 1.

The main difference between (3.15)–(3.16) and (3.2) is in the front between the denaturated and non-denaturated zone of the molecule. In fact, for (3.2) we have that for and it is only for , and for (3.15) and (3.16) we have that holds exactly when and , respectively. Therefore, a clear sharp interface appears at the front of the domain wall.

We point out that the solutions in (3.15) and (3.16) are continuous with first derivative continuous, but second derivative discontinuous in . This means that they are weak solutions of equation (3.5). Indeed, the ordinary differential equation (3.5) does not satisfy the usual requirements of the classical existence and uniqueness theorem for the Cauchy problem. If we consider the initial-value problem defined by(3.18)we have that in each case two types of solutions are possible: and the one defined in (3.15) or (3.16), if the initial point is or , respectively.

The nonlinear coupling in (3.4) introduces a very sharp domain wall and this situation is not possible if a simple harmonic inter-base coupling potential is considered. The reason of the sharp interface is clear: the solution has to connect the two linear asymptotic branches ( and as defined in (3.17)) in a narrow zone.

### 3.2 The general nonlinear case

The extreme case of §3.1 is very special. For this reason, we have to turn back to the case where we have and we investigate what happens in the (more realistic) case of the nonlinear Klein–Gordon equation(3.19)In the static case, the first integral we obtain from (3.19) is given by(3.20)and by considering *C*=0 (because of the asymptotic requirements), from (3.20) we obtain the following real branch:(3.21)where .

We are not able to solve (3.21) explicitly by quadrature, but we can do a simple qualitative discussion to have the necessary information on the nature of the solutions. Let us introduce the function(3.22)and let us point out that *Y*=0 is a *double* root for (3.22). Then a classical theorem on analytical functions implies that: , such that . This remark is sufficient to show that the domain wall solution of (3.21) cannot have the strong *compactification* of the wall we have found in §3.1. In fact, if , from (3.21) we have that the improper integral,(3.23)does not converge and the value *Y*=0 is now only reached for . Since, for *ϵ*=0, we recover the extreme situation (3.8), our goal is to understand what happens when .

In figure 2, we report a numerical investigation of the behaviour of critical quantities associated with the domain wall solution as we approach the extreme case. In figure 2*a*, it is possible to appreciate the appearance of the *boundary layer* phenomena introduced by the nonlinear stacking term. To illustrate this fact, the derivative from (3.21) is shown. As , this derivative affords more and more higher jumps from the null value to the asymptotic value in a smaller and smaller zone. In figure 2*b*, we report the infinity norm of the difference between the derivatives obtained from (3.21) and from (3.7) as a function of *ϵ* and as . This is done in a log–log scale. From this plot, we note the existence of a threshold value , such that the convergence becomes faster for .

The conclusions of our analysis are that a nonlinear stacking interaction is responsible for a sharp localization of the front in the domain wall solution. Moreover, in the extreme case (and only in this case) of a purely nonlinear interaction we have a true compactification. In the general case, we have shown that as the domain wall solutions present a more and more pronounced boundary layer structure at the front between the denaturated and non-denaturated zone: a strongly curved elbow between the two asymptotic values of the wall appears.

The interesting fact is that the convergence of the generic derivative to the compact one as is not uniform. In our opinion, this phenomenon may suggest interesting experiments based on the observation of denaturation rate under mechanical tension. Goel *et al*. (2001) have analysed measures of the replication rate of DNA catalysed by a single enzyme moving along a stretched template strand. They have discovered the dependence on tension of the replication rate and they have presented a model for tuning the replication rate by mechanical tension. We may conjecture that the replication rate may be rather different in the above two ranges where the different convergence rates have been found.

## 4. Breathers and pulses

In the huge taxonomy of the models for DNA dynamics, the possibility that nonlinear effects might focus the vibrational energy of DNA into localized coherent structures is sometimes expressed by considering pulse waves, kinks or breathers.

By means of a small amplitude expansion in the original PB model, the classical NLS equation is retrieved. It is therefore natural to ask what is the effect of the intra-strand anharmonicity in this framework.

As usual, the NLS equation can be derived by considering a *multiple scale expansion* of equation (3.19). By setting (*ϵ* is a small parameter), we look for solutions in the form(4.1)where is determined by the linear dispersion relation and ‘*’ means the complex conjugate. By applying the cumbersome but standard procedure described in appendix B, we obtain in the new variables and , the equation(4.2)where *P* and *Q* are parameters and is the usual group velocity. The simplest localized solution of equation (4.2) is a breather given by(4.3)It is clear from (4.3) that the width of the breather depends on the parameter *P*, and therefore in the approximation used to derive the NLS equation, we cannot appreciate the role of the nonlinear intra-strand elasticity term. In fact, the usual multiscale expansion leading to the NLS equation rescales the term into the substrate potential and then the higher space derivative always appears as a linear term. (The same situation occurs when we consider the possibility of solitons in a nonlinear elastic medium as in Kivshar 1991.) Hence, the NLS equation determined from the PB model in the quasi-continuum limit does not give new insights into the problem.

Since the information obtained from (4.2) is not very interesting, our next step is to consider the possible existence of solitary waves with drop boundary conditions, i.e. pulse travelling waves, and to study the differences between a model with linear stacking and a model with nonlinear stacking. We shall investigate this fact by considering directly the continuous limit of the PB model in equation (3.19).

For this purpose, let us first consider the semilinear equation (3.1). By substituting in (3.1), the usual travelling wave ansatz , we obtain the ordinary differential equation(4.4)where a prime denotes differentiation with respect to . From (4.4), we obtain the first integral of energy as follows:(4.5)By using standard phase plane analysis, it is clear that pulse solitary waves are not possible. In fact, for a pulse existence, we need a double and a simple zero of the function(4.6)and this never happens in the case under scrutiny (see appendix B).

To overcome this problem, we have to once again remember that we are considering a mesoscopic model, where the DNA is highly deformed in order to assume the double helix configuration. For example, the molecule is bent and in the bending it is clear that the bases become closer, and therefore at least a change in the coupling constants has to be taken into account. Moreover, in such configuration, long-range interactions will become more and more important. For these reasons, recently, Archilla *et al*. (2001) have considered in detail the interplay between nonlinearity and geometry in DNA models. By considering a framework as simple as possible, they introduce a Hamiltonian where long-range dipole interactions are considered and the influence of bending may be analysed. In the continuum limit, it is not simple to consider the long-range effect and indeed a quasi-continuum theory is necessary. By oversimplifying the mathematics, it is possible to obtain some analytical results considering a dimer. This approach has been proposed in Larsen *et al*. (2004), where the effect of the bending (and also of the twisting) has been represented by a harmonic modification of the Morse potential. This means that we can modify equation (3.1) as(4.7)(We point out that for the new parameter *J*, a value of 0.05 has been proposed by Larsen *et al*. (2004).) Equation (4.7) has the energy integral(4.8)By considering drop boundary conditions, we must set *E*=0 and then have to consider the zeros of(4.9)It is easy to check that has a double zero in *F*=0 and another simple zero . This means that a pulse soliton is possible. To be more precise, solitary pulse waves exist if . Therefore, in this case, a pulse solitary wave may propagate.

These results are interesting because the formation of a localized coherent structure seems to need a geometrical deformation to modify the substrate interactions. This is exactly what many enzymes do in order to enhance the denaturation of the two strands (Calladine & Drew 1997).

As in §3, the next step is to understand what happens when anharmonicity of the interactions between bases on the same strand is taken into account. Therefore, we modify (3.19) as follows:(4.10)By introducing the ansatz , we obtain(4.11)When , the first integral of energy is given by (again imposing drop boundary conditions)(4.12)and, by introducing the notation , we have that(4.13)

First, we study the special case *ϵ*=0, where the speed of the solitary wave is exactly the linear wave speed:(4.14)In this special case, from (4.12) we get(4.15)and it is clear that the existence of a pulse solitary wave is possible only when .

It is possible to check that sufficient conditions for the existence of a pulse like *compact* solitary wave, as proposed in (Saccomandi 2004), hold for (4.15). This result implies that it is not only possible to find a pulse solitary wave, but also a *compacton*, i.e. a solitary wave with a compact support, exists.

To prove the existence of a solitary pulse wave in the general case , we must consider the zeros of(4.16)It is worth noting that, when (*ϵ*>0), we have to consider , while when (*ϵ*<0), we have to consider .

Of course, *F*=0 is a double zero for both functions , but, for the existence of solitary pulse waves, we have to ensure also the existence of another simple zero and to fulfil the conditions(4.17)for all .

If the case is considered, in order to satisfy (4.17), we must have . Under such condition, the existence of is ensured. In the case , it is possible to show that the first of the requirements in (4.17) is satisfied only when , but in this case never exists.

We may summarize these results as follows:

for , pulse waves exist for

*J*>1 when and for*J*<1 for ; for no pulse waves are possible;for , pulse waves exist only for

*J*<1 when (*ϵ*<0) and for (*ϵ*=0) there is a pulse compacton.

In both cases, the computation of the pulse solitary waves is possible via numerical quadrature applied to (4.13) and (4.15) for different values of *ϵ*. In figure 3, we show the shape of the various waves obtained by fixing for the cases: and ; for *ϵ*=0 (compacton solution); for , such that .

We point out that once again the case is a singular case and the corresponding ordinary differential equation does not satisfy the usual regularity conditions.

The results of this section show that a modification of the Morse substrate potential is necessary to predict pulse solitary waves and an anharmonic stacking potential is necessary to produce pulse compactons.

These results may suggest interesting questions. For example, it is well known that denaturation may be inhibited by constraining the DNA filaments, such that they are not free to bend and twist (e.g. enclosing them in carbon nanotubes). Therefore, our results in conjunction with single molecules experiments can be used as hint for understanding, at least qualitatively, the shape of the substrate potential in order to check if the Morse potential is realistic or not.

## 5. Other frameworks

Our results may be extended in several directions. The possibility of the existence of compactons is not restricted to the PB model. Almost all the models for the nonlinear dynamics of DNA are expressed in the continuum limit by means of a semilinear Klein–Gordon equation (or a system of such equations). These equations are semilinear because stacking interactions are usually considered to be harmonic for the reasons explained in §1.

The consideration of nonlinear stacking potential may introduce the possibility of compacton solutions for all the various models. We point out that our investigation does not depend on the special form of the intra-strand potential or the substrate potential. This is because a first integral of nonlinear Klein–Gordon model equations, given by(5.1)exists for any stacking potential and any substrate potential and this is not only for the choice of *W* in (3.3) and for the Morse substrate potential. The need of more refined substrate potentials in the study of the dynamics of hydrogen-bonded chains, as discussed in §6, is a well-known problem in the physics of condensed matter and it is not clear at all if the Morse potential is the most effective choice (e.g. Konwent *et al*. 1996).

To highlight possible generalizations of our method, in this section we consider another celebrated model for DNA dynamics: the Takeno & Homma (1983) and Yomosa (1983) planar base–rotor model. In this case, the soliton excitation is searched in the form of a kink solution. This allows to complete our discussion by considering another fundamental kind of solitary wave.

In the Yomosa model, we consider the angles of rotation and of the *n*th base in one strand and of the corresponding complementary base, and we take the Hamiltonian(5.2)where *I* is the moment of inertia of the base and, as in (2.1), *W* is the stacking interaction potential, whereas *V* is the hydrogen bonding energy. In their basic paper, Takeno & Homma (1983) consider the possibility of a general intra-strand potential, but they perform all the explicit computations only with a linear stacking interaction.

Yomosa (1983), by considering linear stacking and inter-strand potential given by , obtains, under the assumptions of the continuous limit approach, the two following differential equations:(5.3)Taking the difference of these equations, it is possible to derive, in the new dependent variable , directly the sine-Gordon equation(5.4)where and . On the other hand, if we sum the two equations (5.3), we can derive a linear equation for the phonon mode.

It is well known that for the classical sine-Gordon equation (5.4), it is possible to have a kink (+) or an anti-kink (−) solution in the form(5.5)where and . We point out that (5.5) is a solution of (5.4) also in the static case . In the following, we shall restrict the attention only to the static case.

If in the stacking potential we consider a nonlinear term, we have to modify equation (5.4) in the same spirit of the previous sections, and then we obtain(5.6)The effect of the anharmonic term is also clear in the static case. When , the first integral of (5.6), obtained by requiring that holds, is given by(5.7)It is easy to show that the *compactification* of the static kink holds. This result may be checked by taking into account that is an elliptic integral of the first kind, and therefore there exist two real points and , where and .

On the other hand, if in (5.6) we consider , then we have(5.8)and we obtain(5.9)From (5.9), it is clear that the solution approaches its asymptotic value in a sharper and sharper way as .

In figure 4, the solutions of (5.9) for several values of *ϵ* and the compact solution, for *ϵ*=0, of (5.7) obtained by numerical integration are reported. This plot shows that as , we have . For the aim of comparisons, the classical static kink (5.5) is also provided. In figure 5, we report the error as function of the parameter *ϵ*. It is easy to see that for the error tends to zero and there is a critical value , such that the linear (in the log–log scale) rate of convergence becomes superlinear for .

## 6. Concluding remarks

In our opinion, a detailed study of the complex deformation field of the DNA molecule is very difficult by using a model of masses linked by springs. Indeed, this is the reason why to investigate the complex configurations of the double helix, the elastica theory has been very successful. Inextensible rods are a practical way to synthesize the combined effects of the main bending, shearing and twisting deformations. The more complex theory of extensible rods has been an important step in coupling stretch with twist, and therefore in understanding new important phenomena as, for example, the overstretching transition (Marko 1997).

Of course, rod theory cannot be used easily in studying the denaturation phenomena, because in this framework the double helix is modelled as a single filament, and therefore the possibility of describing the breaking of base pairs associated with transverse deformations is lost. Only recently to overcome this problem, a new theory of the elastica, denoted as the theory of birods, has been proposed by Moakher & Maddocks (2005).

To overcome these difficulties, a multiscale modelling approach may be ideal. For example, Villa *et al*. (2004) propose an approach that combines a coarse-grained model of the DNA loop, based on the classical theory of elasticity, with an atom level model of proteins and protein–DNA interfaces based on molecular dynamics. This approach is still nowadays very heavy from a computational point of view and scarcely synthetic.

The advantages of our approach are the mathematical feasibility and the universality of our results. We have been able to display some exact solutions and, when this was not possible, the numerical treatment of the equations has been simplified and reduced to the numerical approximation of simple integrals. These computations were performed by a simple algorithm based on the Gaussian formulae and developed in the Matlab environment. Moreover, our results, as already pointed out, do not depend on the specific stacking energy but by the balance between nonlinear and dispersive effects. The use of the simple stacking potential (2.6) allows us to introduce the ratio *ϵ* between parameters of the nonlinear and linear terms in the constitutive equations as a direct and simple measure of *compactification*. In Destrade & Saccomandi (2006), a general mechanical theory explaining how to control the balance between nonlinearity and dispersion is presented. This theory may be used to generate more complex models.

We have already discussed why and in which context we think that our findings may be considered biologically relevant. On the other hand, it is clear that several important steps would be necessary to improve our results and to assess if our ideas may be useful in a biological setting. First of all, it is necessary to consider discrete, heterogenous and helicoidal models of DNA. To study the effect of anharmonic stacking potentials in such class of models could be possible only by applying suitable numerical methods. As we have pointed out in §3, the standard mathematical assumptions of regularity are no more valid, and therefore an appropriate choice of the numerical methods for the approximation of the models would be needed. In our opinion, the methods of Gorbach & Flach (2005) dealing with compact-like structures in discrete models are a good starting point in this direction. In such a way, we hope to have the possibility of using realistic constitutive parameters and to generate results more easily comparable with experimental results.

Another step in the pathway of our main philosophy is to consider non-convex stacking potentials. These kinds of energies are useful to model the coexistence of various phases of the molecule during the denaturation phenomenon. The biological relevance of the related models in the framework of DNA dynamics is the subject of a forthcoming paper.

The present work shows that there are strict connections between the different approaches to DNA modelling (molecular modelling, continuum mechanics, statistical physics, dynamical systems), and suggests that all the skills coming from these different fields of expertise are needed to improve our knowledge of the subject.

## Acknowledgments

This research was supported by the 2004–2005 PRIN funding scheme of Italian MIUR under the project Modelli Matematici per la Dinamica del DNA. G.S. is grateful to Claude Reiss for introducing him to the wonderful world of molecular biology. We thank all the anonymous referees for their constructive criticism that enabled us to improve the original version of our manuscript.

## Footnotes

- Received March 13, 2006.
- Accepted March 23, 2006.

- © 2006 The Royal Society