\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 $\nu \equiv v / v_{{\rm max}}$, with the
values of parameters given in the text. Minimization of the Lagrangian
is performed by a Montecarlo algorithm for 35 values of $\nu$; the
corresponding value of $H (\nu )$ is plotted by solid dots. The
continuous line represents the best fit with functions of the form $H
(\nu ) = A / (B - \nu^2 )$. {\it (Section V-8)}
{\bf Figure 13.} Variation of $\langle \psi \rangle $ versus
temperature: the dash line corresponds to the TI results in the
continuum limit, the solid line gives the exact TI results obtained by
numerical solution of the TI operator, and the plus signs correspond
to molecular dynamics simulations. {\it (Section VI-1)}
{\bf Figure 14.} Results of molecular dynamics simulations at three
different temperatures (a) $T=150$ K, (b) $T = 340$ K, (c) $T=450$ K.
The horizontal axis indicates the position along the 256 cells of the
molecule and the vertical axis indicates time. The stretching $\psi_n$
of the base pairs along the molecule is indicated by a grey scale, the
lighter grey corresponding to $ \psi \leq -0.1$ \AA\ and black
indicating $\psi \geq 1$ \AA. Therefore black regions show broken
base-pairs. {\it (Section VI-2)}
{\bf Figure 15.} Variation of the specific heat versus temperature for
the PB model and its extension with anharmonic coupling. The very
narrow peak corresponds to the anharmonic coupling case
$(\delta=0.35,\rho=0.5)$, the dotted curve and the solid broad peak to
harmonic coupling with coupling constants $\beta' = 1.5\beta$
and $\beta''=\beta$, respectively, where $\beta$ is the coupling constant
used for the numerical simulations discussed in section VI-2. {\it
(Section VI-3)}
\bye