\magnification=1200 \font\tafont = cmbx10 scaled\magstep2 \font\tbfont = cmbx10 scaled\magstep1 \font\sixrm = cmr6 \font\petit = cmr9 \clubpenalty=10000 \widowpenalty=10000 \newcount\notenumber \def\clearnotenumber{\notenumber=0\relax} \def\fnote#1{\advance\notenumber by 1 \footnote{$^{\the\notenumber}$}{\petit #1}} \newcount\titbnumber \def\cleartitbnumber{\titbnumber=0\relax} \def\tbn{\global\advance\titbnumber by 1 \the\titbnumber } \def\titlea#1{{~ \vskip 1 truecm \tafont {#1} } \vskip 1 truecm \cleartitbnumber} \def\titleb#1{ \bigskip {\tbfont #1 } \medskip} \def\titletot#1{{~ \vskip 3 truecm \tafont {\centerline {#1}} } \vskip 3 truecm} \parindent=0pt \parskip=10pt \def\Z{{\bf Z}} %integers \def\R{{\bf R}} %reals \def\i{{(i)}} \def\unm{{1 \over 2}} \def\grad{\nabla} \def\Lapl{\triangle} \def\dd{{\rm d}} %"d dritto"; to be used in integrals \def\und{{1 \over \delta}} \def\pa{\partial} \def\dint{\int_{- \infty}^{+ \infty} } \def\arcsinh{{\mathop{\rm arcsinh}\nolimits}} \def\W{{\cal W}} \def\N{{\cal N}} \def\L{{\cal L}} \def\A{\AA} %Angstrom \def\er{~{\rm eV/rad^2}} \def\ea{~{\rm eV/\A^2}} \def\ec{~{\rm erg/cm^2}} \def\eb{~{\rm erg/rad^2}} \def\ps{{\rm ~ ps}} \def\smu{{\rm s}^{-1} } \def\erg{{\rm erg}} \def\a{\alpha} \def\b{\beta} \def\g{\gamma} \def\d{\delta} \def\eps{\varepsilon} \def\k{\kappa} \def\l{\lambda} \def\r{\rho} \def\s{\sigma} \def\th{\theta} \def\om{\omega} \def\Om{\Omega} \def\.#1{{\dot #1}} \def\~#1{{\widetilde #1}} \def\^#1{{\widehat #1}} \def\=#1{{\bar #1}} \def\[{\left[} \def\]{\right]} \def\({\left(} \def\){\right)} \def\RC#1{#1} {\nopagenumbers \font\tafont = cmbx10 scaled\magstep2 \def\titletot#1{{~ \vskip 3 truecm \tafont {\centerline {#1}} } \vskip 3 truecm} \parindent=0pt \titletot{SIMPLE MODELS OF NONLINEAR DNA DYNAMICS} %\vfill \centerline{Giuseppe Gaeta\footnote{$^a$}{Present address: Dipartimento di Fisica, Universita' di Roma, 00185 Roma (Italy)}\footnote{$^{*}$}{Supported in part by CNR grant 203.01.48} } \centerline{\it Centre de Physique Th\'eorique, Ecole Polytechnique} \centerline{\it 91128 Palaiseau (France)} \bigskip \centerline{Claude Reiss\footnote{$^{**}$}{Supported in part by CEE under contract SC1-CT91-0705} } \centerline{\it Institut Jacques Monod} \centerline{\it Tour 43, 2 Place Jussieu, 75251 Paris cedex 05 (France)} \bigskip \centerline{Michel Peyrard{$^{**}$} and Thierry Dauxois{$^{**}$} } \centerline{\it Laboratoire de Physique, Ecole Normale Sup\'erieure de Lyon} \centerline{\it 46 Allee d'Italie, 69364 Lyon cedex 07 (France)} \vfill \vfill \eject \parskip=0pt ~~~{\bf I - Introduction} \dotfill {\bf 1} \bigskip ~~~{\bf II - Structure and properties of DNA} \dotfill {\bf 5} \medskip II-1. Structure of DNA \dotfill 5 II-2. Degrees of freedom of the base pairs \dotfill 6 II-3. Degrees of freedom of the sugar-phosphate backbone \dotfill 7 II-4. Forces acting on the double helix \dotfill 8 II-5. Opening of the double helix \dotfill 9 II-6. Conformational defects in the double helix and their propagation: early proposals for solitons in DNA \dotfill 9 \bigskip ~~~{\bf III - Introduction of simple models for DNA dynamics} \dotfill {\bf 13} \medskip III-1. Classical models of DNA \dotfill 13 III-2. Nonlinear models of DNA dynamics \dotfill 14 III-3. Common features of DNA models \dotfill 17 \bigskip ~~~{\bf IV - Planar models of DNA dynamics} \dotfill {\bf 18} \medskip IV-1. A model of DNA torsion dynamics \dotfill 18 IV-2. The continuum limit of planar Yakushevich model \dotfill 20 IV-3. Solitons in planar Yakushevich model \dotfill 20 IV-4. A model of DNA vibration dynamics \dotfill 23 IV-5. The continuum limit of planar Peyrard-Bishop model \dotfill 24 \bigskip ~~~{\bf V - Helicoidal models of DNA dynamics} \dotfill {\bf 25} \medskip V-1. Helical geometry and helicoidal DNA models \dotfill 25 V-2. Introduction of helicoidal terms in the Peyrard-Bishop model \dotfill 26 V-3. Dispersion relations for the helicoidal Peyrard-Bishop model \dotfill 27 V-4. Nondispersive waves and breather modes in the helicoidal Peyrard-Bishop model \dotfill 28 V-5. Introduction of helicoidal terms in the Yakushevich model \dotfill 32 V-6. Solitons in the helicoidal Yakushevich model \dotfill 33 V-7. Exact solution for the helicoidal (1,0) soliton \dotfill 35 V-8. Numerical results for solitons in the helicoidal Yakushevich model \dotfill 36 \bigskip ~~~{\bf VI - Thermodynamics of the PB model} \dotfill {\bf 37} \medskip VI-1. Statistical mechanics of the PB model \dotfill 37 VI-2. Numerical simulations at constrained temperature \dotfill 38 VI-3. An improved model for denaturation \dotfill 40 \bigskip ~~~~~~~~{\bf Conclusions and discussion} \dotfill {\bf 42} \bigskip ~~~~~~~~{\bf References} \dotfill {\bf 44} \vfill \eject} \parskip=10pt \parindent=0pt \pageno=1 \titlea{I - INTRODUCTION.} \def\tbn{\global\advance\titbnumber by 1 I - \the\titbnumber . } DNA is central to all living beings, as it carries, for any species, all the information needed for birth, development, living and probably sets the average life-time of each species. Physically, DNA is a giant, double stranded linear macromolecule, with length ranging from millimeters (bacteria) to meters (humans) and even approaching kilometers (salamander). Yet DNA is always packed so as to fit a micrometer size space, either as a loose structure in procaryotes (bacteria), or in densely packed, dedicated unit, the nucleus, in eucaryotes (the more evoluted cells). As a rule, each cell carries a copy of the DNA molecule, and in multicellular organisms, the copy is identical in each cell. DNA is partitioned into genes (a few thousands to a few tens of thousands per copy), most of them coding for products needed for the cell life or its contribution to the life of the pluricellular organism. The formal genetic information (or coding sequence) of a gene is always deposited on one strand only of the double helix (but both strands carry genes). DNA may be considered a passive deposit of information, as expression of genetic information needs several steps carried out by dedicated machineries. Gene expression involves successively transcription - the gene is transcribed from the DNA by an enzyme, RNApolymerase (RNAP), onto a "messenger" RNA - and translation, the process by which the messenger is translated into the gene product (protein) by a special "factory", the ribosome. Gene expression is carefully controlled and regulated, to ensure supply of the gene product in amounts at the time and the place required by the actual circumstances (intra- or extracellular). It is perhaps worth describing transcription in more detail, as several steps of this process may depend on nonlinear dynamic effects, which have stimulated studies in this field. The process begins at a particular site of the gene, named promoter, to which RNAP sticks quite strongly. The essential partners (RNAP and promoter-bearing DNA) can be purified and the process carried out in the test-tube, in a simple salt buffer. Following complexation of RNAP with the promoter, the latter becomes partly unwound (1 to 2 helix turns disrupt their base-pairs to form a "transcription bubble"), thus exposing the sequence serving as template for transcription [Hip]. The unwinding, which involves $50-100$ KCal/mol, occurs in the test-tube in the absence of any supply of chemical energy, and may therefore feed on the thermal reservoir by some non-linear mechanism. Once the RNAP has complexed the promoter, productive transcription ("elongation") can start. The messenger RNA is actually synthesized from its monomeric constituents (which are X-triphosphates, with X=A,U,G,C, which are cleaved to monophosphates), in a sequence complementary to the DNA template strand. Interestingly, the transcription bubble, created at the promoter, is conserved and accompanies the transcribing RNAP [Wan1]. The latter is prone to frequent pausing at well-specified places of the template. One of the (many) signals responsible for pausing is borne by a row of six A's, [Wan1] which is known to result in unusual local DNA conformation (permanent bending). In the complex pausing under control of that signal, no transcription bubble could be detected, but the latter is restored as transcription resumes. For DNA, transcription appears as a succession of an initial production of a local defect (bubble) at the promoter, translocation of this defect over some distance, arrest and disappearance, then restauration and again translocation of the defect, etc... Taken together, the results suggest that the transcription bubble, created at the promoter (possibly by a nonlinear process, see above), moves along the template with RNAP and vanishes as the latter pauses. RNAP can resume transcription as the bubble is restored. The famous discovery of the double helix had emphasized a strong relationship between {\it structure} and function in molecular biology. Transcription shows that, in fact, DNA is a very {\it dynamical} entity, undergoing very large deformations and should not be viewed merely as a solid with a particular structure. Another important motion of the DNA molecule is its ``breathing'' or fluctuational opening. In these very large fluctuations, base-pairs are temporarily broken and the two bases are exposed for chemical reaction for a very short time ($10^{-7}$ s). These fluctuational openings can be considered as intrinsic precursors for the denaturation (see below) and they could play a role in carcinogenesis by external molecules [Sob2]. The molecular deformation involved in melting or in fluctuational openings are so large that they cannot be described by linear approximations. Biomolecular dynamics is a fascinating topic for nonlinear science because it is related to basic phenomena of life and we {\it know} that it has to be fundamentally nonlinear. The goal of this paper is precisely to motivate, describe and discuss some of the nonlinear dynamical models proposed so far to study the dynamical processes taking place in DNA, and specifically the processes mentioned above, namely the fluctuational opening, and the motion of the transcription bubble. The interest in nonlinear dynamics of DNA was indeed started, as it will be also recalled later when shortly reviewing relevant literature, by the proposal of Englander et al. [Eng], that the RNAP could take profit from the existence of solitons propagating along the double helix, and "surf" them in order to perform transcription without having to provide energy to open the hydrogen bonds of the basis pairs which have to be transcribed. Such an hypothesis would also imply that it is crucial for RNAP to perform transcription at essentially constant speed (at least between two stops of the process); it has also been speculated that one of the functions of the redundance in the genetic code - i.e. the existence of different triplets of bases coding for the same genetic information - could indeed be to provide flexibility to the sequence in order to allow for constant speed of the transcription bubble [Wan2]. While the first proposal by Englander et al. was essentially qualitative, several models have been proposed in order to substantiate this idea in quantitative terms. While here we discuss in details some recently proposed models of the sort and shortly review the existing literature on the subject, we hope the present paper could serve a more general purpose, namely that of attracting a greater number of physicists - not only biophysicists, but also theoretical and mathematical ones, especially those working in Nonlinear Physics - to this fascinating field. \medskip The plan of the paper is as follows: in chapter II we briefly recall the structure of DNA, with particular emphasis on the determination of the relevant degrees of freedom (at body temperature) for its dynamics, and in particular for the processes we want to study; we also give a short description of these processes using the more precise nomenclature that we introduce in this section. In chapter III we discuss the general structure and common features of the models used to study DNA nonlinear dynamics; we also give a very short overview of the existing literature, both to trace the history of the subject and to allow the interested reader to access works using approaches slightly different from the one we will consider in the following. Chapter IV is dedicated to study of so-called {\it planar models} of DNA dynamics, i.e. models considering only stacking and intrapair base interactions (in other words, only nearest neighbours interactions along the double helix); in particular, we will consider the model proposed by Yakushevich [Yak2] for the torsion dynamics of DNA, and aiming to describe the motion of the transcription bubble, and that proposed by Peyrard and Bishop [Pey1] for the opening - or vibration - dynamics of DNA, and aimed at describing fluctuational openings; in both cases we will consider the discrete model and its continuum limit, and in the case of Yakushevich model we also discuss the soliton solutions. In chapter V we consider so-called {\it helicoidal models} Êof DNA dynamics: these descriptions consider also the interaction between bases which, although being distant along the double helix, are nearby in three-dimensional space due to the helical geometry of DNA; in real DNA this interaction is indeed present and mediated by water filaments connecting the bases; the introduction of this interaction yields {\it qualitative} changes (beside, obviously, quantitative ones) in the dynamics of the model, already at the level of dispersion relations. In particular, we will consider the helicoidal versions of the two models considered in chapter IV, and extract predictions on physically observable quantities from these; it will turn out that these predictions are in good agreement with experiments, both for the breather solutions in the Peyrard-Bishop model and for the soliton solutions in the Yakushevich model (although in the latter case the introduction of helicoidal terms have smaller consequences). Chapter VI is different from the previous ones, in that it is concerned with the {\it thermodynamics} of nonlinear DNA models, and specifically of the Peyrard-Bishop model; this is quite relevant, since after proving that these models correctly predict the existence of breather solutions, it has to be investigated that they also give a correct statistical mechanics description of fluctuational opening; no such results are available for the Yakushevich model, and actually the results that we present here for the Peyrard-Bishop model are temporary ones, as the relative computations are continuing and more complete results will be available in the future; we think it is nevertheless of interest to give also a discussion of such a relevant aspect of the question in its current state. \medskip Throughout the paper, equations are numbered independently for each section, and section for each chapter, while for figures, tables and footnotes the numerotation runs through the whole paper. The quoted reference are referred to by the first three letters of the name of the first author and, when necessary, a number; the references are listed in alphabetical order at the end of the paper. \bigskip The work of G.Gaeta has been supported in part by C.N.R. under grant 203.01.48; the work of C.Reiss, M.Peyrard and T.Dauxois has been supported in part by EEC under contract SC1-CT91-0705. \vfill \eject \titlea{II - STRUCTURE AND PROPERTIES OF DNA} \def\tbn{\global\advance\titbnumber by 1 II - \the\titbnumber . } In order to understand the dynamics of DNA, the first step is to consider its structure and degrees of freedom. DNA polymorphism results mainly from the interplay between the conformational flexibilities of the sugar moiety and the degrees of freedom corresponding to the base-pairs and their mutual positionings. Let us consider these various degrees of freedom. \titleb{\tbn Structure of DNA} DNA is assembled from two large linear macromolecules (see fig.2) sharing a similar chemical structure built from three basis elements: a phosphate- glycosidic bond backbone (fig.1a); a sugar moiety, i.e. the furanose- derivative deoxyribose (fig.1b) with its -C3'-C4'-C5'- segment being part of the backbone; four bases, i.e. adenine (A), guanine (G), cytidine (C) and thymine (T), linked to the C1' of the sugar ring (fig.1c). The linear macromolecule (or "strand") assembled from these three components is characterized by its polarity (there is a 3' end and a 5' end), and the polarity- specified sequence of the bases borne by consecutive deoxyriboses, which is the carrier of the genetic information. Two strands do associate to form DNA (fig.1) according to two "rules": $i)$ the strands run parallel to each other but have opposite polarities; $ii)$ the link between the strands is made by pairing the bases they bear, where as a rule only A-T and G-C pairs can form; these pairs are held together by hydrogen bonds (two for A-T pairing, three for G-C pairing, see fig.1). The double strands wind into the well known Watson- Crick double helix. In solution, the helix is fairly rigid and can be considered rod-like for sequences from tens to a few hundred bases long. Usually the helix is found to be right handed, i.e. looking down the helix axis in either direction, each strand winds clockwise as it moves away from the observer; but it can be left handed under particular environmental circumstances. The helix fits into a cylinder about $20 \AA$ in diameter. A helix has 10 to 11 base pairs per turn with a pitch averaging $32 - 34 \AA$; i.e., consecutive base pairs are some $3.2$ to $3.4 \AA$ apart along the axis of double helix. The poor precision of these data reflects the fact that the DNA double helix is highly polymorphic. This is true not only of DNAs of different sequences, but even for a DNA of given sequence when subjected to different environmental parameters, like hydration or ionic strength. Clearly, the molecule has several, closely spaced, energy minima which become occupied upon minor changes of these parameters. Data from X-ray diffraction, collected from a variety of single oligomeric DNA crystals, show indeed that the precise structural parameters of a double helix depend on its sequence; and vibrational spectroscopy analysis, performed on the same DNA in solution or in crystalline form, show that the double helix parameters can vary significantly from one phase to the other. An in-depth survey of DNA structure can be found in [Sae]. \bigskip \titleb{\tbn Degrees of freedom of the base pairs.} \smallskip They are best visualized in the right- handed orthogonal axial set $Oxyz$ (see fig.2), where $O$ is chosen at the "center" (close to N1 of the purine) of the pair under consideration. The base-pair is assumed for the moment to be planar and perpendicular to the helix axis. $Oz$ is taken along this axis, $Oy$ runs from C6 (pyrimidine) to C8 (purine), so that $Ox$ intercepts the H-bonds of the pair (fig.2). First, the flexibility of the H-bonds enables rotational freedom between the bases of the pair (fig.2), which are thus not necessarily coplanar; "opening", "propeller- twist" and "buckle" correspond to rotations around $Ox$, $Oy$ and $Oz$, respectively. They are measured by the dihedral angles between the planes of the individual bases (looking down the rotational axis as indicated on fig.2, the angle is positive if the nearest base is rotated clockwise relative to the farthest one). The bases of the pair also have translational freedom (stagger, stretch and shear, see fig.2). In addition to the intra-pair degrees of freedom, the pair as a whole has also rotational as well as translational degrees of freedom, with respect to the next pair. First, the mean plane of the base-pair can depart from the position perpendicular to $Oz$ which has been assumed above. Rotation around $Oy$, $Ox$ and $Oz$ have been dubbed "twist", "roll" and "tilt", respectively. These are measured by the dihedral angles made by the mean plane of the base- pair with the plane perpendicular to the helix axis. Translational freedom along the three axis exist and have been named "rise", "slide" and "shift". Of course, the bases and/or base-pairs are not free to rotate or translate according to these degrees of freedom, as these movements may be opposed by hindrances. Two hindrances are predominant. The first results from the presence of nearest neighbour base-pairs, on both sides of the pair under consideration. If the bases are sketched as rectangular blocks (fig.3), it can bee seen that the propeller twist for instance can give rise to steric clashes with the next base-pairs, as could buckle, roll, tilt, twist, etc... Conversely however, the degrees of freedom of each neighbouring pair can be instrumental in removing these steric hindrances, and so on for their next neighbours (see legend to fig.3). Thus the clash can be handled and eased with the mutual help of a row of base-pairs, each modifying its position in accordance to its degrees of freedom, following rules first proposed by Calladine [Cal1]. The second hindrance is imposed by the glycosil bond linking the base to the sugar. Although this bond is covalent hence very strong, the glycosil ring to which it is attached has a rather flexible structure, as we will see below. Hence here again, the hindrance is rather soft and can be accommodated within rather large limits, but in this case mainly by intranucleotide rearrangements. \titleb{\tbn Degrees of freedom of the sugar-phosphate backbone} Six torsional angles characterize the backbone structure between two consecutive phosphorus atoms (see fig.4), i. e. in the 3' to 5' direction. In addition, the conformation of the sugar is characterized by four more torsional angles ($\nu_0$ to $\nu_4$, $\nu_3$ identical to $\delta$). The last important structural parameter is the torsional angle of the glycosidic bond, $\chi$. The values of torsional angles are in general restricted to the common sterical ranges, {\it syn} ($0^0$), {\it anti} ($180^0$), {\it synclinal} ($\pm 60^0$) and {\it anticlinal} ($\pm 120^0$) (fig.4). In DNA, not all these ranges may be accessible; for instance the glycosil torsional angle $\chi$ is restricted to syn and anti ranges mainly (see figs.4a and 4b). Furthermore, most of these torsional angles are correlated (only $\gamma$ is independent). It is instructive to visualize the positions of the base and of the 5' atom relative to the sugar, observed in common DNA conformations (figs.4a to 4e). The five-membered furanose ring is usually not planar but "puckered", either in the "envelope" form (4 atoms of the ring approximately coplanar, the fifth out of plane by less than about $0.5 \AA$), or in the "twisted" form (two adjacent atoms are displaced on opposite sites of the plane defined by the three other atoms; those displaced on the side of C5' are called "endo", those on the opposite side are "exo"; nonsymmetrical twist is characterized by its major departing element (see fig.5). The conformation of the sugar is concisely described by the pseudo-rotation phase angle $P$, where $\tan P = ( (\nu_4 + \nu_1 )-(\nu_3 + \nu_0 ))/(3.08 \nu_2))$, which allows to compute the five torsional angles $\nu_j$, $j=0,...,4$, by $\nu_j = (\nu_0 \cos (P + j 144^o ))/(\cos P)$. The case $P=0^o$ corresponds to $\nu_2$ maximally positive, that is a symmetrically C2'exo-C3'endo twisted conformation, $P=180^o$ to the mirror image of the latter (C2' endo-C3' exo). Structural analysis of nucleic acids indicate that two $P$ ranges are strongly favoured, $-1^o < P < 34^o$ (C3' endo) and $137^o < P < 194^o$ (C2' endo), preferred in DNA. Notice the consequences of puckering on the relative directions of the C1'-Base (glycosyl) bond and the C4'-C5' backbone bond, and on the relative positions of consecutive phosphorus atoms (figs.5a and 5b). The variation of the energy of furanose (in nucleosides) with $P$ is displayed on fig.6. The preferred C3' endo and C2' endo are separated by a barrier of about 1.5 Kcal and the highest barrier found over the whole $P$ range is about 5.5 Kcal. Transitions between the various conformations is therefore easy: the sugar ring moiety of nucleic acids is indeed highly flexible, as already noticed. This flexibility can conciliate (to some degree) constraints imposed to the glycosil bond by the base (or base pair) and constraints imposed by the backbone on the C5'-C4'-C3' link, with minimal expenditure of energy. \bigskip \titleb{\tbn Forces acting on the double helix} \smallskip We distinguish schematically inter and intra base pairs, and inter and intra backbone forces. In the mean plane of the base pair, protons are exchanged between the NH donor groups of one base and N or O acceptors of the other. Although these hydrogen bonds are weak and not highly directional [Vin], they contribute to the stability of the Watson- Crick type pairing, hence have a crucial role in coding the genetic information, its transcription and replication. Notice that an alternative to the Watson-Crick pairing is the Hoogsteen pairing: in the former, the H-bonds involve atoms or groups born by the six-membered rings of the purines only, whereas in the latter, N7 of the five-membered ring can be an acceptor (see fig.2 for numbering). Intra-base paired H-bonds are easily disrupted at physiological temperatures by a variety of chemical agents and physical parameters at concentrations and values commonly encountered in living systems (see below). We mentioned earlier the hindrances in the degrees of freedom, associated with steric clashes between neighbouring base pairs. In addition, base pairs, stacked vertically along the helix axis, interact mainly via p-p electron coupling, dipole-dipole interaction, London dispersion forces and hydrophobic forces (in solution). These forces result in a complex interaction pattern between overlapping base pairs, with minimum energy distance close to $3.4 $ \AA\ in the normal DNA double helix, strong repulsion between both pairs below this distance, and fair attraction above this distance. The interaction depends strongly on the degree of overlap between two consecutive pairs (stacking), which varies greatly according to the precise values of the geometric parameters associated with the degrees of freedom of each base or base pair. Like H-bonds, base pair stacking depends on temperature, the state of protonation of the bases, the local dielectric constant and other parameters external to the nucleic acid, summarized as "environmental" parameters. Long range intra- and inter- backbone forces depend mainly on the presence of the phosphate groups. The distances between phosphates on the two strands is ca. $20 $\AA\ at least, hence their interactions are weak. In contrast, along the same strand, the distance between phosphates can be as small as $5 $\AA , meaning that their mutual repulsion can be rather strong. To remain in its double helical, native form, DNA must be kept in a media having a minimal ionic strength. The phosphate groups are then shielded by the counterions supplied by the media. The shielding is very stable; as the NaCl concentration of the media is changed from $0.5$ mM to $0.5$ M, the number of Na$^+$ ions shielding the phosphates remains constant, since on the average $0.88$ Na$^+$ shield each phosphate group of the backbones, throughout this ion concentration range. Again, environmental parameters in the physiological range can alter the shielding (type of counterion element and valency, pH) and structural transitions of the double helix can significantly modify the interphosphate distances (see fig.5e and 5d). \titleb{\tbn Opening of the double helix} The double helix is not only polymorphic, it can also undergo partial or even total strand separation; the double strand is also referred as the "native" state of DNA, and the process of separation is also called "unwinding", "opening", "melting", or "denaturation". DNA unwinding can be obtained in the test tube, adjusting various physical or chemical parameters, such as the temperature or the ionic strength of the buffer. The opening, which corresponds to the disruption of the H bonds joining the bases in pairs located in the unwound region, can be monitored experimentally. A convenient parameter for this is the DNA absorbance in the near UV (more precisely between $240$ and $280 nm$), which increases by some $40 \% $ as the double helix melts, as a consequence of the disruption of stacking interaction between consecutive base pairs. In figure 7 it is shown the "melting profile" of DNA, i.e. the variation of absorbance at $260 nm$ as the temperature rises from $50$ to $80 {}^o C$ [Gab1]. There are two striking features: $i)$ melting is cooperative, as $80 \%$ of the increase in absorbance occurs in an interval of $ 5 {}^o C$; $ii)$ melting occurs in several discrete steps (peaks). It was shown that each peak corresponds to the sharp melting, or highly cooperative opening, of a well defined segment of DNA, which can be mapped precisely within the double helix; its size is roughly proportional to the peak area. The mean peak temperature characterizes the mean AT/GC composition ratio of the segment\fnote{A linear relationship holds between the melting temperature and the AT/GC composition, with pure AT segments melting well below the pure GC ones; recall the latter exchange three H bonds per base pair, the former only two.}. The width of the peak characterizes the dispersion around the mean composition of the segment. The experiment by which fig.7 was obtained, was performed in a buffer of specified ionic strength; dividing (respectively, multiplying) this by a factor of ten, the whole melting profile would result shifted of about $18 {}^o C$ towards lower (respectively, higher) temperatures [Gab2]; the profile itself would be modified to some extent, as cooperativity restricts to smaller (respectively, extends to larger) segments as ionic strength is reduced (respectively, increased). The native/denaturated phase diagram can be computed from the sequence [Mar1]. \titleb{\tbn Conformational defects in the double helix and their propagation: early proposals for solitons in DNA} Given the numerous, coupled degrees of freedom found in DNA, a full description of its molecular motions is difficult. For modelling, one has therefore to select those degrees most likely to dominate conformational changes. Models must furthermore cope with the very existence of non-linear couplings between molecular units of DNA, or of nonlinear interaction potentials within this molecule, such as Sine-Gordon or multistable potentials (illustrated for instance on fig.6) It was indeed suggested by model building and computation that the change of sugar pucker from C2' endo to C3' endo produces a kink of the helix axis of some $40^o$ and unwinding of the adjacent base pairs by $-10^o$. The change could be the result of a solvent molecule hitting locally the double helix (brownian motion). It was further suggested that, provided the collision occurs with the right momentum, the C2' endo $\to$ C3' endo transition, once created, could propagate as a harmonic wave and travel some distance along the helix as an acoustic wave. It was also found that transient, large amplitude, localized openings of a few base pairs occur in the molecule due to fluctuations of the structure. This DNA "breathing" is thought to be responsible for the exchange in $D_2 O$ (heavy water) of H atoms of the bases normally not accessible to the solvent, and intercalation of drugs with shape and chemical structure resembling those of the bases [Ban, Eng, Lew, Sob1, Sob2, Yom1]. Both effects lead to the proposal [Ban, Sob1, Sob2, Eng, Lew, Yom1] that the transient opening of the double helix and its propagation could be handled formally as a soliton. In particular, the isotopic exchange experiment was analyzed by analogy to solitary waves occurring in a linked (interbase H-bond) pair of chains, carrying coupled (sugar-phosphate backbone) pendula (the bases). In this mechanical analog of the double helix, solitary waves build up due to non-linear coupling between potential energy (gravitational) and restoring torque energy. This proposal arose strong interest among non-linear physicists and stimulated a wealth of theoretical studies, considering the possibility of solitary waves building up and propagating in DNA [Tak1, Kru3, Sco3, Zha, Xia, Tak2, Yak1, Bal, Mut2]. To most molecular biologists however, the presence of solitons in naked DNA deemed a laboratory curiosity, as they know that genetic information is derived (replication, transcription) from DNA by dedicated enzymatic machineries: if solitons are to participate in the genetic processes, they must occur in these DNA- specific protein complexes, which was so far not documented by experimental evidence. To date, detailed structural and dynamic data exist for such complexes. They show that the genetic information is invariably derived from single-stranded DNA, meaning that the double helix must be locally unwound beforehand. The single-stranded stretches are of limited size and move in close synchrony with the enzyme, as it performs its function along the nucleic acid. Remarkably, this movement is not uniform but occurs in fast steps of variable lengthes, followed by pausing of variable durations: the genetic information is derived from the DNA matrix in steps of variable size. Let us describe in more detail these features in the case of DNA transcription. In this process, the genetic information carried by one strand is transcribed by the enzyme RNA polymerase (RNAP) into a single- stranded nucleic acid, the messenger RNA. Transcription can be carried out in vitro by the purified enzyme (extracted from a bacteria for instance), the DNA template and the four nucleotide building blocks, and has been studied in great detail. The actual transcription process is preceded by several preliminary steps; first, RNAP binds to DNA through electrostatic binding to the phosphates. This binding is not static, as RNAP slides along DNA at high speed (on the mean, the displacement amounts to 1000 bases per second). It has been suggested that RNAP displaces at that rate by "surfing" a solitary "kink" wave. Upon meeting a particular site on the matrix, the promoter, located slightly ahead of the transcription starting site, RNAP stops, then complexes to the promoter . Following a rather lengthy activation process, the complex is then ready to start transcription. This complex remains stable and active for long times (equilibrium constant of the order of $10^{10} M^{-1}$). The activation of the complex is a complicated process, during which both partners undergo transconformations. The most dramatic one is experienced by the promoter, which is locally unwound and melted (separation of the strands following breaking of the base pairs) over 15-20 consecutive base pairs. Furthermore, the double helix axis, which was straight before complex formation, is now bent at some point within the melted region. Strikingly, this melting, which involves the expense of $50$ Kcal at least, feeds only on the thermal reservoir and in particular does not require chemical energy. It is plausible that the melted region results from the local accumulation of "breather" solitons. Using various biochemical methods, the structure of the complex can be described at the molecular level. RNAP "footprints" the promoter (DNA section covered by RNAP) over 50-60 base-pairs, including the transcription bubble [Hip]. The conformation of the footprinted, still helical DNA, is deeply modified, as the width of the helical grooves is changed [Duv1] and the helix axis, which was straight before complexation, is now markedly bent [Heu]. Upon complexation, about half of the phosphate-shielding counterions are removed from the footprinted DNA [Sch]. The kinetics of the complexation are characterized by four intermediate states at least, leading from the initial promoter binding of RNAP to the final complex, ready to start transcription [Duv2]. The process is strongly temperature dependent; for instance, the transcription bubble forms sharply in a very small ($1-2 ~^o C$) temperature range. Comparison of many promoter sequences, recognized by RNAP from a given species, showed that the promoters share a common pattern of stability, encoded by a succession of sequences constraint in weak (AT) or strong (GC) base pairs [Jac1]. Synthetic sequences, constructed according to these observations but sharing no homology with natural promoter sequences, proved quite efficient promoters in vivo [Jac2]. The experiment enabled further to identify some of the factors (base composition constraints) controlling the transcription starting site (address and multiplicity), transcription efficiency and some elements involved in the selection of the transcribed strand [Jac1]. It led to a proposal of a mechanism of promoter-RNAP complex formation, based on the formation of specific hydrophobic contacts between the complex partners in both grooves of the promoter DNA. It was also demonstrated that in vivo the efficiency of a given promoter can be tuned over an almost ten-fold range, by modifying the supercoiling of the template [Jac3]. Since the precise helix geometry of DNA is modified by controlled changes in supercoiling (positive or negative), the hydrophobic contacts are modified as well, explaining consistently the change in promoter efficiency. The supercoiling can be introduced by various means at a large distance (thousands of bases) from the promoter, but the precise mechanism of propagation of supercoiling along the template is not yet clear (soliton-like propagation of the topological defect, or equipartition of supercoiling over all units of the template?). Initiation of transcription requires selection of the transcribed DNA strand and of the initiation site, which is located inside this melted region. As transcription proceeds, this melted region (and the associated helix bending) "moves" with the transcription complex until transcription stops. Actually, upon transcription elongation, the base- pairs ahead of this region melt, those behind reform. With respect to the normal helix winding, the double helix is overwound upstream and underwound downstream the melted region. The melted region appears thus as a conformational defect which, once formed, propagates with the progressing transcription complex. During propagation, the size of the region remains probably constant. It is not known whether propagation dissipates energy, since chemical energy is available from the nucleotide tri- to mono-phosphate clivage which takes place at each transcription elongation step. As already mentioned, transcription elongation is not a steady process, but proceeds in a series of fast steps (up to 1000 bases transcribed per second) separates by extensive pausing at well defined places. Pause sites are from a few to a few hundreds bases apart, and the duration of the pause can exceed 10000 fold the time required for a single elongation step at non-pausing sites; over three- fourth of the total elongation time is spend pausing. The melted region pauses and progresses just as does the transcription complex. Several types of signals are deposited on the double helix and force transcription to pause. One of them - the first and only one so far - has been recently identified [Wan1]; it is associated with a conformational singularity of the double helix, introduced by a series of six adenines in a row which is thought to induce a permanent bent in the helix axis. If we assume that RNAP "surfs" a topological soliton associated with the melted region and the bent helix axis, it could be that its encounter with the localized, permanent bent associated with the row of 6 A's annihilates both bents and stops temporarily its progression, until thermal fluctuations for instance have restored proper bending of the helix axis. \vfill \eject \titlea{III - INTRODUCTION OF SIMPLE MODELS FOR DNA DYNAMICS} \def\tbn{\global\advance\titbnumber by 1 III - \the\titbnumber . } \titleb{\tbn Classical models of DNA} Transcription displays several features (localization of the promoter, building up of melted region and helix bending, propagation and pausing of the latter) which could perhaps be better understood within the frame of a non-linear model. DNA replication, though less well explored yet, could deserve a similar approach. However these in-vivo processes are extremely complicated. Moreover they are non equilibrium processes which makes the thermodynamics approach even more complicated. This is why physicists have concentrated their attention on phenomena which looked more amenable to analysis, with the hope that, developing some models for these phenomena would be useful in the long term for understanding the biological processes. Certainly, one of the reasons for the attraction exerted by DNA dynamics on physicists is the fact that DNA is the fundamental molecule of life. A second reason is that the dynamics of DNA transcription or replication involves extremely large motions which go well beyond any linear approach. Therefore, DNA provides a beautiful example of a system in which {\it nonlinear} effects are essential, and it has undoubtedly served as a ``toy model'' for fundamental studies in nonlinear dynamics. The problems which have been considered so far have been selected because they are well defined. In these approaches, the structure of DNA is supposed to be perfectly regular, although some calculations included inhomogeneities to represent the genetic code, and the environment of the molecule is ignored. More precisely, it is simply indirectly introduced through its effects on the model parameters, or through thermal fluctuations in the calculations involving temperature. Two types of questions have been asked: {\it i)} what are the possible nonlinear excitations of the molecule? In these approaches the models are designed to have solutions which can mimic some of the aspects of DNA dynamics known from experiments, and in particular the formation of local conformation changes which can propagate along the molecule as in the transcription, or very large base motions that would represent the fluctuational openings. {\it ii)} how do DNA models behave at non-zero temperature, and, in particular, can they describe the thermal denaturation of the molecule observed experimentally? This is also a problem which is well defined in terms of physics because the thermal denaturation appears as a kind of ``melting'' phase transition. Since denaturation has been extensively investigated experimentally when it was considered as a possible way to reach information on the sequence, the results of the theoretical analysis of the thermal denaturation can be compared to experience, thus providing some means to refine and calibrate the DNA models. Nonlinear dynamics of biological molecules was first considered for proteins in the famous Davydov model of the $\alpha-$helix [Dav1,Dav2] which couples a deformation of the protein molecule with an internal vibration of a C=O bond. Nonlinearity acts in this case to trap the energy of the vibration into a localized state that can propagate along the molecule. A proper description of this phenomenon can only be achieved within a quantum theory. The situation in DNA is different since we are interested in large conformational changes of the molecule, occurring at low frequencies and involving large groups of atoms in the molecule. Therefore a classical description could suffice. The difficulty however is to select the degrees of freedom which are relevant and must be involved in the model, and to give appropriate effective potentials for them. \titleb{\tbn Nonlinear models of DNA dynamics} The pioneering work pointing out the possible role of nonlinear excitations of soliton type in DNA is due to Englander {\it et al.} [Eng]. Since solitons are extremely stable localized excitations in nonlinear systems [Sco1, Sco2], they are good candidates to represent local openings in DNA. In their analysis, Englander {\it et al.} assumed the existence of states which stay permanently open and can move by a random diffusive motion. They compared the characteristic frequency of the openings, which can be estimated experimentally, to the typical size of the open states and concluded that a reasonable size would be about 10 base-pairs. A simple model, analogous to the sine-Gordon model for a chain of coupled pendulums, was used to estimate the energy of the open states. The interest of this work goes well beyond the first description of the open states it proposed, because it initiated the study of nonlinear dynamics of DNA. The various models that appeared then can be classified according to the properties they attempt to describe and the degree of accuracy they use in the representation of the molecule. A very detailed structural study of a locally deformed state of DNA assumed to be a precursor state of denaturation was proposed by Sobell [Sob3]. However this approach is simply a qualitative description of this ``premelton'' which does not explicitly consider the nonlinear effects that could be associated to the large distorsion of the double helix. It contains nevertheless many of the ideas found in subsequent works, and in particular the idea that the premelton can move along the molecule while preserving its shape. Another study which considered the structure of DNA in rather great details is the work developed by Krumhansl and Alexander [Kru2] to describe the $B$ to $A$ and $B$ to $Z$ DNA conformational changes. The model couples a displacement of the nucleosides to the pucker of the sugar ring, which is in turn coupled to the twist of the bases. Additional coupling due to the hydration spine of the molecule were introduced in the hamiltonian, but neglected in deriving equations of motion. In spite of the large number of degrees of freedom per base-pair, using some simplifying assumptions, and in particular the continuum limit approximation which assumes that the variables change only slowly from one site to the next, the authors were able to derive analytical solutions describing the region which interpolates between two conformations of DNA. Although these approaches have the interest of providing a description of the molecule which is sufficiently close to the real structure, they yield equations which are in general not tractable analytically. Therefore, most of the analyses of nonlinear excitations in DNA rely on much simpler models. The case of longitudinal solitonlike equations has been treated first with the continuum model of a rigid rod [Mut2] and then with a Toda model [Mut4]. This work concluded that, at physiological temperature, a significant number of longitudinal solitons can be generated. Although these approaches stress the importance of nonlinearity in DNA dynamics, their conclusion are difficult to probe experimentally because the measurements on longitudinal motions of DNA can hardly go beyond the determination of the speed of sound along the molecule. On the contrary transverse motions are easier to observe because they are related to the denaturation of the base-pairs which can be detected by a variation of the UV absorbance [Sae] or exchange of the imino proton [Gue]. Moreover they are directly related to the biological activity of the molecule. Therefore they have been the subject of a larger number of studies. Yomosa [Yom] proposed a plane base-rotator model in which the bases are assumed to rotate within a plane which is orthogonal to the axis of the double helix. Denoting by $\chi_n (t) $ and ${\chi'}_n (t)$ the angular positions of the two bases of the $n$th base-pair, a dipolar coupling is assumed between adjacent bases while the interaction between the two bases in a pair is written in the general form of a double Fourier serie which is then truncated. Thus the hamiltonian of the model is $$ \eqalign{ H = \sum_n & \biggl\{ {1 \over 2} [\dot \chi_n^2 + \dot {\chi'}_n^2] + A ( 1 - \cos \chi_n) + A (1- \cos {\chi'}_n) \biggr. \cr & \biggl. + B ( 1 - \cos \chi_n \cos \chi'_n) + S [ 1 - \cos (\chi_n - \chi_{n-1})] + S [ 1 - \cos ({\chi'}_n - {\chi'}_{n-1})] \biggr\} \; \cr }$$ where $A$, $B$, $S$, are potential parameters. Then a continuum limit approximation is used together with an expansion of the dipolar coupling terms which are reduced to an harmonic interaction. As a result, considering separately in-phase or out-of-phase motions of the bases, the equations of motion can be reduced to sine-Gordon or double sine-Gordon equations [Yom2]. The model has been later refined by Homma and Takeno [Hom] who were able to find soliton solutions with the sinusoidal coupling coming from the dipolar interaction. Another plane-base rotator model has been developed by Yakushevich [Yak2] who considered more specifically the $H-$bonds connecting the bases in a pair to derive another type of interaction potential between paired bases. We discuss this Y model and some extensions in details below. The variety of the rotator models shows clearly that the appropriate degree of description of the rotations of the bases pairs, and the associated potential models, is far from being currently mastered. The question of the selection of potential parameters is especially difficult and controversial since the potentials which appear in these models lump a large variety of microscopic interactions and solvent effects in a single expression\fnote{The question of the choice of appropriate potential parameters gave rise to a sharp debate in the literature. See [Tec1, Tec2, Zan1, Zan2].}. However, it is reassuring to notice that, in spite of their differences, the different models give rise to similar types of localized excitations which can describe open states either in terms of single kinks or as kink-antikink pairs. These states are similar to the local openings introduced by Englander et al. [Eng]. They stay permanently open, and, when they move along the helix, they open-up bases pairs in front while others in their wake close. An extension of the same ideas has been proposed by Zhou and Zhang [Zho] in a model that combines conformational, rotational, longitudinal and transverse motions. Another approach was introduced by Peyrard and Bishop [Pey1] who analyzed the nonlinear dynamics and statistical mechanics of a model which originated from the work of Prohofsky [Pro2] on thermal denaturation of DNA. Rather than the base rotations, this model - which is discussed further in the next sections - uses the stretching of the hydrogen bonds as the main variable, but the essential difference with the approaches discussed above is that the bases in a pair are connected by a Morse potential so that {\it they have no stable open state}, except when the molecule is completely denaturated. Therefore, in this view, the local openings of DNA are not treated as permanently open states that move along the helix, but as fluctuational openings, or ``breathing modes''. As temperature increases, these breathers grow in amplitude while their frequency drops dramatically as the denaturation temperature is approached. They appear as precursor states of the denaturation. A similar view was then chosen by Muto et al. [Mut5] and Techera et al. [Tec1]. These simple models ignore some important aspects of the molecule, and particularly the role of the discreteness which can affect significantly the propagation of lattice solitons, or the inhomogeneities of the molecule. These effects have been considered in some recent works [Dau2, Sal, Tec2, For1]. In this review, we shall focus our attention on two simple models for rotation and stretching transverse motions of DNA (also called torsional and vibrational dynamics) which are representative of the two types of models we have discussed: the Yakushevich model (Y-model) which has permanently open states, and the Peyrard Bishop model (PB-model) which has only breathing excitations. Both these models assume that the molecule is homogeneous, i.e. the differences between the base-pairs in a sequence are ignored; moreover, they consider separately the vibrational degrees of freedom and the torsional dynamics of the molecule. It must be stressed that while the reduction to two degrees of freedom per base discussed above - namely to rotations of the bases in their plane and stretching of the H bonds associated to sugar puckering - can be rather well justified physically, the present splitting of the dynamics in two noninteracting parts is admittedly due to an attempt to deal with a manageable problem rather than to physical considerations. We will see anyway that, somehow rather surprisingly, these "one degree of freedom models" are able to capture relevant experimentally observed features of DNA dynamics; this will give a justification {\it a posteriori} to the splitting. \titleb{\tbn Common features of DNA models} Let us now shortly describe some features common to all the models we are going to discuss in the following sections. We write down the models in term of kinetic ($T$) plus potential ($V$) energy\fnote{We stick now to the Hamiltonian formulation, but at a later time we will also need the Lagrangian formulation of the models}. If $\phi_n^\i$ ($i=1,2$, $n \in \Z$) is the generalized coordinate describing the base at site $n$ on the $i$-th chain of the double helix - which we suppose of infinite length - the Hamiltonian is written as $$ H = T (\.\phi ) + V (\phi ) \eqno(1) $$ where we assume $$ T = \sum_{i,n} \unm \~m (\.\phi_n^\i )^2 \eqno(2) $$ with $\~m$ the mass or the moment of inertia of the bases, according to the degree of freedom taken in consideration. As for the potential energy, it can be divided into three parts: a pairing - or {\it transverse} - interaction energy $V_T$, a stacking interaction energy $V_S$, and a third kind of energy (helicoidal interaction energy) $V_H$ to be discussed later on. In chapter IV, $V_H$ will just be set to zero, while in chapter V we will consider how the introduction of $V_H$ affects the models. So we have $$ V = V_T + V_S \eqno(3) $$ (models which satisfy (3), i.e. without the "helicoidal" interaction, will be called {\it planar}), and we assume that $$ V_T = \sum_n U_T (\phi_n^{(1)} , \phi_n^{(2)} ) ~~;~~ V_S = \sum_{i,n} U_S ( \phi_{n+1}^\i , \phi_n^\i ) \eqno(4) $$ We use coordinates such that $\{ \phi_n^\i = 0 \}$ correspond to the equilibrium stable configuration. It should be stressed that, since the stacking is considerably stronger than H bonding, the nonlinear effects will appear essentially only in the direct pairing (transverse) interactions. This amounts to say that the dynamics stretches the stacking interactions much less that the H bonds. These observations suggest to consider an harmonic potential approximation for stacking interactions. However, a recent work considering anharmonic stacking interactions [Dau3] has shown that nonlinear stacking can have a dramatic effect on the {\it thermodynamics} of the transition; anyway, it was observed from numerical simulations that it has only a small effect on the {\it dynamics} of the nonlinear excitations, except in the immediate vicinity of the denaturation transition. \vfill \eject \titlea{IV - PLANAR MODELS OF DNA DYNAMICS} \def\tbn{\global\advance\titbnumber by 1 IV - \the\titbnumber . } We will discuss here two models of torsional and vibrational DNA dynamics, proposed respectively by Yakushevich and by Peyrard and Bishop. As already recalled in the previous discussion, earlier proposals of similar models were formulated by several authors; we will not discuss these here, and concentrate instead on these two models, denoted respectively as Y and PB model, as representative of the two classes (torsional and vibrational) of models. \titleb{\tbn A model of DNA torsion dynamics} In the model proposed by L. Yakushevich [Yak2], the pairing interaction is considered in the approximation of harmonic potential, depending on the distance between the endpoints of the bases, i.e. the atoms bridged by H bonds. Although the force is assumed to be linear in this distance (which is not too realistic for large amplitude motions), when we describe the position of the bases by means of an angle $\th$ of rotation around the N-C bond the force becomes a nonlinear function of such an angle. The details of pairing interactions in such a model are shown in fig.8. Here the distance $\ell$ between points $p_1$ and $p_2$ satisfies $$ \ell^2 = \[ d + 2r - r( \cos \th_1 + \cos \th_2 ) \]^2 + \[ r (\sin \th_1 + \sin \th_2 ) \]^2 \eqno(1) $$ It is then convenient to pass to coordinates $$ \phi_\pm = {\th_1 \pm \th_2 \over 2 } \eqno(2) $$ and, following Yakushevich, consider the approximation $ d \to 0$. In other words we are considering the stretching of H bonds due to small rotations of the bases and extend this to the full range $[0, 2 \pi ]$ of angular coordinates. Notice that we are not considering the directionality effects, as only distance is taken into account. This feature is nonphysical, as we know that hydrogen bonds are directional (although, as mentioned in sect.II.4, those between pairing bases are not strongly directional); on the other side, in this way we expect to give a description which is fairly accurate, also quantitatively, for small oscillations around the equilibrium position, and which is still qualitatively correct for larger deviations (the stretching of the H bonds requires much more energy than a change in their orientation). Once again, the justification of such a model will lie in the fact it reproduces the phenomena - and physical scales - we aim to model. The potential describing the pairing interaction will then be $$ U_T = \unm \a [ 1 + \cos^2 \phi_- - 2 \cos \phi_+ \cos \phi_- ] \eqno(3) $$ (notice that this is $2 \pi$-periodic in both angular coordinates, as it should. As for stacking interactions, we will also keep the harmonic approximation for them; if $\th_n^\i$ is the angular coordinate relative to the $n$-th base along the $i$-th chain ($i=1,2$), the stacking interaction potential between bases at sites $n$ and $n+1$ will be $$ U_S = \unm \b \[ ( \th_{n+1}^{(1)} - \th_n^{(1)} )^2 + (\th_{n+1}^{(2)} - \th_n^{(2)} )^2 \] \eqno(4) $$ Passing again to the coordinates (2), the above reads $$ V_S =\unm \b \[ ( \phi_{n+1}^{(+)} - \phi_n^{(+)} )^2 + ( \phi_{n+1}^{(-)} - \phi_n^{(-)} )^2 \] \eqno(5) $$ In the following, in order to avoid confusion and to lighten notation, we will prefer to write $\psi$ for $\phi_+$ and $\chi$ for $\phi_-$; with this notation the planar Yakushevich model is defined by the Hamiltonian $$ H = T + V_S + V_T \eqno(6) $$ $$ T = \sum_n \unm I \[ ( \.\th_n^{(1)} )^2 + ( \.\th_n^{(2)} )^2 \] = \sum_n \unm I \( \.\psi_n^2 + \.\chi_n^2 \) \eqno(7) $$ $$ V_S = \sum_n \unm \b \[ ( \psi_{n+1} - \psi_n )^2 + ( \chi_{n+1} - \chi_n )^2 \] \eqno(8) $$ $$ V_T = \sum_n \unm \a \[ 1 + \cos^2 \chi_n - 2 \cos \psi_n \cos \chi_n \] \eqno(9) $$ where we introduced the $\psi , \chi$ coordinates in the kinetic term $T$ as well. The equilibrium (stable) configuration is represented by $\{ (\psi_n , \chi_n ) = (0,0) ~ \forall n \}$; $I$ is the moment of inertia of the bases around the N-C bond connecting them to the sugar backbone. The above hamiltonian gives the equations of motion $$ \eqalign{ I {\ddot \psi_n} =& \b ( \psi_{n+1} - 2 \psi_n + \psi_{n-1} ) - \a \sin \psi_n \cos \chi_n \cr I {\ddot \chi_n} =& \b (\chi_{n+1} - 2 \chi_n + \chi_{n-1} ) - \a \sin \chi_n ( \cos \psi_n - \cos \chi_n ) \cr } \eqno(10) $$ which for small oscillations reduce to $$ \eqalign{ I {\ddot \psi_n} =& \b ( \psi_{n+1} - 2 \psi_n + \psi_{n-1} ) - \a \psi_n \cr I {\ddot \chi_n} =& \b (\chi_{n+1} - 2 \chi_n + \chi_{n-1} ) \cr } \eqno(11) $$ \titleb{\tbn The continuum limit of planar Yakushevich model} Let us now consider the continuum limit of this model: the arrays $\{ \psi_n (t) \}$, $\{ \chi_n (t) \}$ will then be substituted by fields $\psi (x,t) , \chi (x,t)$ which take values in the circle $S^1 = [0,2 \pi ]$ and such that $$ \psi ( n \d , t ) = \psi_n (t) ~~;~~ \chi ( n \d , t ) = \chi_n (t) \eqno(1) $$ where $\d$ is the unit distance between successive bases along the axis of the double helix, $\d \simeq 3.4 \A$. It should be stressed that in considering the continuum approximation, it makes no sense to consider wavelength $\l$ shorter than the interbase distance, i.e. than $\d$. We can now rewrite the various terms of the hamiltonian in the field formulation; we get $$ \eqalignno{ T =& \und \int \unm I ( \psi_t^2 + \chi_t^2 ) \dd x &(2) \cr V_T =& \und \int \unm \a [1 + \cos^2 \chi - 2 \cos \psi \cos \chi ] \dd x &(3) \cr V_S =& \und \int \unm \d^2 \b \[ ( \grad \psi )^2 + ( \grad \chi )^2 \] + O (\d^4 ) \dd x &(4) \cr } $$ where $\grad \equiv \pa / \pa x$ and the $(1/\d )$ factors are only for normalization purposes, i.e. to keep the normalization of the field hamiltonian equal to that of the discrete one. The evolution equations for the fields $\psi , \chi $ issuing from this hamiltonian (or, equivalently, the continuum limit of the equations of motion for $\psi_n , \chi_n$) are $$ \eqalign{ I \psi_{tt} =& \b \d^2 \Lapl \psi - \a \sin \psi \cos \chi + O(\d^4 ) \cr I \chi_{tt} =& \b \d^2 \Lapl \chi - \a \sin \chi ( \cos \psi - \cos \chi ) + O(\d^4 ) \cr } \eqno(5) $$ where $\Lapl$ is the spatial Laplacian, $\Lapl \equiv \pa^2 / \pa x^2 $ We stress that here $\d$ is a physical parameter, fixed by nature, so that we cannot simply consider the limit $\d \to 0$ in order to get rid of the $O(\d^4 )$ terms. We will not discuss this point in detail here, but just state that if $\l_0$ is the shorter wavelength appearing in the Fourier decomposition of $\psi$, we have to ask not only $\l_0 > \d$, but $$ \l_0 >> 2 \d \eqno(6) $$ \titleb{\tbn Solitons in planar Yakushevich model} We want now to focus our attention on soliton solutions for the Yakushevich model. The soliton solutions propagate without changing their form, i.e. we look for $$ \psi (x,t) = \psi (x - vt) \equiv \psi (z) ~~;~~ \chi (x,t) = \chi (x - vt) \equiv \chi (z) \eqno(1) $$ This ansatz, once used into the field evolution equations, produces the {\it equations for Yakushevich solitons}: $$ (v^2 - \k ) \psi_{zz} = - \sin \psi \cos \chi ~~;~~ (v^2 - \k ) \chi_{zz} = - \sin \chi ( \cos \psi - \cos \chi ) \eqno(2) $$ Corresponding, the hamiltonian and lagrangian for Yakushevich solitons read $$ \eqalignno{ H_S =& {I \over \d} \int \unm (v^2 + \k ) \[ (\grad \psi )^2 + (\grad \chi )^2 \] + W(\psi , \chi ) \dd x &(3) \cr L_S =& {I \over \d} \int \unm (v^2 - \k ) \[ (\grad \psi )^2 + (\grad \chi )^2 \] - W(\psi , \chi ) \dd x &(4) \cr } $$ where now $\grad \equiv \pa / \pa z$. On physical grounds, we would require that for $z \to \pm \infty$ the fields go to the equilibrium (zero-energy) configuration; due to the $2 \pi$ periodicity of $W$ in both $\psi$ and $\chi$, this means $$ \psi_z ( \pm \infty ) = \chi_z ( \pm \infty ) = 0 ~~;~~ \psi ( \pm \infty ) = 2 \pi n_\pm ~~;~~ \chi (\pm \infty ) = 2 \pi m_\pm \eqno(5) $$ With no loss of generality we can assume $ n_- = m_- = 0 $ so that the finite energy solutions are classified by $(n_+ , m_+ ) \equiv (n,m)$ alone, i.e. by two integers. This corresponds indeed to the topological classification of solitons $\R \to S^1 \times S^1$, as $n$ and $m$ are {\it rotation numbers} for the fields $\psi$ and $\chi$. Solitons in the $(n,m)$ class will be topologically stable, and we can look for them by minimizing the energy with prescribed boundary conditions. We will not, anyway, use the topological theory of solitons [Dub], and keep the discussion at a simple level. By interpreting $z$ as a time, and $\~\mu \equiv (v^2 - \k )$ as an effective mass, (2) aquires the meaning of a mechanical (Newton) equation for a particle of mass $\~\mu$ moving in the two dimensional plane (with coordinates $\psi , \chi$) under the effects of the potential $W$. The boundary conditions given above fix then the total energy to be $E=0$; it is then immediate to see that nontrivial solutions can exist only if the effective mass is negative; this also gives an upper bound on the speed of solitons: $$ \~\mu = v^2 - \k < 0 ~~;~~ |v| < \sqrt{\k } = v_{{\rm max}} \eqno(6) $$ which in turn implies $$ \mu \equiv - \~\mu \le \k = \mu_{{\rm max}} \eqno(7) $$ From now on we will use $\mu = - \~\mu = \k - v^2 $; notice that in terms of this effective mass, the effective potential will be $-W$. Considering again equations (2), and performing standard manipulations [Gae2], we get now $$ \unm \mu { \dd ~ \over \dd z} \( \psi_z^2 + \chi_z^2 \) = { \dd ~ \over \dd z} W( \psi , \chi ) \eqno(8) $$ so that integrating - and using the boundary conditions (5), which sets the integration constant to zero - we get $$ \unm \mu \[ ( \grad \psi )^2 + ( \grad \chi )^2 \] = W (\psi , \chi ) \eqno(9) $$ This relation is satisfied by soliton solutions $( \=\psi , \=\chi )$; using it in (3), and writing $\nu = \k + v^2$, we have for the energy of such a solution $$ E = \unm { I \over \d} \int (\nu + \mu ) \[ ( \grad \=\psi )^2 + ( \grad \=\chi )^2 \] \dd z ~~ = {I \k \over \d } \int \[ ( \grad \=\psi )^2 + ( \grad \=\chi )^2 \] \dd z \eqno(10) $$ In the simplest cases, i.e. for the $(1,0)$ and the $(0,1)$ soliton solutions, it is actually possible to give an analytic expression for the soliton solutions. This is due to the fact that in this case only one of the fields is nonzero, so that we are dealing with a one-dimensional mechanical problem; a glance at the effective potential $-W$ suffices to get convinced that $(n,0)$ and $(0,n)$ soliton solutions with $|n| > 1$ require {\it both} fields to be not identically zero\fnote{Obviously, this is {\it a fortiori} true if both the soliton numbers $(n,m)$ are not zero.}, so that we can not reduce to a one dimensional motion. We will not give here the details of the (standard) computation of these explicit solutions, but just their final form; the reader is referred to [Gae2] for details. In the case of the (1,0) soliton we have, with $ A = \sqrt{2 \a / \mu }$, the explicit solution $$ \psi (z) = \arccos \[ { \sinh^2 (Ax) - 1 \over \sinh^2 (Ax) + 1 } \] \eqno(11) $$ For the (0,1) soliton we have, with $B = 2 \sqrt{\a / \mu}$, $$ \chi (z) = \arccos \[ {B^2 z^2 -1 \over B^2 z^2 + 1 } \] \eqno(12) $$ It should be remarked that both (11) and (12) obey simple scaling relations under changes of $\a , \mu$ and of the length scale for $z$; this is not a fortuitous case, but indeed corresponds to general scaling properties of soliton solutions for the {\it planar} Yakushevich model, as discussed elsewhere [Gae2]. \vfill \eject \titleb{\tbn A model of DNA vibration dynamics} We now consider the "vibrational" dynamics of DNA, associated to sugar puckering modes; as recalled before, this corresponds to radial motion of the bases with respect to the axis of the double helix; the coordinate $x_n^\i$ describing the position of the base will then be its distance from the axis of the double helix, or better - in order to keep $\{ x_n^\i = 0 \}$ as the equilibrium configuration - this distance minus the distance corresponding to the equilibrium position. We will consider such a model, which was proposed by Peyrard and Bishop [Pey1]; they suggested to use as pairing potential the Morse potential: namely, if $x_n^\i$ is the coordinate defined above for the $n$-th base on the $i$-th chain, we have $$ U_T = D \[ e^{- \a (x^{(1)} + x^{(2)} )/2} - 1 \]^2 = D \[ e^{- \a \psi} - 1 \]^2 \eqno(1) $$ where in the last expression we introduced again the coordinates $$ \phi_\pm = { x^{(1)} \pm x^{(2)} \over 2 } ~~;~~ \phi_+ \equiv \psi ~,~ \phi_- \equiv \chi \eqno(2) $$ As for stacking interactions, we keep the harmonic potential approximation, which gives $$ U_S = \unm \b \[ \( x_{n+1}^{(1)} - x_n^{(1)} \)^2 + \( x_{n+1}^{(2)} - x_n^{(2)} \)^2 \] = \unm \b \[ \( \psi_{n+1} - \psi_n \)^2 + \( \chi_{n+1} - \chi_n \)^2 \] \eqno(3) $$ Finally, the hamiltonian $$ H = T + V_S + V_T \eqno(4) $$ for this (Peyrard-Bishop, or PB) model will be given by $$ \eqalignno{ T =& \sum_{i,n} \unm m \( x_n^\i \)^2 = \sum_n \unm m \( \.\psi_n^2 + \.\chi_n^2 \) &(7) \cr V_S =& \sum_n \unm \b \[ \( \psi_{n+1} - \psi_n \)^2 + \( \chi_{n+1} - \chi_n \)^2 \] &(8) \cr V_T =& \sum_n \unm D \[ e^{-\a \psi_n } -1 \]^2 &(9) \cr } $$ Here $m$ is the mass of the bases; the parameters $\a , \b$ are {\it not} the same as in the Yakushevich model, as we deal with different motions. Notice that the quadratic part of $V_T$ is simply $$ V_T^{(2)} = \sum_n \unm ( D \a^2 ) \psi_n^2 \eqno(5) $$ The above hamiltonian gives the equations of motion $$ \eqalign{ m {\ddot \psi_n} =& \b ( \psi_{n+1} - 2 \psi_n + \psi_{n-1} ) - \a D \[ e^{- \a \psi_n } - 1 \] \cr m {\ddot \chi_n} =& \b (\chi_{n+1} - 2 \chi_n + \chi_{n-1} ) \cr } \eqno(6) $$ which for small oscillations reduce to $$ \eqalign{ m {\ddot \psi_n} =& \b ( \psi_{n+1} - 2 \psi_n + \psi_{n-1} ) - \a^2 D \psi_n \cr m {\ddot \chi_n} =& \b (\chi_{n+1} - 2 \chi_n + \chi_{n-1} ) \cr } \eqno(7) $$ \titleb{\tbn The continuum limit of planar Peyrard-Bishop model} In order to take the continuum limit of this model, we proceed similarly to the case of the Yakushevich model; i.e., the arrays $\psi_n (t) , \chi_n (t)$ are substituted by fields $\psi (x,t) , \chi (x,t)$ such that $$ \psi (n \d ,t ) = \psi_n (t) ~~;~~ \chi ( n \d ,t) = \chi_n (t) \eqno(1) $$ with $\d$ the distance between successive bases along the axis of the double helix, $\d \simeq 3.4 \A$. Notice that $\chi$ takes value on the whole real line $\R$, while for $\psi$ the physical range is $\psi \ge -2 d_0$, where $d_0$ is the distance of the bases from the axis of the double helix at equilibrium. Actually, for both $\psi$ and $\chi$ it makes no sense to consider too big absolute values; the relevant point here for later discussion is that, contrary to the Yakushevich case, the topology of the target space is trivial (this means in particular that we have no topological solitons). With the same remarks as in the Yakushevich case, we can write the hamiltonian in the field formulation as $$ \eqalignno{ T =& \und \int \unm m \( \psi_t^2 + \chi_t^2 \) \dd x &(2) \cr V_T =& \und \int \unm D \( e^{- \a \psi} - 1 \)^2 \dd x &(3) \cr V_S =& \und \int \unm \d^2 \b \[ ( \grad \psi )^2 + ( \grad \chi )^2 \] \dd x &(4) \cr } $$ and the corresponding field evolution equations are $$ \eqalign{ m \psi_{tt} =& \b \d^2 \Lapl \psi + \a D \( e^{-\a \psi} - 1 \) \cr m \chi_{tt} =& \b \d^2 \Lapl \chi \cr } \eqno(5) $$ We recall that these are valid if and only if the solutions have wavelengths large compared to $\d$. \vfill \eject \titlea{V - HELICOIDAL MODELS OF DNA DYNAMICS} \def\tbn{\global\advance\titbnumber by 1 V - \the\titbnumber . } \titleb{\tbn Helical geometry and helicoidal DNA models} The models discussed up to now have been qualified as {\it planar}; this is due to the fact that they can be schematized as in fig.9, which also means that they don't take into account the peculiar geometry of the double helix. It turns out that taking this into account, even in a "minimal" way, produces qualitative changes in the dynamics, quite relevant also in relation with "observable" properties of the model (especially in the case of the PB model). The modification we consider here was first proposed in the context of the Yakushevich model, and then applied also to the Peyrard- Bishop one [Gae1,Dau1]; the models modified in this way will be qualified as {\it helicoidal}. In the following, we will write hY for "helicoidal Yakushevich", and hPB for "helicoidal Peyrard- Bishop". In a double helix (with pitch of the helix $2h$ in base units) it happens that bases which are half-wind of the helix apart, i.e. $h$ sites away, come to be spatially near, so that they can interact; this fact is shown schematically in fig.10. It should be stressed that such an interaction is actually present in real DNA, where it is realized by bridging water filaments between the concerned bases. From the point of view of our models, and referring to the notation previously employed, this means that the hamiltonian will now be of the form $$ H = T + V_T + V_S + V_H \eqno(1) $$ where $V_H$ represents these new or - as they will be called - {\it helicoidal} interactions. These will appear through terms of the form, referring again to section IV.1 for notation, $$ U_H = k \( \phi_n^{(1)} , \phi_{n+h}^{(2)} \) + k \( \phi_n^{(2)} , \phi_{n+h}^{(1)} \) \eqno(2) $$ It should be noticed that helical geometry poses some constraint on the choice of coordinates for helicoidal models [Gae3]. The choice of coordinates made above for both models we consider took this problem into account, and is therefore correct. As mentioned above, the helicoidal interaction is due in real DNA to bridged water filaments and therefore again to H bonds; these bonds will be stretched much less than the H bonds pairing bases at corresponding sites on different helices (i.e. than the bonds modelled by $V_T$), so that it is reasonable to assume they remain in the harmonic regime, and the helicoidal interaction can be modelled by an harmonic potential (as we do below). As it is the case for all the interactions in our simple models, also the $V_H$ introduced here should be seen as effective ones; therefore one should be most careful in extending the above discussion, which justified the harmonic approximation, to quantitative matters, i.e. to an estimation for the ratio of the coupling constants relative to $V_T$ and $V_H$ which should be $\sim 10$ by naive arguments based on counting the bonds involved in the water filament which mediates the helicoidal interaction in real DNA [Gae5]. Indeed, in the context of the hPB model, Dauxois gave quite different estimates (a ratio of $\sim 3$) based on less naive arguments; with this estimate, he could obtain agreement with experimental findings [Dau1], and we will use his values in the following computations. We will leave the coupling constant for $V_H$, as already done for the other constants appearing in the model, as an effective constant to be adjusted. We also stress that, given the relatively long range of this effective interaction, it is potentially even more sensitive than the other ones to the effect of the complicate interactions in the molecule and with its environment. \titleb{\tbn Introduction of helicoidal terms in the Peyrard- Bishop model} We will begin by shortly considering the effect of introducing "helicoidal" terms in the PB model, as was first done by Dauxois [Dau1]. We add to the PB hamiltonian a term $V_H$ which, coherently with the discussion of the previous section, is simply an harmonic potential: $$ U_H = \unm \g \[ \( x_{n+h}^{(1)} - x_n^{(2)} \)^2 +\( x_{n+h}^{(2)} - x_n^{(1)} \)^2 \] \eqno(1) $$ Upon introduction of the $\psi, \chi$ coordinates, so that $$ x_n^{(1)} = \psi_n + \chi_n ~~;~~ x_n^{(2)} = \psi_n - \chi_n \eqno(2) $$ the interaction energy reads, with straightforward algebra, $$ U_H = \unm \g \[ \( \psi_{n+h} - \psi_n \)^2 + \( \chi_{n+h} + \chi_n \)^2 \] \eqno(3) $$ Notice the different signs in the $\psi$ and $\chi$ parts of this expression. The total hamiltonian of the helicoidal Peyrard- Bishop (hPB) model reads therefore $$ H = T + V_S + V_T + V_H \eqno(4) $$ with $T, V_S , V_T $ as in the PB model and $$ V_H = \sum_n \unm \g \[ \( \psi_{n+h} - \psi_n \)^2 + \( \chi_{n+h} + \chi_n \)^2 \] \eqno(5) $$ with $h$ the half-pitch of the helix in base units (i.e. $h=5$ for ordinary DNA). The hPB hamiltonian gives as equations of motion $$ \eqalign{ m {\ddot \psi_n} =& \b ( \psi_{n+1} - 2 \psi_n + \psi_{n-1} ) + \g ( \psi_{n+h} - 2 \psi_n + \psi_{n-h} ) - \a D \[ e^{- \a \psi_n } - 1 \] \cr m {\ddot \chi_n} =& \b ( \psi_{n+1} - 2 \psi_n + \psi_{n-1} ) - \g ( \psi_{n+h} - 2 \psi_n + \psi_{n-h} ) \cr } \eqno(6) $$ which for small oscillations reduce to $$ \eqalign{ m {\ddot \psi_n} =& \b ( \psi_{n+1} - 2 \psi_n + \psi_{n-1} ) + \g ( \psi_{n+h} - 2 \psi_n + \psi_{n-h} ) - \a^2 D \psi_n \cr m {\ddot \chi_n} =& \b (\chi_{n+1} - 2 \chi_n + \chi_{n-1} ) - \g ( \psi_{n+h} - 2 \psi_n + \psi_{n-h} ) \cr } \eqno(7) $$ Let us now consider the continuum limit of this model. Clearly the helicoidal interaction leads to the appearance of nonlocal (in $x$, the coordinate along the axis of the double helix) terms; the field evolution equations will be, understanding that $\psi \equiv \psi (x,t)$, $\chi \equiv \chi (x,t) $ unless otherwise indicated, and writing $w = h \d$, $$ \eqalign{ m \psi_{tt} =& \b \d^2 \Lapl \psi + \g \W_+ (\psi ) + \a D \( e^{-\a \psi} - 1 \) \cr m \chi_{tt} =& \b \d^2 \Lapl \chi - \g \W_- (\chi ) \cr } \eqno(8) $$ with the linear nonlocal operator $\W_\pm$ defined as $$\W_\pm (\phi ) = \phi (x+w , t) \mp 2 \phi (x,t) + \phi (x-w,t) \eqno(9) $$ \titleb{\tbn Dispersion relations for the helicoidal Peyrard- Bishop model} We set $$\k = \b \d^2 / m ~~,~~ \s = D \a^2 / m ~~,~~ \r = \g / m ~;$$ with these the equations for small oscillations in the hPB model are $$\eqalign{ \psi_{tt} =& \k \Lapl \psi + \r \W_+ (\psi ) - \s \psi \cr \chi_{tt} =& \k \Lapl \chi - \r \W_- (\chi ) \cr} \eqno(1) $$ Moving to Fourier components, i.e. to $$ \eqalign{ \psi (x,t) =& F_{q,\om} e^{i qx} e^{i \om t} \equiv F_{q,\om} | q, \om > \cr \chi (x,t) =& G_{q,\om } e^{i qx} e^{i \om t} \equiv G_{q,\om } | q, \om > \cr} \eqno(2) $$ we have to take care of the action of the nonlocal operators $\W_\pm$ on $| q, \om >$: it is given by $$ \W_\pm \( | q, \om > \) = \( e^{iqw} \mp 2 + e^{-iqw} \) | q, \om > = 2 \[ \cos (qw) \mp 1 \] | q, \om > \eqno(3) $$ Consequently, the dispersion relations are immediately seen to be $$ \eqalign{ \om_\psi^2 =& \k q^2 + 2 \r \[1 - \cos (qw) \] + \s \equiv \k q^2 + 4 \r \sin^2 \({qw \over 2 } \) + \s \cr \om_\chi^2 =& \k q^2 + 2 \r \[1 + \cos (qw) \] \equiv \k q^2 + 4 \r \cos^2 \({qw \over 2 } \) \cr } \eqno(4)$$ The dispersion relations (4) are plotted in fig.11 (with values of the parameters following the estimates given by [Dau1]). It is immediately clear, by comparing with the dispersion relations of the planar PB model, that the introduction of helicoidal terms has drastically changed the (small oscillations) dynamics. Notice that in the PB and hPB models the $\psi$ and $\chi$ fields do exactly decouple, so that we are actually in presence of separate dynamics; we stress that this is true not only for small oscillations, but at the full nonlinear level as well. Let us remark in particular some features of these dispersion relations: {\it a)} Now there is no mode with $\om (q) =0 $ for $q \to 0$; this means that we have a threshold for creation of excitations of given (small) amplitude. Notice that the acoustic branch in the planar model corresponded to the $\chi$ field; long wavelength $\chi$ excitations are energetically costly due to the plus sign in the $\chi$ part of $V_H$, see eq.(5) in section V.2. {\it b)} The crossing of the two branches\fnote{As already mentioned, in the PB model the two branches are non interacting, so that the crossing points have no physical meaning; the situation will be different in the context of the Y model.}, and their relative ordering for $q=0$, are controlled by the sign of the quantity $$ s = 1 - {\a \over 2 \g } $$ {\it c) } the $\psi$ branch has absolute minimum at $q=0$, while the $\chi$ branch has an absolute minimum at $q_0 \not= 0$ ($q_0 \simeq 4 \pi / w$, which corresponds to a wavelength $\l_0 \simeq w/2$). This latter fact is an obvious consequence of the introduction of the helicoidal interaction and of the signs appearing in the $\chi$ part of $V_H$: $\chi$ waves with $\chi_{n+h} = - \chi_n$ cost less energy. \titleb{\tbn Nondispersive waves and breather modes in the helicoidal Peyrard- Bishop model} The waves corresponding to extremal points of the dispersion relations have zero group velocity; i.e., they are {\it standing waves} or, in other words, correspond to {\it breather modes}: there are oscillatory local openings of the double helix which remain localized without spreading or travelling. The PB model is aimed at describing the denaturation of DNA; this process arises indeed by local opening of "bubbles" in the double helix, i.e. opening of consecutive pairs of bases (see the discussion in previous sections). The experiments tell us that the size of the bubbles is of a few base units, and the time frequencies involved are of the order of picoseconds; it is hoped that the breather modes of the hPB model can indeed model this process. We discuss here numerical findings\fnote{More complete and refined numerical simulations of this model are being performed; these also concern estimates of parameters [Dau4].} based on parameter estimates given in [Dau1]. It turns out that the hPB model is rather well compatible with experimental findings. Notice also that, due to the virial theorem, the energy of waves of small amplitude $A$ and obeying the dispersion relation $\om = \om (q)$ is just $$ E = A^2 m \om^2 (q) \eqno(1) $$ Let us now write for ease of notation $$ \xi = { qw \over 2} \eqno(2) $$ and rewrite the dispersion relations, in terms of the physical coupling parameters, as $$ \eqalign{ \Om_\psi (\xi ) \equiv & m \om^2_\psi ( q(\xi )) = 4 \g \[ B \xi^2 + C + \sin^2 (\xi ) \] \cr \Om_\chi (\xi ) \equiv & m \om^2_\chi ( q(\xi )) = 4 \g \[ B \xi^2 + \cos^2 (\xi ) \] \cr } \eqno(3) $$ where the parameters $B,C$ are given by $$ B = {\b \over h^2 \g } ~~;~~ C = {\a D \over 4 \g } \eqno(4) $$ The extremal points $\xi = \xi_0 $ of these are obtained by solving $\pa \Om / \pa \xi = 0$, namely $$ 2 B \xi \mp \sin (2 \xi ) = 0 \eqno(5) $$ with the minus sign for the $\psi$ branch and the plus one for the $\chi$ branch. Notice that necessarily $\xi_0 \le 1/B = h^2 \g /\b$. The wavelength corresponding to $\xi_0$ is simply given by $$ \l (\xi_0 ) = { 2 \pi \over q (\xi_0 ) } = { \pi h \over \xi_0 } \d \eqno(6) $$ (Notice that the above bound on $\xi_0$ ensures then $\l_0 \ge [( \pi \b ) / (h \g ) ] \d $; with the value of $\b / \g$ considered in this section, this also ensures $\l_0 >> \d$, as required for the validity of our continuum approximation). The time frequency $\om$ and period $T$ corresponding to $\xi_0$ are $$ \om ( \xi_0 ) = \sqrt{{\Om ( \xi_0 ) \over m }} ~~~~;~~~~ T ( \xi_0 ) = {2 \pi \over \om ( \xi_0 ) } \eqno(7) $$ >From (1), and expecting each field mode to have an average energy $ = \unm k T$, we have ($T \simeq 37 ~^0 C$) the expected amplitude of oscillations corresponding to $\xi_0$: $$ A ( \xi_0 ) = \sqrt{ { k T /2 \over \Om ( \xi_0 ) }} \eqno(8) $$ As for the parameters $\a , \b , \g , m$, we stick to the estimates given by Dauxois [Dau1] and recently refined in [Dau4]; adapting the notation to our present one, these are given in table I (we report ${\widetilde \a} \equiv \a^2 D$ rather than $\a$ since the former is the coefficent of the linear term in the equations of motion). The corresponding numerical results are reported, for the $\psi$ and $\chi$ branches, respectively in tables II and III. As anticipated, the orders of magnitude that we obtain in this way for $\l$ and $T$ are in good agreement with experimental observations (notice that opening of the double helix correspond to the $\psi$ mode). \vfill \eject ~ \bigskip $$ \matrix{ {\widetilde \a} & \simeq & 0.79 \ea & \simeq & 1.3 \cdot 10^4 \ec \cr \b & \simeq & 0.06 \ea & \simeq & 9.8 \cdot 10^2 \ec \cr \g & \simeq & 0.02 \ea & \simeq & 3.3 \cdot 10^2 \ec \cr m & \simeq & 300 ~ {\rm m.u.} & \simeq & 5 \cdot 10^{-22} ~ {\rm g} \cr} $$ \medskip \centerline{\petit Table I - Parameter choice for hPB model} \bigskip $$\matrix{ \xi_0 & 0. & 1.4 & 3.7 & 4.0 \cr \Omega (\xi_0 ) & 3.2 ~10^{4} \erg & 8.5 ~ 10^{3}\erg & 1.5 ~10^{4} \erg & 2.2 ~10^{4} \erg \cr \om (q_0 ) & 8.0 ~10^{12} \smu & 4.1 ~10^{12} \smu & 5.5 ~10^{12} \smu & 6.6 ~10^{12} \smu \cr T (q_0 ) & 0.8 \ps & 1.5 \ps & 1.1 \ps & 0.9 \ps \cr A (q_0 ) & 0.08 \A & 0.16 \A & 0.12 \A & 0.10 \A \cr \lambda (q_0 ) & \infty & 11.2 ~\d & 4.2 ~\d & 3.9 ~\d \cr} $$ \medskip \centerline{\petit Table II - Physical quantities ($\chi $ branch, hPB model)} \bigskip $$ \matrix{ \xi_0 & 0. & 1.8 & 2.8 \cr \Omega (q_0 ) & 2.6 ~10^4 \erg & 6.8 ~10^4 \erg & 5.9 ~10^4 \erg \cr \om (q_0 ) & 7.2 ~10^{12} \smu & 1.2 ~10^{13} \smu & 1.1 ~10^{13} \smu \cr T (q_0 ) & 0.9 \ps & 0.5 \ps & 0.6 \ps \cr A (q_0 ) & 9.0 ~10^{-2} \A & 5.6 ~10^{-2} \A & 6.0 ~10^{-2} \A \cr \lambda (q_0 ) & \infty & 8.7 ~\d & 5.6 ~\d \cr}$$ \medskip \centerline{\petit Table III - Physical quantities ($\psi$ branch, hPB model)} \bigskip \vfill \eject \titleb{\tbn Introduction of helicoidal terms in the Yakushevich model} Let us now consider the effects of introducing helicoidal interactions in the Yakushevich model; here again we add to the Yakushevich hamiltonian an harmonic potential modelling the helicoidal interaction; i.e. we have $$ U_H = \unm \g \[ \( \th_{n+h}^{(1)} - \th_n^{(2)} \)^2 +\( \th_{n+h}^{(2)} - \th_n^{(1)} \)^2 \] \eqno(1) $$ Here we change again to normal coordinates $\psi , \chi$, with $$ \th_n^{(1)} = \psi_n + \chi_n ~~;~~ \th_n^{(2)} = \psi_n - \chi_n \eqno(2) $$ in terms of which (1) reads $$ U_H = \unm \g \[ \( \psi_{n+h} - \psi_n \)^2 + \( \chi_{n+h} + \chi_n \)^2 \] \eqno(3) $$ Notice here as well the different signs with which $\psi$ and $\chi$ enter in $U_H$. The total hamiltonian of the helicoidal Yakushevich (or hY) model reads therefore $$ H = T + V_S + V_T + V_H \eqno(4) $$ where $T, V_S , V_T $ are as in the Y model and $$ V_H = \sum_n \unm \g \[ \( \psi_{n+h} - \psi_n \)^2 + \( \chi_{n+h} + \chi_n \)^2 \] \eqno(5) $$ with $h$ the half-pitch of the helix in base units ($h=5$). This hamiltonian gives the equations of motion $$ \eqalign{ I {\ddot \psi_n} =& \b ( \psi_{n+1} - 2 \psi_n + \psi_{n-1} ) + \g ( \psi_{n+h} - 2 \psi_n + \psi_{n-h} ) - \a \sin \psi_n \cos \chi_n \cr I {\ddot \chi_n} =& \b ( \psi_{n+1} - 2 \psi_n + \psi_{n-1} ) - \g ( \psi_{n+h} - 2 \psi_n + \psi_{n-h} ) - \a \sin \chi_n ( \cos \psi_n - \cos \chi_n ) \cr } \eqno(6) $$ and small oscillations are described by the equations $$ \eqalign{ I {\ddot \psi_n} =& \b ( \psi_{n+1} - 2 \psi_n + \psi_{n-1} ) + \g ( \psi_{n+h} - 2 \psi_n + \psi_{n-h} ) - \a \psi_n \cr m {\ddot \chi_n} =& \b (\chi_{n+1} - 2 \chi_n + \chi_{n-1} ) - \g ( \psi_{n+h} - 2 \psi_n + \psi_{n-h} ) \cr } \eqno(7) $$ In order to consider the continuum limit, we introduce again the nonlocal linear operators $$\W_\pm (\phi ) = \phi (x+w , t) \mp 2 \phi (x,t) + \phi (x-w,t) \eqno(8) $$ and write the full equations as $$\eqalign{ \psi_{tt} =& \k \Lapl \psi + \r \W_+ (\psi ) - \s \sin \psi \cos \chi \cr \chi_{tt} =& \k \Lapl \chi - \r \W_- (\chi ) - \s \sin \chi ( \cos \psi - \cos \chi ) \cr} \eqno(9) $$ where we have used $$\k \equiv \b \d^2 / I ~~,~~\r \equiv \g / I ~~;~~ \s = \a / I ~~;~~ w = h \d \eqno(10) $$ with $\d$ the distance between consecutive bases along the axis of the double helix ($\d \simeq 3.4 \A$). As for the small oscillation equations, they are simply given by $$\eqalign{ \psi_{tt} =& \k \Lapl \psi + \r \W_+ (\psi ) - \s \psi \cr \chi_{tt} =& \k \Lapl \chi - \r \W_- (\chi ) \cr} \eqno(11) $$ Notice that these equations are formally analogous to the ones obtained for the hPB model, so that we will not discuss the dispersion relations for the hY model, but just refer to the corresponding discussion of the hPB model. We stress that in the case of the Y and hY models the fields $\psi , \chi$ are decoupled at the level of small oscillations (so that they are again "normal coordinates"), but are instead coupled through nonlinear terms, contrary to what happens in PB and hPB models. Since the $\psi$ and $\chi$ fields are now interacting through nonlinear terms, the "crossing points" (i.e. the points at which the two branches meet) do have physical meaning, corresponding to resonances between the two fields [Gae4]. \titleb{\tbn Solitons in the helicoidal Yakushevich model} As already recalled, the Yakushevich and helicoidal Yakushevich models have been primarily developed to describe travelling solitons, which are supposed to play a central role in the transcription of DNA [Eng]. Therefore, we will concentrate on solitons, and will not consider standing nondispersive waves in the hY model. These waves corresponding to the extremal points of the dispersion relations, would correspond again to breather modes giving local pulsating and nonpropagating openings of the double helix, i.e. breaking of the H bonds pairing bases at same site on different helices. Numerical computations [Gae3] show that indeed experimental data on spontaneous denaturation of DNA are best fitted by the hPB model. The full field equations for the hY model are given by eq.(9) of the preceding section. The corresponding hamiltonian $H$ and lagrangian $L$ are given by $$ \eqalign{ L =& {I \over \d } \int \unm \( \psi_t^2 + \chi_t^2 \) - \unm \k \[ ( \grad \psi )^2 + ( \grad \chi )^2 \] - W (\psi , \chi ) \dd x \cr H =& {I \over \d } \int \unm \( \psi_t^2 + \chi_t^2 \) + \unm \k \[ ( \grad \psi )^2 + ( \grad \chi )^2 \] + W (\psi , \chi ) \dd x \cr } \eqno(1) $$ where now $W( \psi , \chi )$ is made up of a local and a nonlocal ($\N$) part: $$ \eqalign{ & W (\psi , \chi ) = \unm \s \[ 1 + \cos^2 \chi - 2 \cos \psi \cos \chi \] + \N (\psi , \chi ) \cr & \N ( \psi , \chi ) = \unm \r \[ \[ \psi (x+w) - \psi (x) \]^2 + \[ \chi (x+w) + \chi (x) \]^2 \] \cr } \eqno(2) $$ The appearance of nonlocal terms in (1),(2) is the main new qualitative feature encountered in the helicoidal version of the model. Physically, this correspond to the introduction of a second length scale in the model: beside the lattice spacing $\d$, we now have also introduced the (half) pitch of the helix, by the scale $w=h \d$. We want now to concentrate on solitons, i.e. on solutions of the form $$ \psi (x,t) = \psi (x - vt) \equiv \psi (z) ~~;~~ \chi (x,t) = \chi (x - vt) \equiv \chi (z) \eqno(3) $$ Notice that the discussion of boundary conditions conducted in the planar case does still apply here, as well as the related classification of soliton solutions by topological (or soliton) integer numbers $(n,m)$. Inserting the ansatz (3) into (1), we get the hamiltonian and the lagrangian for helicoidal Yakushevich solitons, where $\grad \equiv \pa / \pa z$: $$ \eqalign{ H =& {I \over \d} \int \unm (v^2 + \k ) \[ (\grad \psi )^2 + (\grad \chi )^2 \] + W(\psi , \chi ) \dd z \cr L =& {I \over \d} \int \unm (v^2 - \k ) \[ (\grad \psi )^2 + (\grad \chi )^2 \] - W(\psi , \chi ) \dd z \cr } \eqno(5) $$ We write again $\~\mu = v^2 - \k = - \mu$, and the bounds on $v,\mu $ given in the planar case are still valid here. By the same procedure as in the planar case, we get the equations for solitons as [Gae2] $$ \unm \mu \[ ( \grad \psi )^2 + ( \grad \chi )^2 \] = W (\psi , \chi ) \eqno(6) $$ which again shows, once used in (5), that the energy of soliton solutions $\=\psi , \=\chi$ is given by $$ E = {I \k \over \d } \int \[ ( \grad \=\psi )^2 + ( \grad \=\chi )^2 \] \dd z \eqno(7) $$ \vfill \eject \titleb{\tbn Exact solution for the helicoidal (1,0) soliton} The nonlocal terms make it quite difficult to obtain explicit soliton solutions for the hY model; one can nevertheless get some scaling relations obeyed by these solutions, which are also useful to obtain the explicit form of the simplest of them. The idea is to scale down the length $w = h \d$ to $\d$, which produces an "effective planar lagrangian", and then apply the continuous approximation. The limitation to this approach comes from the special sign in the nonlocal interaction of $\chi$ field, which prevents from obtaining an effective lagrangian of the Yakushevich type, and would introduce terms of the type $\chi \grad \chi$ (see the discussion in [Gae5]). The limitations to work with an "effective planar lagrangian" for the rescaled hY model do not apply if we just consider the (1,0) solitons, due to $\chi \equiv 0$. In this case the explicit soliton solution is just derived from the planar Y model solution by scaling arguments. Indeed, by scaling arguments [Gae2], one proves that by redefining an effective mass $$ \mu_H \equiv \mu - \r \d^2 \equiv (\k - \r ) \d^2 - v^2 \eqno(1) $$ the equations are formally analogous to the equations encountered in the planar case; i.e. we get $$ \unm \mu_H ( \grad \psi )^2 = \a (1 - \cos \psi ) \eqno(2) $$ This can be integrated as shown earlier; notice that in these terms the introduction of helicoidal interactions amounts to a rescaling of the constant $A$ appearing in the explicit solution, which is now $$ A = \sqrt{ { 2 \a \over \mu_H }} = \sqrt{{2 \a \over \mu - \r \d^2 }} \eqno(3) $$ The dependence of (1,0) soliton energy on speed is given by $$ E ( \l , \g ) = { A' \over B' - \g^2 } E ( \d , 0 ) \eqno(4) $$ where $$ A' = 1 + { \g \over \b } ~~;~~B' = 1 + {\g \over h^2 \b } \eqno(5) $$ (see [Gae2] for details). This means that the scaling of energy with speed is qualitatively analogous to the one obtained for the planar Y model; similar considerations apply to the relation of size versus speed of soliton solutions. It should be stressed that the (1,0) solitons correspond to those whose role in DNA transcription was first suggested by Englander [Eng], so that the discussion of the present section is of particular relevance, although concerning only one specific soliton solution. \titleb{\tbn Numerical results for solitons in the helicoidal Yakushevich model} For solitons other than the (1,0) ones, we are not able to determine explicit analytical solutions, nor exact scaling laws at fixed $\l$; we must therefore resort to numerical computations. We fix the parameters $\a , \b , \g , I$ appearing in the model, and let $\mu$ vary between $0$ and $\mu_{{\rm max}}$ (i.e. we consider different values of $\g$; it should be recalled that physically $\g \simeq 10^{-8}$). For any given value of $\mu$, one can determine numerically, on a sufficiently large discrete lattice, the configuration $\=\psi , \=\chi$ which extremalizes - among those with prescribed boundary conditions - the Yakushevich soliton lagrangian (see eq.(2) of section 22). Let $S(\mu )$ be this extremal value of the lagrangian, $$ S(\mu ) = - { I \over \d} \int \unm \mu \[ ( \grad \=\psi )^2 + ( \grad \=\chi )^2 \] + W ( \=\psi , \=\chi ) \dd z \eqno(1) $$ The corresponding energy, by comparing the soliton lagrangian and hamiltonian and recalling that for soliton solutions $$ \unm \mu \[ ( \grad \=\psi )^2 + ( \grad \=\chi )^2 \] = W ( \=\psi , \=\chi ) \eqno(2) $$ is given simply by $$ H ( \mu ) = {\mu_{{\rm max}} \over \mu } S (\mu ) \eqno(3) $$ In terms of $\g$, this reads $$ H (\g ) = {1 \over 1 - \g^2 } S \( \mu ( \g ) \) \equiv {1 \over 1 - \g^2 } S \( (1 - \g^2 ) \mu_{{\rm max}} \) \eqno(4) $$ The minimization of the lagrangian is performed by a Montecarlo algorithm; the values of the parameters are given in table IV and discussed in [Gae2]; the results of the numerical computations are given in fig.12 for the (1,1) soliton (as recalled above, the $\psi$ field is supposed to matter for the solitons we are looking for: we want therefore a nonzero $\psi$ soliton number); here we plot $H(\g )$ in terms of $\g$, together with the best fit by functions of the form $$ H(\g ) = { A \over B - \g^2 } \eqno(5) $$ This shows good qualitative agreement with the situation seen above for the (1,0) soliton, which is in turn qualitatively equivalent to the one for planar Y solitons. \vfill \eject \titlea{VI - THERMODYNAMICS OF THE PB MODEL} \def\tbn{\global\advance\titbnumber by 1 VI - \the\titbnumber . } \def\ket#1{| #1 \rangle} \def\bra#1{\langle #1 |} \medskip \titleb{\tbn Statistical mechanics of the PB model} \smallskip The PB model was designed not only to study nonlinear dynamics of DNA, but also its thermodynamics, and particularly the thermal denaturation of the molecule. Due to the one-dimensional character of the system, and because the interactions are restricted to nearest neighbor interactions, it can be treated exactly, including fully the nonlinearities, with the transfer operator method [Pey1, Dau4]. As discussed above, the hamiltonian of the model is the sum of two decoupled contributions: an hamiltonian of an {\it harmonic} lattice for variable $\chi$ and an hamiltonian containing nonlinear terms for the variable $\psi$. The first part describes the acoustic motions of the molecule (bending modes) and it is well behaved. Since, for denaturation, we are interested in the stretching motion, the statistical mechanics can consider only the part of the hamiltonian which depends on $\psi$, i.e. $$ {H} = \sum_n \Bigr[ {1\over 2} m{\dot \psi_n^2} +{\beta\over 2} (\psi_n-\psi_{n-1})^{2} +D(e^{-\alpha\psi_n}-1)^2\Bigr] \; . \eqno(1) $$ For a chain containing $N$ units with nearest neighbor coupling and periodic boundary conditions, the classical partition function, given in terms of the hamiltonian (1), can be expressed as : $$ {\cal Z} = \int_{-\infty}^{+\infty} { \prod_{n=1}^N d\psi_n dp_n e^{- \beta_B {H}} } = {\cal Z}_p \times \int_{- \infty}^{+ \infty} \prod_{ n= 1}^N d\psi_n e^{- \beta_B f(\psi_n , \psi_{n - 1} ) } = {\cal Z}_p \; {\cal Z}_\psi \; , \eqno(2) $$ where $f(\psi_n , \psi_{n - 1} )$ is the potential part of the Hamiltonian and $\beta_B = 1/(k_B T)$ ($k_B$ beeing the Boltzmann constant). The momentum parts are readily integrated to give the usual kinetic factor for $N$ particles $ {\cal Z}_p= \left( 2 \pi m k_B T \right)^{N/2}$. The potential part can be evaluated exactly in the thermodynamic limit of a large system $(N \to \infty )$ using the eigenfunctions and eigenvalues of the transfer integral operator [Sca, Kru1, Cur] $$ \int d\psi_{n-1} e^{- \beta_B f(\psi_n, \psi_{n-1}) } \phi_i (\psi_{n-1}) = e^{-\beta_B \epsilon_i} \phi_i (\psi_{n}) . \eqno(3) $$ The calculation is similar to the one performed by Krumhansl and Schrieffer [Kru1] for the statistical mechanics of the $\phi^4$ field. It yields $ {\cal Z}_\psi = \exp ({- N \beta_B \epsilon_0 }) $ , where $\epsilon_0$ is the lowest eigenvalue of the operator. The free energy of the model can then be computed as ${\cal F} = - k_B T \ln{\cal Z}=-(Nk_BT/ 2)\ln (2\pi m k_BT)+N \epsilon_0 \; $ and the specific heat $C_v=-T(\partial^2 {\cal F} / \partial T^2)$. The quantity which gives a measure of the extent of the denaturation of the molecule is the mean stretching $\langle \psi_m\rangle $ of the hydrogen bonds, which can also be calculated with the transfer integral method [Pey1], and yields $$ \langle \psi\rangle =\langle \psi_m\rangle = { {\sum}_{i=1}^N \bra{\phi_i(\psi)} \psi \ket{\phi _i(y)} e^{-N \beta \epsilon_i} \over {\sum}_{i=1}^N \bra{\phi_i(\psi)} \phi_i(\psi) \rangle e^{-N \beta_B \epsilon_i} } =\bra{\phi_0(\psi)} \psi \ket{\phi_0(\psi)} = \int \phi_0^2(\psi) \psi d\psi \quad , \eqno(4) $$ since in the limit of large $N$ the result is again dominated by the lowest eigenvalue $\epsilon_0$ associated with the normalized eigenfunction $\phi_0(\psi)$. In the continuum limit, the TI eigenvalue problem can be solved analytically [Pey1], but experiments on proton exchange in DNA [Gue] show some evidence of exchange limited to a single base pair which suggests that discreteness effects can be extremely large in DNA. Therefore the eigenvalue equation of the transfer operator has been solved numerically without approximations [Dau4]. The TI operator is symmetrized and the integral is replaced by sums of discrete increments, using summation formulas at different orders. The problem is then equivalent to finding the eigenvalues and eigenvectors of a symmetric matrix. Figure 13 compares the thermal evolution of $\langle \psi_m \rangle$ obtained with the continuum approximation and the exact numerical calculation, for the model parameters discussed below. Both methods show a divergence of the hydrogen bond stretching over a given temperature, but the melting temperature given by the numerical treatment is significantly higher, pointing out the large role of discreteness in DNA dynamics if one uses realistic parameters for the model. The TI calculation shows also that the specific heat has a broad maximum around the denaturation temperature. \titleb{\tbn Numerical simulations at constrained temperature} Therefore, the thermodynamics of the PB model shows that it exhibits a thermal evolution that is qualitatively similar to the denaturation of the molecule observed experimentally. But this statistical approach does not give informations on the {\it mechanism} of the denaturation, and in particular, if it does start locally by the formation of denaturation bubbles in agreement with the experiments. In order to study this aspect, the dynamics of the model in contact with a thermal bath has been investigated by molecular dynamics simulation with the Nos\'e scheme [Nos]. Beginning with the hamiltonian ${H}$ (eq.(1) of previous section) and the 2N-dimensional phase space of a chain of $N$ base-pairs with periodic boundary conditions, the fixed temperature canonical ensemble can be simulated by the addition of a single variable $s$, which regulates the energy flows, and an additional parameter $M$, which fixes the scale of the temperature fluctuations. Nos\'e demonstrated that in this phase space of an extended hamiltonian $H'$, the microcanonical ensemble of $H'$ is precisely the canonical ensemble of ${H}$ at temperature $T$. This property is only exact for equilibrium properties, but investigations currently in progress [The] show that, provided that the characteristic time of the Nos\'e thermostat, controlled by $M$, is properly chosen, it can also give reliable results for the dynamical properties. Most of the simulations [Dau4] have been performed with a chain of 256 base pairs with periodic boundary conditions, but in order to achieve better statistics, some simulations have been performed on a Connection Machine-200 with 16384 base pairs. Equations are integrated with a 4-th order Runge Kutta scheme with a time-step chosen to conserve $H'$ to an accuracy better than 0.001\% during a run. They have been performed with a system of units adapted to the energy and time scales of the problem. Energies are expressed in eV, masses in atomic mass unit ($a.m.u.$) and length in \AA . The resulting time unit is $1\ t.u. = 1.0214\; 10^{-14}$~s. The choice of appropriate model parameters is a very controversial topic [Zan]. There are well established force fields for molecular dynamics of biological molecules, but they have been designed to provide a good description of the {\it small} amplitude motions of the molecule and are not reliable for the very {\it large} amplitude motions involved in the denaturation. In the model, the Morse potential is an effective potential which links the two strands. It results from a combination of an attractive part due to the hydrogen bonds between two bases in a pair and the repulsive interaction between the charged phosphate groups on the two strands. The potential for the hydrogen bonds can be rather well estimated [Gao] but the repulsive part is harder to determine because the repulsion is partly screened by ions of the solvent. The parameters that are chosen rely on estimations to give realistic properties for the model in terms of vibrational frequencies, size of the open regions, etc, but future work will be needed to confirm the choice. We do not expect however that a better choice would change {\it qualitatively} the results presented here. The parameters used in the simulations are : a dissociation energy $D=0.04$ eV, a spatial scale factor of the Morse potential $\alpha = 4.45$ \AA$^{-1}$, a coupling constant $\beta=0.06$ eV/\AA, a mass $m=300$ a.m.u. The constant of the Nos\'e thermostat has been set to $M=1000$. A first scan of the dynamics of the model is obtained by imposing a slow temperature ramp (200K to 540K) that generates sets of states which are approximately equilibrated and are then used as initial states for simulations at constant temperature. Figure 14 shows a time evolution of the dynamics of the model at three temperatures. The stretching of the base pairs is indicated by a grey scale, darker dots corresponding to larger stretching. Looking at this figure, one notices immediately two major features. First, as one moves along an horizontal direction, i.e. along the molecule for a given time, the amplitude of the stretching varies very much from site to site. This is especially true at high temperature, but it is still noticeable at 150 K, well below the melting temperature. This shows that there is no equipartition of energy in this nonlinear system, but on the contrary a tendency for the energy to localize at some points, which is more and more pronounced as temperature increases. At high temperature, the figure shows large black regions which correspond to denaturated regions of the molecule. These black areas are the denaturation bubbles observed experimentally. At the highest temperature shown here (fig.14c) they extend over 20 to 50 base pairs and their boundaries are sharp. If the temperature is raised slightly above 540 K, the bubbles grow even more and finally extend over the whole chain: the molecule is completely denaturated. The second remarkable feature on fig.14 can be observed by moving along a vertical line on the figure i.e. following the time evolution of a given base pair. If one choses one region of the molecule in which the energy is concentrated, one can see alternating black and light-grey dots. This is due to an {\it internal breathing} of the localized excitations that oscillate between a large amplitude (black dots in the figure) and a small amplitude state (light dots) in a regular manner. These motions are the fluctuational openings of DNA discussed above. They exist even well below the denaturation temperature and coexist with denaturated bubbles in the high temperature range. Figure 14b shows that they play the role of precursor motions for the formation of the bubbles. The existence of these breathing modes shows that the breathing solutions obtained from nonlinear dynamics, discussed above, are likely to play some role in DNA since they form spontaneously. Numerical simulations can also reveal some interesting features of the process which can lead to the formation of large amplitude fluctuational openings in DNA. Simulations of the collisions of breathing modes [Dau5] show that discreteness effects are very important in the mechanism that leads to the formation of the bubbles. First, discreteness can stabilize existing bubbles. In the continuum limit breathers are not stable and radiate their energy as small amplitude phonon modes. On the contrary, in a discrete lattice, they can persist during a very long time without radiating. But, more importantly, discreteness plays a direct role in the growth of the breathers because, in the lattice, collisions between breathers are not elastic. When two excitations collide, they exchange energy, and the simulations show that, statistically, the larger excitations increase in energy at the expense of the smaller ones. This provides an efficient way for building large amplitude modes because, as soon as they start to grow, they become able to extract energy from the smaller modes with which they collide, which in turn increases their growth rate. Although this is still very tentative, one could speculate that a similar mechanism could explain the formation of the transcription bubble, which would be able to collect energy from the thermally excited modes of the molecule. \medskip \titleb{\tbn An improved model for denaturation} \smallskip The theoretical studies of denaturation provide data which can be compared with experiments because the thermal denaturation of DNA has been largely investigated experimentally. This comparison can therefore be used to test and improve the model. Although the PB model is able to describe some of the main features of DNA melting, and in particular the divergence of the mean stretching $\langle \psi \rangle$ and the formation of ``bubbles'', there is a crucial point in which this model gives incorrect results, it is the {\it sharpness} of the phase transition. For an homopolymer, the experiments show that the melting occurs very abruptly over a temperature interval which is only a few K or even less. This poses a very fundamental question since DNA is basically a one dimensional system, which is not expected to have a phase transition. It is interesting to notice that, {\it within a one-dimensional model with short range interactions}, a sharp transition seems to be possible if one takes into account properly the nonlinearity of the base stacking interaction [Dau3]. The possibility of a phase transition in one-dimensional DNA was already examined within the Ising-model approach by Poland and Scheraga [Pol] and Azbel [Azb2] who concluded that it can be attributed to cooperativity effects and to the role of the winding entropy released when the two strands separate. A simple extension of the DNA model defined by the hamiltonian (1) of the previous section can describe the dynamics of these effects an gives a very sharp transition in agreement with the experiments. The stacking energy between two neighboring base pairs must be described by the anharmonic potential: $$ W(\psi_n,\psi_{n-1})={\beta\over 2} \Bigl(1+\rho e^{-\delta(\psi_n+\psi_{n-1})} \Bigr)\ (\psi_n-\psi_{n-1})^{2} \; . \eqno(1) $$ This new intersite coupling, replacing the simple harmonic coupling of the previous approach, is responsible for qualitatively different properties. The choice of this potential has been motivated by the observation that the stacking energy is not a property of {\it individual} bases, but a character of the base {\it pairs} themselves. When the hydrogen bonds connecting the bases break, the electronic distribution on bases are modified, causing the stacking interaction with adjacent bases to decrease. And simultaneously the overlapping between the electronic distributions of adjacent bases decreases as the bases move from their equilibrium positions. In eq.(1), this effect is enforced by the prefactor of the usual quadratic term $(\psi_n-\psi_{n-1})^{2}$. This prefactor depends on the {\it sum} of the stretchings of the two interacting base pairs and decreases from ${1\over2}\beta(1+\rho)$ to ${1\over2}\beta$ when either one (or both) base pair is stretched. Although its expression was chosen for analytical convenience, the shape of potential (1) is in agreement with the properties of chemical bonds in DNA. They also provide the cooperativity effects that were introduced phenomenologically in the Ising models. A base pair that is in the vicinity of an open site has lower vibrational frequencies, which reduces its contribution to the free energy. Simultaneously a lower coupling along the strands gives the bases more freedom to move independently from each other, causing an entropy increase which drives a sharp transition. This phenomenon can be compared to recent views on structural phase transition in elastic media which stress that intrinsic nonlinear features characterize the physics of these transformations, and extend the standard soft mode picture [Kru4, Ker]. It is important to notice that, although cooperativity in introduced through purely {\it nearest neighbor} coupling terms, it has a remarkable effect on the 1D transition. Figure 15 shows the drastic change introduced by the anharmonic coupling on the specific heat of the model calculated by the transfer integral method [Dau3]. The full curve corresponding to the anharmonic stacking interaction shows a sharp peak very similar to that one would expect from a first-order phase transition, whereas the harmonic coupling investigated before gives only a smooth maximum. This result suggests that, although DNA structure is very complicated, a simple nonlinear model is able to reproduce with a good agreement the main experimental features of its dynamics for the fluctuational openings as well as the melting curves. Indeed such a model does not pretend to give exact quantitative results that would fit exactly the experimental melting curve of an heteropolymer with all its small structure, but it can perhaps provide a basis for the understanding of DNA dynamics and thermodynamics. \vfill \eject \titlea{CONCLUSIONS AND DISCUSSION} \medskip We have considered in some detail two specific models describing DNA nonlinear dynamics, the one proposed by Yakushevich [Yak2] to describe the motion of the transcription bubble along the proposal of Englander et al. [Eng], and the one proposed by Peyrard and Bishop to describe fluctuational openings; these are representatives of larger classes of models, but we have focused our attention on them as several qualitative and quantitative results are available, showing good agreement with experiment. These models are, as already stressed, extremely simple; so simple that actually only one degree of freedom per base is considered. This fact should not be seen as a drawback: quite on the contrary, the interest of such simple models resides mainly in this simplicity. Indeed, they point out what is the dominant degree of freedom for the process we want to study (i.e. fluctuational openings in the case of PB model, motion of the transcription bubble in the case of the Y model), and therefore provide a fundamental understanding. While adding other details, and maybe considering other degrees of freedom, could expectedly give a more refined description of DNA dynamics, the results reported here show that the dynamics of some fundamental biological processes is well described by very simple nonlinear models; this makes us confident that this approach is a way toward an understanding of the basic dynamical features of these processes, and that it gives a hint on the relevant degrees of freedom. Quite clearly, the very existence of a single dominant degree of freedom for such processes is by itself a nontrivial result, in view of the complexity of the DNA molecule. May-be the future is into slightly less simplified models in which a better description of the geometry would allow a more exact description of the helicoidal character of DNA. Here we have presented a first step to include this helicoidal aspect, but a more complete description is probably important: in an helicoidal structure, the bending of the molecule is strongly coupled to the local unwinding of the helix. Therefore one can expect that the bending would be strongly coupled to the opening. This seems to be true for the transcription bubble created by the RNA polymerase. One interesting aspect of the statistical mechanics is the way a big excitation can grow by collecting the energy of smaller ones when discreteness effects are strong enough. This observation, derived from numerical simulations, can perhaps provide a basis for the understanding of the formation of the transcription bubble in the presence of RNAP, without an external input of energy. Finally, we would like to comment on the mathematical tools used in this paper, or better in the study of the models considered here. Indeed, one could be quite surprised that we confidently reduce a complex dynamics to a very specific reduced form; it should be recalled, anyway, that in the PB model we are led to study a nonlinear Schroedinger equation (NLS), while in the Y model we are confronted with a (double) sine-Gordon equation (SG). Now, these are actually "universal equations", in the sense that under very general assumptions a large class of nonlinear field equations in one spatial dimension admit either the NLS either the SG equation as slow variation / small amplitude limit; this substantiates also the expectation that, whatever the more refined models we could consider, their prediction should not differ too much from those of the simple models we have been considering here, and that could be seen as phenomenological models. Another comment which is in order, especially considering that we aim at attracting physicists to this fascinating area, is that we have been using a rather restricted box of tools from nonlinear dynamics. Dynamical systems theory has undergone an impressive development in recent years, and it is to be expected that the available results from the general theory can provide precious insights into other, more intricate, features of DNA nonlinear dynamics. \vfill \eject \def\RA#1{[#1]} \def\RN#1{[[#1]]} \newcount\refnumber \def\clearrefnumber{\refnumber=0\relax} \def\RN{\advance\refnumber by 1 [{\the\refnumber}] } \parindent=0pt \parskip=5pt \titlea{References.} \RN \RA{Azb1} M.Ya. Azbel, {\it J. Chem. Phys.} {\bf 62} (1975), 3635 \RN \RA{Azb2} M.Ya. Azbel, {\it J. Phys. A} {\bf 12} (1979), L29 \RN \RA{Bal} E. Balanovski, {\it Int. J. Theor. Phys.} {\bf 26} (1987), 49 \RN \RA{Ban} A. Banerjee and H.M. Sobell, {\it J. Biomol. Struct. and Dyn.} {\bf 1} (1983), 253 \RN \RA{Cal1} C.R. Calladine, {\it J. Mol. Biol.} {\bf 161} (1982) 343 \RN \RA{Cal2} C.R. Calladine and H.R. Drew, {\it "Understanding DNA"}, Academic Press, London, 1992 \RN \RA{Chr} P.L. Christiansen and V. Muto, {\it Physica D} {\bf 68} (1993), 93 \RN \RA{Cur} J F. Currie, J A. Krumhansl, A.R. Bishop and S E. Trullinger, {\it Phys. Rev. B} {\bf 22} (1980), 477 \RN \RA{Dau1} T. Dauxois, {\it Phys. Lett. A} {\bf 159} (1991), 390 \RN \RA{Dau2} T. Dauxois, M. Peyrard and C.R. Willis, {\it Physica D} {\bf 57} (1992), 267 \RN \RA{Dau3} T. Dauxois, M. Peyrard and A.R. Bishop, {\it Phys. Rev. E} {\bf 47} (1993), R44 \RN \RA{Dau4} T. Dauxois, M. Peyrard and A.R. Bishop, {\it Phys. Rev. E} {\bf 47} (1993), 684 \RN \RA{Dau5} T. Dauxois and M. Peyrard, {\it Phys. Rev. Lett.} {\bf 70} (1993), 3935 \RN \RA{Dav} A.S. Davydov, {\it J. Theor. Biol.} {\bf 38} (1973), 559 \RN \RA{Dav} A.S. Davydov, {\it "Biology and Quantum Mechanics"}, Pergamon, New York, 1982 \RN \RA{Dic1} R.E. Dickerson, Scientific American, (december 1983), 87 \RN \RA{Dic2} R.E. Dickerson et al. {\it Nucl. Acids Res.} {\bf 17} (1989), 1797 \RN \RA{Dub} A. Dubrovin, S.P. Novikov and A. Fomenko, {\it "Geometrie Contemporaine I, II \& III"}, Mir, Moscow, 1982 \& 1987; {\it "Modern Geometry I \& II"}, Springer, New York, 1984 \RN \RA{Duv1} G. Duval-Valentine and R. Ehrlich, {\it Nucl. Acid Res.} {\bf 14} (1986), 1967 \RN \RA{Duv2} G. Duval-Valentine and R. Ehrlich, {\it Nucl. Acid Res.} {\bf 15} (1987), 575 \RN \RA{Eng} S.W. Englander, N.R. Kallenbach, A.J. Heeger, J.A. Krumhansl and S. Litwin, {\it Proc. Natl. Acad. Sci. (USA)} {\bf 777} (1980), 7222 \RN \RA{Fed} V.K. Fedyanin and L.V. Yakushevich, {\it Studia Biophys.} {\bf 103} (1984), 171 \RN \RA{For1} K. Forinash, A.R. Bishop and P.S. Lomdahl, {\it Phys. Rev. B} {\bf 43} (1991), 10743 \RN \RA{For2} K. Forinash and J. Keeney, {\it J. Biol. Phys.} {\bf 18} (1991), 19 \RN \RA{Gab1} J. Gabarro-Arpa and C. Reiss; in {\it "Progress in molecular and subcellular biology"}, F.E. Hahn ed. (vol.5, p.1), Springer, New York 1977 \RN \RA{Gab2} J. Gabarro-Arpa, P. Tougard and C. Reiss, {\it Nature} {\bf 280} (1979), 515 \RN \RA{Gae1} G. Gaeta, {\it Phys. Lett. A} {\bf 143} (1990), 227 \RN \RA{Gae2} G. Gaeta, {\it Phys. Lett. A} {\bf 168} (1992), 383 \RN \RA{Gae3} G. Gaeta, {\it Phys. Lett. A} {\bf 172} (1993), 365 \RN \RA{Gae4} G. Gaeta, {\it Phys. Lett. A} {\bf 179} (1993), 167 \RN \RA{Gae5} G. Gaeta, to appear in the Proceedings of Leningrad DNA conference (June 1992) \RN \RA{Gao} Y. Gao, K. V. Devi Prasad, E. W. Prohofsky, {\it J. Chem. Phys.} {\bf 80} (1984), 6291 \RN \RA{Gue} M. Gu\'eron, M. Kochoyan and J.L. Leroy, {\it Nature} {\bf 328} (1987), 89 \RN \RA{Heu} H. Heumann, M. Ricchetti and W. Werel, {\it EMBO J.} {\bf 7} (1988), 4379 \RN \RA{Hip} P. von Hippel, D.G. Baer, W.D. Morgan and J.A. McSwiggen, {\it Ann. Rev. Biochem.} {\bf 53} (1984), 389 \RN \RA{Hom} S. Homma and S. Takeno, {\it Prog. Theor. Phys.} {\bf 72} (1984), 679 \RN \RA{Jac1} M.A. Jacquet, R. Ehrlich and C. Reiss, {\it Nucl. Acid Res.} {\bf 17} (1989), 2933 \RN \RA{Jac2} M.A. Jacquet and C. Reiss, {\it Nucl. Acid Res.} {\bf 18} (1990), 1137 \RN \RA{Jac3} M.A. Jacquet and C. Reiss, {\it Mol. Microbiol.} {\bf 6} (1992), 1681 \RN \RA{Ker} W.C. Kerr, A.M. Hawthorne, R.J. Gooding, A.R. Bishop and J.A. Krumhansl, {\it Phys. Rev. B} {\bf 45} (1992), 7036 \RN \RA{Kru1} J.A. Krumhansl and J.R. Schrieffer,{\it Phys. Rev. B} {\bf 11} (1975), 3535 \RN \RA{Kru2} J.A. Krumhansl and D.M. Alexander; in {\it Structure and dynamics: Nucleic acids and proteins}, Eds. E. Clementi and R.H. Sarma, Adenine Press, New York, 1983 \RN \RA{Kru3} J.A. Krumhansl, J.A. Wysin, D.M. Alexander, A. Garcia, P.S. Lomdahl, A.C. Scott and P. Layne; in {\it Structure and Motion: Membranes, Nucleic Acids and Proteins}, Clementi, E., Corongiu, G., Sarma, M.H. and Sarma, R.H. editors, Adenine Press, New York, 1985 \RN \RA{Kru4} J. A. Krumhansl and R. J. Gooding, {\it Phys. Rev. B} {\bf 39} (1989), 3047 \RN \RA{Lad} J. Ladik and J. Cizek, {\it Int. J. Quant. Chem.} {\bf XXVI} (1984), 67 \RN \RA{Lew} M. Lewitt, {\it Cold Spring Harbor Symp. Quant. Biol.} {\bf 46A} (1982), 251 \RN \RA{Mar1} H. Marcaud, J. Gabarro-Arpa, R. Ehrlich and C. Reiss, {\it Nucl. Acid Res.} {\bf 14} (1986), 551 \RN \RA{Mar2} F. Marchesoni and C. Lucheroni, {\it Phys. Rev. B} {\bf 44} (1991), 5303 \RN \RA{Mut1} V. Muto, J. Halding, P.L. Christiansen and A.C. Scott, {\it J. Biomol. Struct. Dyn.} {\bf 5} (1988), 837 \RN \RA{Mut2} V. Muto, J. Halding, P.L. Christiansen and A.C. Scott, {\it J. Biomol. Struct. Dyn.} {\bf 5} (1988), 873 \RN \RA{Mut3} V. Muto, P.S. Lomdahl and P.L. Christiansen, {\it Phys. Lett. A} {\bf 136} (1989), 33 \RN \RA{Mut4} V. Muto, A.C. Scott and P.L. Christiansen, {\it Physica D} {\bf 44} (1990), 75 \RN \RA{Mut5} V. Muto, P.S. Lomdahl and P.L. Christiansen, {\it Phys. Rev. A} {\bf 42} (1990), 7452 \RN \RA{Nos} S. Nose, {\it Mol. Phys.} {\bf 52} (1984), 255; and {\it J. Chem. Phys.} {\bf 81} (1984), 511 \RN \RA{Pey1} M. Peyrard and A.R. Bishop, {\it Phys. Rev. Lett.} {\bf 62} (1989), 2755 \RN \RA{Pey2} M. Peyrard, T. Dauxois, H. Hoyet and C.R. Willis, {\it Physica D} {\bf 68} (1993), 104 \RN \RA{Pol} D. Poland, H. Scheraga,{\it J. Chem. Phys.} {\bf 81} (1984), 511 \RN \RA{Pro1} E.W. Prohofsky,{\it Phys. Rev. A} {\bf 38} (1988), 1538 \RN \RA{Pro2} E.W. Prohofsky,{\it Phys. Rev. A} {\bf 38} (1988), 5332 \RN \RA{Sae} W. Saenger, {\it "Principles of Nucleic Acid Structure"}, Springer, New York, 1984; second corrected printing 1988 \RN \RA{Sal} M. Salerno, {\it Phys. Rev. A} {\bf 44} (1991), 5292 \RN \RA{Sca} D. J. Scalapino, M. Sears, R. A. Ferrel, {\it Phys. Rev. B} {\bf 6} (1972), 3409 \RN \RA{Sch} B. Schmitt and C. Reiss, {\it Biochem. J.} {\bf 277} (1991), 435 \RN \RA{Sco1} A.C. Scott, {\it Am. J. Phys.} {\bf 37} (1969), 52 \RN \RA{Sco2} A.C. Scott, F.V.F. Chu and D.W. Mac Laughlin, {\it Proc. IEEE} {\bf 61} (1980), 1443 \RN \RA{Sco3} A.C. Scott, {\it Phys. Rev. A} {\bf 31} (1985), 3518 \RN \RA{Sob1} H.M. Sobell, C.C. Tsai, S.G. Gilbert, S.C. Jain and T.D. Sacore, {\it Proc. Natl. Acad. Sci. USA} {\bf 73} (1976), 3068 \RN \RA{Sob2} H.M. Sobell, C.C. Tsai, S.C. Jain and S.G. Gilbert, {\it J. Mol. Biol.} {\bf 114} (1977), 333 \RN \RA{Sob3} H.M. Sobell, {\it Proc. Natl. Acad. Sci. USA} {\bf 82} (1985), 5328 \RN \RA{Tak1} S. Takeno and S. Homma, {\it Prog. Theor. Phys.} {\bf 70} (1983), 308 \RN \RA{Tak2} S. Takeno and S. Homma, {\it Prog. Theor. Phys.} {\bf 77} (1987), 548 \RN \RA{Tec1} M. Techera, L.L. Daemen and E.W. Prohofsky, {\it Phys. Rev. A} {\bf 40} (1989), 6636 \RN \RA{Tec2} M. Techera, L.L. Daemen and E.W. Prohofsky, {\it Phys. Rev. A} {\bf 42} (1990), 1008 \RN \RA{The} N. Theodorakopoulos and M. Peyrard, unpublished. \RN \RA{Vin} S.V. Vinogradov and R.H. Linnell, {\it "Hydrogen Bonding"}, Van Nostrand Reinhold, New York 1971 \RN \RA{Wan1} F.J. Wang and C. Reiss, {\it Biochem. and Mol. Biol. Int.} {\bf 30} (1993), 983 \RN \RA{Wan2} F.J. Wang, Thesis, Universit\'e Paris VI 1991 \RN \RA{Xia} J. Xiao, J. Lin and G. Zhang, {\it J. Phys. A.} {\bf 20} (1987), 2425 \RN \RA{Yak1} L.V. Yakushevich, {\it Studia Biophys.} {\bf 121} (1987), 201 \RN \RA{Yak2} L.V. Yakushevich, {\it Phys. Lett. A} {\bf 136} (1989), 413 \RN \RA{Yom1} S. Yomosa, {\it J. Phys. Soc. Jap.} {\bf 51} (1982), 3318 \RN \RA{Yom2} S. Yomosa, {\it Phys. Rev. A} {\bf 27} (1983), 2120 \RN \RA{Yom3} S. Yomosa, {\it Phys. Rev. A} {\bf 30} (1984), 474 \RN \RA{Zan1} L.L. Van Zandt, {\it Phys. Rev. A} {\bf 40} (1989), 6134 \RN \RA{Zan2} L.L. Van Zandt, {\it Phys. Rev. A} {\bf 42} (1990), 5036 \RN \RA{Zha} C.T. Zhang, {\it Phys. Rev. A} {\bf 35} (1987), 886 \RN \RA{Zho} G.F. Zhou and C.T. Zhang, {\it Physica Scripta} {\bf 43} (1991), 347 \vfill \eject \titleb{Figure captions} \footnote{}{{\it After each figure caption, we indicate the section in which the figure is first referred to.}} {\bf Figure 1.} Schematic drawing of the double helix (B-form). the sugar-phosphate backbone is represented as a ribbon, with detailed chemical structure shown in figure 1a and figure 1b (sugar moiety). The bases are sketched as rectangular plates, with detailed chemical structures shown in figure 1c. {\it (Section II-1)} {\bf Figure 2.} Purine (Pu) and pyrimidine (Py) sketched as large and small rectangular plates respectively. The table depicts the main intra-pair and inter-pair, rotational and traslational degrees of freedom and their names. {\it (Section II-1)} {\bf Figure 3.} Steric inter-base hindrance and how it can possibly be relieved: $(a)$ the purines are in closz contact (hatched area) in the major groove (see figure 2); $(b)$ relief provided by reducing propeller-twist in the upper base-pair rotation of the purine around $Oy$ (see figure 2); $(c)$ additional relief provided by enhanced twist of the upper base pair, relative to its lower neighbour. {\it (Section II-2)} {\bf Figure 4.} The rotational angles of the sugar-phosphate backbone, with a thymine attached. Figures 4a and 4b depict the observed ranges of values of $\chi$ ({\it anti} and {\it syn} ranges), with the T base plane shown perpendicular to the "average" sugar plane. Figures 4c, 4d, 4e show the position of the O5' atom over the sugar, for the three observed value ranges of the torsional angle $\gamma$, i.e. trans-left, left-trans and left-left. {\it (Section II-3)} {\bf Figure 5.} Sugar pucker. In figures 5a and 5b, only the twisted forms are shown (the observer is in the plane defined by C1', O4' and C4'), in th emost common conformations, i.e. C3' endo (figure 5a) and C2' endo (figure 5b). The effect of C2' endo and C3' endo conformations on the distance between two consecutive phosphates is emphasized. {\it (Section II-3)} {\bf Figure 6.} Nucleoside free energy difference as a function of the pseudo-rotational angle $P$. The ribose pucker corresponding to selected values are displayed on top. Notice the shallow energy barrier between the favoured C3' endo and C2' endo puckers. {\it (Section II-3)} {\bf Figure 7.} Thermal melting profile of the linearized $5.5$ kbase long DNA from the simian virus SV40 in $19$mM ${\rm Na}^+$. The variation of the absorbance at $260$nm with temperature $T$ is plotted versus $T$. {\it (Section II-5)} {\bf Figure 8.} Interaction between pairing bases in the Yakushevich model. Bases are considered as rigid and free only to rotate around a fixed axis (compare with the original proposal of DNA solitons in [Eng]); moreover one usually considers $d \to 0$. Here the hydrogen bond pairing between the bases is considered to give origin to a restoring elastic force, depending linearly on the distance $\ell$ between the atoms on the two bases bridged by the H bond; the nonlinearity appears when considering such a force as a function of the two torsion angles $\theta_1 , \theta_2$. {\it (Section IV-1)} {\bf Figure 9.} Abstract representation of planar models of DNA. Each base is considered to interact only with first neighbours on its chain (stacking interactions, modelled by harmonic potential), and on the pairing base at the same site on the other chain (pairing forces due to H bonds, modelled by non-harmonic potentials). The abstract model can be considered planar, which expalins the name; in actual DNA, the geometry of the backbone is at the origin of the double helix structure, with a turn of about $36 {}^0$ per base. {\it (Section V-1)} {\bf Figure 10.} Abstract representation of helicoidal models of DNA. In addition to stacking and pairing interactions indicated in figure 9, these consider also interaction between bases on opposite chains which are half-pitch of the helix away in the sequence, and therefore nearby in three-space due to helical geometry; two bases connected by such an "helicoidal" interaction are indicated. In real DNA, these interaction are mediated by water filaments joining the concerned bases. {\it (Section V-1)} {\bf Figure 11.} Dispersion relations for the helicoidal Peyrard- Bishop model, see equation(V-3.4). The functional form of the dispersion relations is also the same for the helicoidal Yakushevich model, bu the values (and the meaning) of the parameters appearing in these is different, see section V-3. Here we use the value of parameters appropriate for the hPB model as given in [Dau1, Dau4], see Table I. {\it (Section V-3)} {\bf Figure 12.} Numerical results for the energy of the (1,1) soliton as a function of