\section{Introduction}
Understanding of the detailed molecular mechanisms of conformational
dynamics of the double helical DNA is an important goal on the long
way to control of genom functions. Computer simulations are the main
theoretical instrument in such studies. An important recent progress
in the methodology \cite{AMBER94:,MacKerell:95,Darden:93} has made
possible molecular dynamics (MD) simulations of short DNA fragments
in realistic water environment with explicit counterions
(comprehensive surveys of the literature can be found in the recent reviews
\cite{Ravishanker:97,Auffinger:98}). In should be noted, however,
that although in future these methods will
probably be able to describe in full detail
all interactions involved in the numerous biological functions of DNA,
at present their capabilities are very limited. There are many
important domains where they can hardly be applied, in the foreseeable
future, for instance, dynamics of long linear DNA fragments or
plasmids. Even for small oligomers these methods
require supercomputer resources because the solute
molecule must be placed in a large enough water box to accommodate all
neutralizing counterions without exceeding reasonable levels of
effective DNA and salt concentrations. This makes any methodological
questions difficult, especially those concerning statistical
reproducibility of results because, until now, almost all reports on
the subject described observations made for single trajectories.

The foregoing arguments explain why it is necessary to find
alternative conditions for MD simulations of B-DNA which would be less
computationally demanding, but yet acceptably accurate. The utility of
such DNA model would be twofold. On the one hand, it can give a
unique opportunity to probe the properties of large systems of real
biological importance where more rigorous approaches would be
prohibitively expensive. On the other hand, for smaller molecules they
can provide for an exhaustive sampling of the configurational space,
which is a necessary prerequisite of a more systematic studies of detailed
mechanisms involved in the conformational dynamics of double helical
DNA. We have recently found that surprisingly stable dynamics of B-DNA
duplexes can be obtained with semi-explicit treatment of long-range
electrostatics effects. \cite{Mzjacs:98} Paradoxically, it appeared
that this simplistic approach results in B-DNA conformations that
are significantly closer to experimental data than ever before,
including the recent more rigorous and expensive calculations. This
relative success raised several intriguing questions concerning its
origin. On the one hand, a better agreement with experimental data can
always result from a fortunate cancelation of errors. On the other
hand, it is possible that the earlier reported {\em a priori} more
rigorous calculations failed to reveal the full accuracy of the recent
atom parameter sets, notably AMBER94, \cite{AMBER94:} because of
insufficient duration of trajectories or some
problem in calculation of long range electrostatic interactions
by the particle-mesh Ewald technique. \cite{Darden:93}
Both these possibilities should be checked, but the
former is certainly easier to analyze because this does not require
expensive computations.

In this paper we address some of the above questions by further
exploring the dynamic behavior of B-DNA duplexes with partial explicit
hydration shell. Our main objective is to justify what we call a
``minimal B-DNA'', that is the least computationally demanding
simulation conditions where the essential experimental features of
this molecule can be reliably reproduced. It is shown that a good
candidate is a B-DNA duplex with explicitly hydrated minor groove,
empty major groove, and the electrostatic effects of the bulk solvent
accounted for implicitly. We report also about relatively successful attempts
to obtain stable dynamics without explicit hydration, which, however,
give much worse agreement with experimental data. In contrast,
dynamics of the minimal model reliably reproduces certain well-known sequence
dependent effects, notably, modulation of the minor groove width and
spontaneous formation of specific ``spine'' hydration pattern observed
in many experimental structures.  At the same time, for two very
different sequences it converges to significantly different
conformations which, however, are both close to the canonical B-DNA
form.  These results appear relatively insensitive to variations of
parameters involved in the simplified treatment of long-range
electrostatic interactions.

\section{Methods and Simulation Protocols}
All new MD simulations reported here have been performed with
the internal coordinate molecular dynamics (ICMD) method.
\cite{Mzjbsd:89,Jain:93,Mzjcc:97,Mzprep:98,Mzjchp:99}
Here we employ a modified implicit leapfrog integrator detailed in the Appendix.
We consider dynamics of two different DNA duplexes with
sequences d(CGCGAATTCGCG)$_2$ and d(GCGCGCGCGCGC)$_2$.
Most of the calculations were made with the first molecule often
referred to in the literature as ``{\em Eco}RI dodecamer''.
The second dodecamer is used for comparison.
The {\em Eco}RI sequence was the first to
crystallize in the B-form of DNA \cite{Drew:81} and since then it
became perhaps the most studied DNA fragment both experimentally and
theoretically.
\cite{Levitt:83b,Rao:90,Srinvasan:90,Swaminathan:91,Miaskiewicz:93,Kumar:94,McConnell:94}
In the recent years it has been often used in benchmark tests of new
force fields and algorithms.
\cite{Mzjacs:98,York:95,Yang:96,MacKerell:97,Duan:97,Young:97,Young:97a,Cieplak:97,Sprous:98,Jayaram:98,Winger:98}
The DNA model was same as in our recent report,
\cite{Mzjacs:98} namely, all torsions were free as
well as bond angles centered at sugar atoms, while other bonds and
angles were fixed, and the bases held rigid. The dynamic properties of
this model, notably, the fastest motions and maximal possible time
steps have already been studied. \cite{Mzjacs:98} AMBER94
\cite{AMBER94:} force field and atom parameters were
used with TIP3P water \cite{TIP3P:} and no cut off schemes.

The solvent effects were taken into account by the
mixed strategy \cite{Mzjacs:98} which we continue to explore and
improve here. In this method, following to a long known approximate
approach, \cite{Tidor:83,Mazur:91,JUMNA:95}
long range effects are taken into account
implicitly by reducing phosphate charges and using a linear distance
dependent dielectric function $\varepsilon = r$. On the other hand,
the minor groove of the DNA molecule is filled up by explicit water
with the water cloud somewhat protruding in space. In the initial
hydration procedure the DNA molecule is first covered by a 5~\AA\
thick water shell by placing it in a water box and eliminating
overlapping and distant solvent molecules. To allow water molecules to
penetrate into the minor groove, rather small cut-off distances are
used as the overlap criteria, namely, 1.8 and 1.3~\AA\ for oxygen and
hydrogen water atoms, respectively. After that cylinder-like volumes
around each strand are built from spheres centered at phosphorus atoms
with radii of 12~\AA. All water molecules that appear outside the
intersection area of the two volumes are removed. The solvent
remaining is next relaxed by energy minimization first with the solute
held rigid and then with all degrees of freedom. This procedure gives
a partially hydrated duplex, with the minor groove completely filled
and a few solvent molecules in the major groove. In some calculations
we additionally cleaned the major grove from the remaining water after
the second minimization. The motivations for this additional step are
discussed in the text.

The standard heating and equilibration protocols have been slightly
modified compared with our previous studies. Dynamics was initiated
with zero solvent temperature by giving the solute a kinetic energy
corresponding to 300~K. During the initial short 2 ps run the
temperature was weakly coupled to 250~K by the Berendsen algorithm
\cite{Berendsen:84} with a relaxation time of 1 ps. During the next 5
ps the temperature coupling was switched off and at the end of the
whole 7 ps period the system temperature normally reached 200~K.
Starting from this state the production trajectory was computed with
the temperature coupled to 300~K with a relaxation time of 10 ps. This
algorithm provides for a softer start and helps to reduce dissociation
of hot water molecules during the early non-equilibrated phases of
dynamics. \cite{Mzjacs:98} For trajectories starting
from the canonical B-DNA form the
energy minimization was applied both before and after hydration.
Hydration of the Xray structure \cite{Drew:81} was performed with
all crystallographic water molecules kept in their places.

Normally we used time steps of 5 and 10 fs in heating and production
runs, respectively. The possibility and the necessary conditions for
this large time steps have been discussed elsewhere.
\cite{Mzjacs:98,Mzjcp:97,Mzjpc:98}
Duration of trajectories varied depending upon
specific conditions discussed in the text. In calculations with the
solvent shell a minimum duration of 5 ns was usually sufficient to
arrive at stable levels of average helical parameters as well as
rmsd's from the reference conformations. In production runs
conformations were saved with a 2.5 ps interval. Structures from the
last nanosecond were used for computing an average conformation
referred to as the final MD state of the corresponding trajectory. The
reference canonical A and B forms were constructed with NUCGEN utility
of AMBER \cite{AMBER:} or JUMNA program. \cite{JUMNA:95} Curves
\cite{Curves:} and ``Dials and Windows'' procedures \cite{CDW:} were
used for analysis of DNA conformations. XmMol program \cite{XmMol:}
was used for visual inspection and animation of trajectories.

\begin{table*}[t]\caption{\label{Ttjli} General Information about
Trajectories}
\begin{tabular}[t]{|ddddddd|}
Code & Sequence & From  & Water\tablenotemark[1] & Clean\tablenotemark[2]
& Q$_p$\tablenotemark[3] (aeu) & Duration (ns) \\\hline
TJA\tablenotemark[4]
    &CGCGAATTCGCG& B73\tablenotemark[5]   & 134  &   -   &   -0.5   &  5.0 \\
TJB\tablenotemark[4]
    &CGCGAATTCGCG& {\em Eco}RI\tablenotemark[6]
                                          & 114  &   -   &   -0.5   &  5.1 \\
TJC\tablenotemark[4]
    &CGCGAATTCGCG& B73                    & 134  &   -   &   -0.5   &  5.0 \\
TJD &CGCGAATTCGCG& Xray\tablenotemark[7]  & 131  &   -   &   -0.5   &  8.8 \\
TJE &CGCGAATTCGCG& B73                    &  -   &   -   &   -0.25  &  2.7 \\
TJF &CGCGAATTCGCG& B73                    &  -   &   -   &   -0.5   &  5.0 \\
TJG &CGCGAATTCGCG& B73                    &  -   &   -   &   -0.75  &  2.9 \\
TJH &CGCGAATTCGCG& Xray                   & 175  &   +   &   -0.25  &  5.0 \\
TJI &CGCGAATTCGCG& Xray                   & 175  &   +   &   -0.5   &  5.5 \\
TJJ &CGCGAATTCGCG& Xray                   & 175  &   +   &   -0.6   &  8.0 \\
TJK &CGCGAATTCGCG& Xray                   & 175  &   +   &   -0.75  &  6.9 \\
TJL &GCGCGCGCGCGC& B73                    & 166  &   +   &   -0.5   &  8.1 \\
\end{tabular}
\tablenotetext[1]{
The total number of explicit water molecules
left after the hydration procedures described in the text.}
\tablenotetext[2]{
Major groove water removed before equilibration.}
\tablenotetext[3]{
Reduced total phosphate charge.}
\tablenotetext[4]{
Trajectories reported and discussed in detail in Ref. \onlinecite{Mzjacs:98}
}
\tablenotetext[5]{
The canonical B-DNA form with Arnott B73 parameters. \cite{Arnott:72} }
\tablenotetext[6]{
The `kinked' duplex conformation found in the {\em Eco}RI endonuclease complex.
\cite{McClarin:86} (file PDE001 in Nucleic Acids Database \cite{NDB:}) }
\tablenotetext[7]{
The crystal conformation. \cite{Drew:81}
(file 1bna in Protein Database. \cite{PDB:}) }
\end{table*}

\section{Results and Discussion}
\subsection{General Outline of Simulations and Results}
Tables \ref{Ttjli}, \ref{Thlpa} and Figs. \ref{Falstr} and \ref{Falmg}
characterize in general the numerical experiments and results
considered in detail in sections below. Table \ref{Ttjli} summarizes
variations of computational protocols for different trajectories. For
brevity each trajectory is referred to by using a single-letter code
with a prefix ``TJ''. Trajectories from A to K computed for the
{\em Eco}RI dodecamer differ either by the starting state, or by the
details in the application of hydration procedures, or by the reduced
phosphate charge. TJA, TJB, and TJC have been the subject of our
previous report. \cite{Mzjacs:98} They are discussed together with the
rest and are included in Table \ref{Ttjli} for completeness. It is
worth reminding here that TJC was computed by using the conventional
Cartesian coordinate MD method starting from the same state as TJA.
TJD and TJI also started from the same conformation, namely, the Xray
structure, \cite{Drew:81} with
significantly different initial distribution of water
molecules obtained by varying parameters of the hydration algorithm.
TJE, TJF, and TJG were computed in the dehydrated state of the duplex
to clarify the role of the explicit water and study the effect of
phosphate screening. Similarly, TJH, TJI, TJJ, and TJK were calculated
for the same model system with different phosphate screening
parameters. Finally, TJL computed for another sequence checks the
possibility of observing sequence dependent structural effects with
this DNA model.

\begin{table}[h]\caption{\label{Thlpa} Selected Helical Parameters of
Different DNA Conformations}
\begin{tabular}[t]{|dddddd|}
         & Xdisp & Inclin & Propel & Rise & Twist\\\hline
A\tablenotemark[1]
         & -5.43 &  19.12 &  13.70 & 2.56 & 32.70\\
B\tablenotemark[1]
         & -0.70 &  -5.98 &   3.87 & 3.38 & 36.01\\
Xray\tablenotemark[2]
         & -0.54 &   0.19 & -13.45 & 3.37 & 35.89\\
TJA      & -1.10 &   2.06 & -12.54 & 3.31 & 35.67\\
TJB      & -0.95 &   0.62 & -12.80 & 3.37 & 36.36\\
TJC      & -1.40 &   4.04 & -16.30 & 3.36 & 34.69\\
TJD      & -1.41 &   1.83 & -12.55 & 3.34 & 35.82\\
TJE      & -3.28 &   2.98 &  -5.62 & 3.41 & 31.13\\
TJF      & -3.14 &   2.66 &  -5.16 & 3.44 & 31.65\\
TJG      & -3.08 &   1.73 &  -3.96 & 3.49 & 31.79\\
TJH      & -2.05 &   3.94 &  -7.39 & 3.40 & 32.08\\
TJI      & -1.65 &   2.43 & -10.21 & 3.39 & 34.15\\
TJJ      & -1.24 &   -.32 & -11.50 & 3.37 & 34.84\\
TJK      & -0.17 & -11.91 & -10.54 & 3.51 & 34.89\\
TJL      & -0.37 &  -2.24 &  -4.53 & 3.42 & 34.56\\
\end{tabular}
\tablenotetext[1]{
Canonical DNA forms \cite{Arnott:72} constructed with
NUCGEN procedure of AMBER \cite{AMBER:} for the {\em Eco}RI dodecamer. }
\tablenotetext[2]{
The crystal conformation of the {\em Eco}RI dodecamer. \cite{Drew:81}
(file 1bna in Protein Database. \cite{PDB:}) }
\end{table}

\begin{figure}
\centerline{\psfig{figure=Falstr.ps,height=13cm,angle=0.,%
bbllx=300bp,bblly=80bp,bburx=560bp,bbury=720bp,clip=t}}
\caption{\label{Falstr}
The final states of representative trajectories in Table
\protect\ref{Ttjli}, together with the canonical B-form
and the Xray structure of the {\em Eco}RI dodecamer.
The optimal helical axis \protect\cite{Curves:} is shown by continuous
roughly vertical lines. On the right of each structure
the rmsd's from the canonical B-DNA and the Xray conformation
\protect\cite{Drew:81} are given.}
\end{figure}

\begin{figure}
\centerline{\psfig{figure=Falmg.ps,width=7cm,angle=0.,%
bbllx=120bp,bblly=200bp,bburx=480bp,bbury=570bp,clip=t}}
\caption{\label{Falmg}
The profiles of the minor grooves for the structures
shown in Fig. \protect\ref{Falstr}. The local minor groove width
is given in angstr{\"o}ms as computed with
the spline algorithm implemented in Curves.
\protect\cite{Stofer:94,Curves:} The dotted line indicates the
level corresponding to the canonical B-DNA.
}\end{figure}

Figure \ref{Falstr} presents the final states of representative
trajectories in Table \ref{Ttjli} together with the canonical B-form
and the Xray structure of the {\em Eco}RI dodecamer. The drawings are the
outputs of Curves program \cite{Curves:} which also computes the
optimal helical axis shown by a continuous roughly vertical line in
each plate. Figure \ref{Falmg} shows variations of the minor groove
widths for the same structures computed by the spline algorithm
implemented in Curves. \cite{Stofer:94} On the right of each
structure in Fig. \ref{Falstr} the rmsd from the canonical B-DNA is
given. For {\em Eco}RI dodecamer the second value shows the rmsd from
the Xray conformation. This structure is asymmetrical and, since the
sequence is palindromic, the rmsd can be computed in two orientations.
Both orientations were tried and the lower rmsd values are given in
Fig. \ref{Falstr}. It is seen that the structures are generally close
to the two experimental B-DNA conformations, with a few exceptions to
be discussed hereafter. The majority of the computed structures of the
{\em Eco}RI dodecamer have a characteristic profile of the minor groove with
a narrowing in the center, which is a well-known feature of the
experimental crystallographic structure. \cite{Kopka:83}
These results confirm our
initial report on these simulation conditions. \cite{Mzjacs:98} Time
dependencies of rmsd's are not analyzed here systematically. We note
only that for trajectories with the final rmsd values below 2~\AA\
they were qualitatively similar to those reported for TJA, TJB and
TJC. Namely, after the first 100 ps the rmsd from
the canonical B-form and the Xray structure are usually found between
2.5 and 3.5~\AA. In our earlier calculations \cite{Mzjacs:98} this
large initial shift was either already present in the initial
conformation or resulted from {\em in vacuo} energy minimization of
the canonical B-DNA form, but it also occurs in trajectories
starting from the Xray conformation. During the first two nanoseconds
the rmsd values gradually decrease, with the 2~\AA\ level commonly
passed within the second nanosecond. The lowest rmsd values from
experimental B-DNA conformations are usually reached within the first
three nanoseconds and next the drift becomes below the level of
fluctuations in the nanosecond time scale.

Selected averaged helical parameters corresponding to the structures
in Fig. \ref{Falstr} are assembled in Table \ref{Thlpa}. The complete
set of these parameters have been discussed for TJA, TJB and TJC,
\cite{Mzjacs:98} and here we note only that parameters not included in
Table \ref{Thlpa} do not show systematic deviations from experimental
values. In contrast, Xdisp and inclination obtained in different
simulations for {\em Eco}RI sequence show a small but systematic deviation.
Propeller is generally well reproduced, and it is included in Table
\ref{Falstr} because it plays an important structural role for AT-rich
sequences \cite{Coll:87} and is rather sensitive
to simulation
conditions. Finally, twist and rise are the key helical parameters
responsible for the major part of rmsd's. The time variation of the
helical parameters in the present B-DNA model has already been
characterized. \cite{Mzjacs:98} It is similar to the above described
behavior of rmsd's in the sense that it usually takes about 3 ns for
the main parameters to come close to the standard B-DNA values. In all
calculations reported here the initial large rmsd's are accompanied by
low average twist which usually raise from 32--33$^{\rm\circ}$\  to
34--35$^{\rm\circ}$\  during calculations.

The data presented in Figs. \ref{Falstr} and \ref{Falmg}, and Tables
\ref{Ttjli} and \ref{Thlpa} are discussed in detail in sections below.

\subsection{B-DNA Dynamics in the Dehydrated State}
It is known since long ago that reasonably good duplex DNA
conformations can be obtained by minimization of potential energy with
implicit modeling of solvent effects, notably, by reducing phosphate
charges and scaling of Coulomb forces with a distance dependent
electrostatic function. \cite{Lavery:94} On the other hand, it is also
known that these structures are not dynamically stable. \cite{Swaminathan:91}
The dehydrated state presents non-negligible theoretical interest because
the ``naked DNA'' is the core of the whole system and understanding
its mechanics is important for distinguishing between the trends that
originate from intra-duplex interactions and water effects. We briefly
report here about a few interesting lessons that we have learned from
ICMD simulations of B-DNA in the dehydrated state.

We hoped initially that, due to the reduced flexibility of the
standard geometry DNA model, its structure could be more stable than
in the earlier Cartesian coordinate MD simulations. It was found,
however, that trajectories generally tend to be trapped in earlier
described deformed states characterized by a collapsed minor groove.
\cite{Swaminathan:91} Such transformations invariably start from a phosphate $\rm
B_I\to B_{II}$ flip, which can next result in a significant deformation
of the standard B-DNA shape, with partial local closing of the minor
groove. This behavior seems to be an intrinsic defect of {\em in
vacuo} conditions. However, the probability of collapsing transitions
was high only during the first nanosecond when, in animation, the DNA
molecule exhibited very violent periodical low frequency bending and
stretching motions. This intermediate phase is caused by the heating
procedures that commonly last less then 10 ps and always use some
velocity scaling which, however, affects only normal modes that have
non-zero instantaneous velocities. The slowest normal modes of a
dodecamer B-DNA fragment have periods beyond 10 ps, \cite{Lin:97}
therefore, during
heating, those that happen to be in their fast phase are strongly
overheated. {\em In vacuo} the uneven distribution of the kinetic
energy is rather persistent and it apparently takes hundreds of
picoseconds to reach equipartition. In many test calculations we
observed that if the structure managed to survive the first
nanosecond, it continued to stay in a reasonably good B-DNA form for
rather long time so that its average dynamic properties could well be
characterized. $\rm B_I\to B_{II}$ transitions still occur, but the
minor groove does not close immediately and many such transitions
appear reversible.

Trajectories TJE, TJF, and TJG in Table \ref{Ttjli} have been computed
with special precautions that take into account the foregoing
qualitative features. Dynamics were initiated with random velocities
and during the first nanosecond the temperature was raised to 250~K by
the Berendsen algorithm \cite{Berendsen:84} with a relaxation time of
100~ps. Calculations were next continued with the standard
parameters, that is with the bath temperature of 300~K and a 10~ps
relaxation time. If a transition to a collapsed state was observed
the calculations were restarted from one of the preceding
conformations with random velocities at 300~K. Such restarts were
repeated until a sufficiently long continuous trajectory have been
obtained. In the case of TJF one restart was necessary to get a
continuous 5~ns trajectory. TJG was also stable during 5~ns after one
restart, but only the first 2.9 ns were used for analysis because a
few reversible $\rm B_I\to B_{II}$ transitions that happened in its
second part caused non-negligible temporary distortions. In contrast,
in the case of TJG three restarts were necessary to obtain a 2.7~ns
continuous run at the end of which the structure collapsed.

\begin{figure}
\centerline{\psfig{figure=Fnowt.ps,width=7cm,angle=0.,%
bbllx=180bp,bblly=300bp,bburx=430bp,bbury=470bp,clip=t}}
\caption{\label{Fnowt}
The time variations of rmsd from several
reference structure for TJF. A, B, Ini and Fin correspond to
the canonical A and B forms, the initial minimum energy structure,
and the final MD state, respectively. All rmsd values are
given in angstr{\"o}ms.
}\end{figure}

Figure \ref{Fnowt} shows the time variations of rmsd from several
reference structure for TJF. As wee see, the {\em in vacuo} dynamics are
qualitatively different from that with explicit hydration. Namely, the
short-time rmsd fluctuations are similar or even larger,
but the trajectory apparently samples from a broad energy
valley around one and the same state in which it appeared
after equilibration. This system has an
essentially constant list of atom-atom contacts.
The local energy minima can occur due to complex energy profiles of
base stacking, but they seem to be rather shallow and the trajectory
passes above the barriers so that no distinct conformational
transitions within the B-DNA family are observed.
A few sharp peaks in the traces in Fig. \ref{Fnowt} correspond
to short-living $\rm B_I\to B_{II}$ transitions. The rmsd between
the structures averaged over the five consecutive nanosecond intervals
was as low as 0.1~\AA. The final TJF state shown in Fig. \ref{Falstr}
has rmsd of only 0.5~\AA\ from the initial minimum energy
structure.

Thus, a naked B-DNA duplex presents a relatively simple object
with a single significant energy minimum which is rather broad.
This state resembles the canonical B-DNA and is characterized
by the rmsd's of 3.1 and 4.5~\AA, respectively, from B and A forms
and some of its helical parameters shown in Table \ref{Thlpa}
are also between the canonical A and B-DNA values.
The other two {\em in vacuo} trajectories were very similar.
This means, in particular, that the shorter durations of TJE
and TJG are quite sufficient for correct averaging.
A visual analysis of the final TJF state superposed with the canonical
B-DNA shows that the relatively high rmsd
value results from the following main contributions: (i) lower
negative Xdisp, (ii) much lower average twist, and (iii) small but
non-negligible bending of the helical axis towards the center of the
minor groove. The latter effect is very clear in superposition figures
and it is nicely detected by the Curves algorithm. \cite{Curves:} The
first two features give this structure some A-like character although
one notes that in fact only these two helical parameters clearly
deviate towards the A form. There was only very minor difference
between the final states of these three {\em in vacuo} trajectories,
which is discussed hereafter.

\subsection{Negative Role of Water Films in Major Groove}

\begin{figure}
\centerline{\psfig{figure=FtjX.ps,width=7.5cm,angle=0.,%
bbllx=130bp,bblly=280bp,bburx=490bp,bbury=500bp,clip=t}}
\caption{\label{FtjX}
Comparison of the time variations of rmsd from reference
conformations for TJD and TJI. Ini and Fin mark the traces of
rmsd from the initial energy minimum and the final
MD state, respectively. The lower plates exhibit variation
of the ``proper'' and ``improper'' rmsd's from the crystal
structure \protect\cite{Drew:81} shown by the solid and broken
lines, respectively. All rmsd values are
given in angstr{\"o}ms.
}\end{figure}

The hydration procedure described above has been designed to
completely fill up the minor groove of B-DNA and also left in the
major groove some small number of water molecules to which we
initially did not pay much attention. The effect of
this water is not negligible, however. Figure \ref{FtjX} compares the
time variations of rmsd from reference conformations for TJD and TJI.
In both cases the initial state of was prepared by hydrating the
crystallographic structure which, because of the water crowding, has
only slightly changed during energy minimization (the rmsd 1.25~\AA). In the
crystal conformation, the two opposite phosphate traces come close to
each other in the middle, which is usually described as narrowing of
the minor groove. \cite{Kopka:83} For this specific shape our
hydration procedure left in the major groove a relatively large
number of water molecules.
The left two plates in Fig. \ref{FtjX} show that, in dynamics, the
structure slowly returns to the initial state after first leaving it
rapidly during heating. The lower plate shows the time variation of
the rmsd from the crystal structure. Due to the symmetry of the
sequence the rmsd comparison can be made in two ways giving the
``proper'' and ``improper'' rmsd values, the latter corresponding to
the opposite assignment of the two strands. The crystal conformation
itself is strongly asymmetric and characterized by a large
``improper'' rmsd of 1.7~\AA. In dynamics, this initial asymmetry
should relax, and the same level must be reached by the two rmsd
values after an ultimate equilibration. Figure \ref{FtjX} shows that
although the TJD trace of the ``proper'' rmsd reaches very low values
it is systematically shifted downward from the ``improper'' trace. In
other words, during these nine nanoseconds this symmetric DNA duplex
always ``remembers'' the asymmetry of its initial state.

The right two plates in Fig. \ref{FtjX} exhibits analogous traces for
TJI where all major groove water molecules were removed at the
beginning. The total number of water molecules appears larger
due to small variation of the hydration procedure, notably,
the initial 5~\AA\ shell was constructed with all crystallographic
waters considered as part of the duplex. Note that the structure again
moves out from the initial energy minimum and back, but the maximal
rmsd value reaches 2.8~\AA\, which is much larger than in the previous
case. The deviation from the crystal conformation is also one angstrom
larger, and, importantly, the two rmsd traces level after 3 ns
indicating that the ``memory'' of the initial asymmetric conformation
is essentially lost.

This example illustrates the general effect observed in many other
simulations. When the number of water molecules is not sufficient to
form a competent shell they cannot diffuse and tend to occupy fixed
positions thus stabilizing some random conformations and slowing down
the sampling. Very often they form films that connect several bases
and phosphate groups and can induce local bending of the helical axis.
To reduce these effects, in the last few trajectories in Table
\ref{Ttjli} the major grove was cleaned from the remaining water
molecules before heating. This does not mean that we consider
interactions in the major groove unimportant. The present DNA model
essentially performs the structural analysis of water surrounded by
the mobile DNA walls in the minor groove. The relatively small number
of these molecules makes possible nearly exhaustive sampling. In doing
so we implicitly assume that, due to much larger number of solvent
molecules involved in hydration of the major groove, its water shell
can easily take any necessary shape corresponding to a given DNA conformation.
Our calculations show that an empty major groove is a much better
approximation of such behavior than a faintly hydrated one.

\subsection{Characteristic Hydration Patterns}

Figure \ref{Fstereo} shows the final states of TJA and TJI with the
most stable hydration sites in the minor groove. Atom coordinates from
the trajectory points saved during the last nanosecond were
superimposed, and the position fluctuations were computed for all
water molecules. Figure \ref{Fstereo} shows average oxygen positions
for 20 least mobile water molecules in the minor groove, with the
sites close enough for hydrogen bonding joined by the thick lines. The
amplitudes of the position fluctuations range from 0.5 to 1.0~\AA\ and
are close to that of the neighboring DNA atoms indicating a tight
binding.  In case of TJA similar sites are also found in the major
groove, but as we noted above they are not representative and we do
not discuss them here. In TJA the slowest water molecules form a
continuous cluster at the bottom of the minor groove. In contrast, in
TJI several such clusters are separated by areas where water is more
mobile. In both cases, however, in the middle of the minor groove a
linear trace is found corresponding to the characteristic hydration
pattern often observed in crystals and referred to as the ``hydration
spine''.  \cite{Dickerson:82}

\begin{figure}
\centerline{\psfig{figure=Fstereo1.ps,width=7.5cm,angle=0.,%
bbllx=150bp,bblly=200bp,bburx=430bp,bbury=450bp,clip=t}}
\centerline{\psfig{figure=Fstereo2.ps,width=7.5cm,angle=0.,%
bbllx=150bp,bblly=280bp,bburx=430bp,bbury=530bp,clip=t}}
\caption{\label{Fstereo}
The final states of TJA and TJI with the most stable hydration
sites in the minor groove. The average oxygen positions
of 20 least mobile water molecules in the minor groove
are joined by the thick lines if they are close enough for
hydrogen bonding.
}\end{figure}

\begin{figure}
\centerline{\psfig{figure=Fspine1.eps,width=7.5cm,angle=0.,%
bbllx=10bp,bblly=10bp,bburx=375bp,bbury=175bp,clip=t}}
\centerline{\psfig{figure=Fspine2.eps,width=7.5cm,angle=0.,%
bbllx=10bp,bblly=10bp,bburx=375bp,bbury=175bp,clip=t}}
\caption{\label{Fspine}
Schematic representation of the time evolution of the
minor groove hydration patterns in TJA and TJI. The
hydrogen bonding contacts shown for the central eight base
pairs were taken from individual snapshots as explained in
the text. The filled circles denote water molecules that keep their
orientations in all subsequent snapshots.
}\end{figure}

Figure \ref{Fspine} shows schematically the time evolution of the
minor groove hydration patterns in these two trajectories. The
hydrogen bonding contacts shown for the central eight base
pairs were taken from individual snapshots as follows.
First, water molecules directly
bound to bases were considered and the shortest water bridges between
them were found. Only one and two member bridges are shown in the
figure. The filled circles denote water molecules that keep their
orientation in all subsequent snapshots. The orientation is
considered to be maintained if there is at least one chain to a
correct water-base contact. For clarity the opposite strands are
shifted by one step with respect to the Watson-Crick pairing.

In TJA most hydration sites were vacant after heating while already at
1.25 ns they are all occupied, with the canonical spine structure
formed in the middle which spans over six consecutive base pairs. The
total number of water molecules selected by the foregoing rules slowly
increases during the whole trajectory; the spine shortens and is
gradually replaced by less economical patterns near both ends. The
number of sites reaches 21 after 5 ns, but we note that the TJA
plate in Fig.  \ref{Fspine} actually shows 32 different water
molecules, which manifests a certain level of exchange with the
environment. It is readily seen that the three water molecules that
form the first layer of the hydration spine are the most stable in
this system. At the same time one notes that a transition from spine
to ribbon type of hydration \cite{Edwards:92} occurs at the upper
border of the spine.

The initial state for the TJI was prepared from the Xray structure
with new water molecules added to already occupied experimental
hydration cites. That is why in TJI the minor groove hydration is
rather dense already after heating. The total number of solvent
molecules selected by the chosen rules fluctuates at the level reached
in TJA after 5 ns. As a whole the evolution of hydration patterns is
rather different in TJI compared to TJA, but still water molecules
at the bottom of the hydration spine are among the least movable. Note
also that in this case a reversible transition from the spine to the
ribbon type of hydration occurs at the upper spine border.

Spontaneous formation of the hydration spine
in the middle of the minor groove has
been observed in all trajectories computed for the {\em Eco}RI dodecamer.
This characteristic water arrangement has been first encountered in
the crystal structure of this duplex \cite{Dickerson:82}
and is generally considered as a characteristic
feature of certain AT-rich sequences that tend to adopt an unusual
conformation with a particularly narrow minor groove. \cite{Dickerson:92} We
may conclude, therefore, that the present B-DNA model nicely
reproduces this well established experimental observation, which
is the most probable origin of a surprisingly good correspondence between the
computed and experimental DNA conformations. It can be
noted also that the water structure in the minor groove is rather
persistent and changes very slowly. As we discussed above, due to the
initial ``cleaning'' of the major groove in TJI the DNA molecule
lost the ``memory'' of the initial state much faster than in TJD.
Figure \ref{Fspine} demonstrates, however, that throughout the
trajectory the hydration patter is more spine-like in the lower part
of the minor groove and more ribbon-like in the upper part, the
feature inherited from the initial Xray structure where the spine
expands from the center two times farther towards the narrower end of
the minor groove. \cite{Drew:81a}

It should be noted, finally, that the spine may be also a model
dependent feature. We have noticed that, although trajectories always
start from a well minimized structure, the dynamics usually begins
with rapid unwinding and extending of the DNA molecule, which narrows
the minor groove and pushes much of the water out. This occurs
regardless of the sequence and the spine is normally the first
hydration pattern that sets up afterwards. The subsequent slow
diffusion of water molecules into the minor groove usually takes
several nanoseconds and often goes faster near the ends than in the
middle. Nevertheless, the high stability of this specific hydration
pattern in the center of the {\em Eco}RI dodecamer as well as for some other
sequences like A-tracts, which are not considered here, is a sequence
dependent feature. To support the last assertion we included in Tables
\ref{Ttjli} and \ref{Thlpa} and also in Figs. \ref{Falstr} and
\ref{Falmg} trajectory TJL computed in similar conditions for a
GC-alternating dodecamer. Table \ref{Thlpa} shows that, except for the
twist, the final TJL state has average helical parameters even closer
to the canonical B-DNA than {\em Eco}RI dodecamer. Notably, its inclination
is negative and Xdisp is much less than in any {\em Eco}RI trajectory. In
terms of rmsd it is similarly close to the canonical B-DNA. At the same
time, the hydration patterns in the minor groove were very different
from those discussed above and it is seen in Figs. \ref{Falstr} and
\ref{Falmg} that the final TJL state has the minor groove evenly wider
than in the canonical B-DNA form. The minor groove was narrow in the
initial phase of the trajectory and then gradually opened. The
distinct properties of regular sequences in this DNA model will be
discussed in detail elsewhere.

\subsection{Long Range Electrostatics}
The most surprising property of the present DNA model is that its
trajectories converge to conformations that much better compare with
experimental data than in more expensive calculations with explicit
counterions and rigorous treatment of electrostatic interactions.
\cite{Mzjacs:98} This would be less surprising if the adjustable
parameters of the simple {\em ad hoc} treatment employed here were
specifically fitted to experimental data, but, in fact, we have
arbitrarily chosen one of the options earlier used in conformational
analysis of DNA, not even the best recommended. \cite{Flatters:97}
These observations suggest that the precision of calculation of long
range electrostatic interactions may be generally less important for
DNA structure than it is sometimes supposed. A comprehensive analysis
of this fundamental issue is beyond the capabilities of this simple
approach, nevertheless, it is interesting to check how this B-DNA
model responds to variation of parameters involved in calculation of
electrostatic interactions.

The {\em in vacuo} trajectories TJE, TJF, and TJG were computed with
phosphate charges varied from -0.25 to -0.75. They produce very
similar structures with rmsd between the final TJE and TJG states
around 0.5~\AA\ and the final TJF state at roughly 0.3~\AA\ from the
both. Table \ref{Thlpa} exhibits very small, but regular variation of
helical parameters. Namely, growing phosphate repulsion moves Xdisp,
inclination and twist closer to the canonical B-DNA values, but
propeller and rise in the opposite direction. The rmsd's of the
final MD states from the reference experimental DNA conformations are
practically identical for all three trajectories.

The same range of phosphate charges is covered by trajectories TJH,
TJI, TJJ and TJK where the minor groove has been hydrated. Table
\ref{Thlpa} shows that in this case variations of helical parameters
are stronger and less regular, but some of them, notably, Xdisp,
inclination and twist exhibit the same trends as {\em in vacuo}.
Reduction of phosphate repulsion accentuates already existing small
deviations from canonical B-DNA values towards A-DNA. The increased
repulsion has an opposite effect and one may expect that, in terms of
helical parameters, the best fitting of the computed and experimental
conformations can be achieved with phosphate charges slightly below
-0.6 aeu. Unexpectedly, however, variation of phosphate charges in an
explicitly hydrated B-DNA model significantly affects the overall
bending of the molecule which is seen in Fig. \ref{Falstr}. In fact,
all computed conformations of the {\em Eco}RI dodecamer except TJB are
slightly bent towards the center of the minor groove. It appears that
the increase in the phosphate repulsion also increases this bending,
which is a non-trivial effect of the minor groove hydration since no
such trend is observed in the {\em in vacuo} conditions.
This effect deserves careful investigation because it in fact corresponds
well to certain experimental observations. Gel retardation experiments
with sequence repeats indicate that, in solution, the {\em Eco}RI
dodecamer should be bent towards the center of the minor groove.
\cite{Diekmann:88} It is also known that DNA bending is generally reduced in
high salt \cite{Diekmann:87} which also reduces the persistence length of
DNA, \cite{Hagerman:88} with the latter effect attributed to reduced
phosphate repulsion due to better screening. Comparison of the final
states of TJH, TJI, TJJ and TJK suggests that, in the minimal B-DNA,
bending is also somehow connected with phosphate repulsion
and increased backbone stiffness also favors bending.

Because of the foregoing effect the rmsd's from experimental B-DNA
conformations does not follow the corresponding deviations
of helical parameters. In the case of TJJ, for instance,
the increased bending overshadows the small improvement in helical
parameters with respect to the canonical B-DNA.
These results generally show, however, that significant variations of the most
arbitrary parameter involved in our simplistic treatment of the
electrostatic effects still result in structures that are close
to experimental B-DNA both in terms of helical parameters and atom rmsd's.
Note that the largest rmsd in Fig.
\ref{Falstr} observed for TJH is lower than in all reported
simulations with the particle-mesh Ewald method. \cite{Mzjacs:98}
Although its helical parameters
noticeably deviate from canonical values it still qualitatively
reproduces the characteristic shape of the minor groove with a short
hydration spine in the center.

\subsection{Concluding Discussion}
The minimal model of B-DNA explored here consists of two qualitatively
different components, namely, the DNA duplex
and a water cloud docked in the minor groove.
Dynamics simulations of a naked B-DNA duplex show that it presents a
relatively simple object with a broad energy valley around a single
state with an intermediate conformation between A and B DNA forms.
This state is ``fragile'' in the sense that $\rm B_I\to
B_{II}$ transitions, which are generally believed to introduce only
subtle nuances in the B-DNA structure, \cite{Berman:97} cause gross
structural transitions to an essentially different DNA form with a
collapsed minor groove. \cite{Swaminathan:91} With a water cloud docked in
the minor groove the system acquires qualitatively different properties.
Its potential energy landscape becomes rough, with multiple
significant local minima corresponding to
rearrangement of atom-atom contacts in water and at
the DNA-water interface. Dynamics simulations of this system can be
viewed as structural analysis of water near a mobile surface that
carries a net negative charge and has a specific distribution of
hydrogen bonding contacts.

The main advantage of these simulation conditions is the
possibility of practically exhaustive sampling with the presently
available computer resources. It should be noted that the
sampling and the affordable duration of trajectories are not strict
synonyms. The number of essential degrees of freedom is another
important parameter that should not be ignored. We expect, for
instance, that, in otherwise similar conditions, the rate of
convergence would be lower for a B-DNA duplexes with both minor and
major grooves hydrated explicitly. That is why, in our opinion, the
minimal B-DNA should be distinguished and used, notably,
for equilibration of the surface water in preparative phases of
more expensive DNA simulations.

The properties of this system are non-trivial and sometimes
counter-intuitive. We have shown here that the canonical B-DNA form
presents a strong point of attractions in its conformational space, so
that numerous trajectories of the {\em Eco}RI dodecamer and one long
trajectory of the GC-alternating dodecamer converge to conformations
with nearly canonical helical parameters. We have shown that the
sequence specific features of the {\em Eco}RI dodecamer are nicely
reproduced, namely, the profile of the minor groove and spontaneous
formation of the characteristic ``spine'' hydration pattern
which is generally considered as an indispensable part of this
structure. \cite{Dickerson:82} To our knowledge, none of the numerous earlier
theoretical studies of this molecule could reach a comparable
level of agreement with experimental data.
All these observations justify the minimal model of B-DNA, confirm its
credibility and suggest that it has a potentially wide scope of application.


\section{Introduction}

Language identification is an example of a general class of problems in which
we want to assign an input data stream to one of several categories as quickly
and accurately as possible. It can be solved using many techniques, including
knowledge-poor statistical approaches. Typically, the distribution of n-grams
of characters or other objects is used to form a model. A comparison of the
input against the model determines the language which matches best. Versions
of this simple technique can be found in Dunning~(1994) and Cavnar and
Trenkle~(1994), while an interesting practical implementation is described by
Adams and Resnik~(1997).

A variant of the problem is considered by Sibun and Spitz~(1994), and Sibun
and Reynar~(1996), who look at it from the point of view of Optical Character
Recognition (OCR). Here, the language model for the OCR system cannot be
selected until the language has been identified. They therefore work with
so-called shape tokens, which give a very approximate encoding of the
characters' shapes on the printed page without needing full-scale OCR. For
example, all upper case letters are treated as being one character shape, all
characters with a descender are another, and so on. Sequences of character
shape codes separated by white space are assembled into word shape
tokens. Sibun and Spitz then determine the language on the basis of linear
discriminant analysis (LDA) over word shape tokens, while Sibun and Reynar
explore the use of entropy relative to training data for character shape
unigrams, bigrams and trigrams. Both techniques are capable of over 90\%
accuracy for most languages. However, the LDA-based technique tends to perform
significantly worse for languages which are similar to one another, such as
the Norse languages. Relative entropy performs better, but still has some
noticeable error clusters, such as confusion between Croatian, Serbian and
Slovenian.

What these techniques lack is a measure of when enough information has been
accumulated to distinguish one language from another reliably: they examine
all of the input data and then make the decision. Here we will look at a
different approach which attempts to overcome this by maintaining a measure of
the total evidence accumulated for each language and how much confidence there
is in the measure. To outline the approach:
\begin{enumerate}
\item The input is processed one (word shape) token at a time. For each
language, we determine the probability that the token is in that language,
expressed as a 95\% confidence range.
\item The values for each word are accumulated into an overall score with a
confidence range for the input to date, and compared both to an absolute
threshold, and with each other. Thus, to select a language, we require not
only that it has a high score (probability, roughly), but also that it is
significantly better scoring than any other.
\item If the process fails to make a decision on the data that is available,
the subset of the languages which have exceeded the absolute threshold can be
output, so that even if a final decision has not been made, the likely
possibilities have been narrowed down.
\end{enumerate}
We look at this procedure in more detail below, with particular emphasis on
how the underlying statistical model provides confidence
intervals. An evaluation of the technique on data similar to that used by
Sibun and Reynar follows\footnote{Some ideas related to the use of confidence
limits can also be found in Dagan {\it et al.}~(1991), applied in a different
area.}.

\section{The Identification Algorithm}

The essential idea behind the identification algorithm is to accumulate the
probability of the language given the input tokens for each language, treating
each token as an independent event. To obtain the probability of a language
$l$ given a token $t$, $p(l|t)$, we use Bayes' rule:
\[ p(l|t) = \frac{p(t|l) p(l)}{p(t)} \]
where $p(t|l)$ is the probability of the token if the language is known,
$p(t)$ is the {\em a priori} probability of the token, and $p(l)$ is the {\em
a priori} probability of the language. We will assume that $p(l)$ is constant
(all languages are equi-probable) and drop it from the computation; in the
tests, we will use the same amount of training data for each language. The
other two terms are estimated from training data, using the procedure
described in section~\ref{train}.

\subsection{The language model and the algorithm}

The input to the algorithm consists of a stream of tokens, such as word shape
tokens (as in Sibun and Spitz, or Sibun and Reynar) or words themselves. The
model for each language contains the probability of each known token given the
language, expressed as three values: the basic probability, and the lower and
upper limits of a range containing this probability for a specific level of
confidence. We will denote these by $p_{B}(t|l)$, $p_{L}(t|l)$, $p_{H}(t|l)$,
for base, low and high values. The probability that a token which has never
been seen before is in a language is also present in the model of the
language. In addition, there is a language independent model, containing the
$p(t)$ values. No confidence range is used for them, although this would be a
simple extension of the technique.

The algorithm proceeds by processing tokens, building up evidence about each
language in three accumulators. The accumulators represent the overall
probability of the language given the entire stream of tokens to date, again
as base, low and high values, denoted $a_{B}(l)$, $a_{L}(l)$, $a_{H}(l)$. They
are set to zero at the start of processing, and the logarithms of the
probabilities are added to them as each token is processed. By taking
logarithms of probabilities, we are in effect measuring the amount of evidence
for each language, expressed as information content. From a practical point of
view, using logarithms also helps keep all the values in a reasonable range
and so avoids numerical underflow.

After processing each token, two tests are applied. Firstly, we examine the
base accumulator for the language which has the highest accumulated total, and
test whether it is greater than a fixed threshold, called the activation
threshold. If it is, then we conclude that enough information has been
accumulated to try to make a decision. The low value for this language
$a_{L}(l)$ is then compared against the high value $a_{H}(l')$ for the next
best language $l'$, and if $a_{L}(l)$ exceeds $a_{H}(l')$ language $l$ is
output and the algorithm halts. Otherwise, the process continues with the next
token, until the best choice language is a clear ``winner'' over any other.

Finally, if we reach the end of the input data without a decision being made,
several options are possible, depending on the needs of the application. We can
simply output the language with the highest base score, even if the second
test is not satisfied. Alternatively, we can output the highest scoring
language, and all other languages whose high probability is greater than the
low probability of this language.

\subsection{Training the model}
\label{train}

The model is trained using a collection of corpora for which the correct
language is known. For a given language $l$ and token $t$, let $f(t,l)$ be the
count of the token in that language and $f(l)$ be the total count of all tokens
in that language. $f(t)$ is the count of the token $t$ across all the
languages, and $F$ the count of all tokens across all languages. The
probability of the token occurring in the language $p(t|l)$ is then calculated
by assuming that the probabilities follow a binomial distribution. The idea
here is that token occurrences are binary ``events'' which are either the
given token $t$ or are not. For large $f(t,l)$, the underlying probability
can be calculated by using the normal approximation to the binomial, giving
the base probability
\[p_{B}(t|l) = \frac{f(t,l)}{f(l)}\]
The standard deviation of this quantity is
\[\sigma(t,l) = \sqrt{f(l)p_{B}(t|l)(1-p_{B}(t|l))}\]
The low and high probabilities are found by taking a given number of standard
deviations $d$ from the base probability.
\[p_{L}(t|l) = \frac{f(t,l)-d\sigma(t,l)}{f(l)}\]
\[p_{H}(t|l) = \frac{f(t,l)+d\sigma(t,l)}{f(l)}\]
In the evaluation below, $d$ was set to 2, giving 95\% confidence limits.

For lower values of $f(t,l)$, the calculation of the low and high
probabilities can be made more exact, by substituting them for the base
probability in the calculation of the standard deviation, giving
\[p_{L}(t|l) = \frac{f(t,l)-d\sqrt{f(l)p_{L}(t|l)(1-p_{L}(t|l))}}{f(l)}\]
\[p_{H}(t|l) = \frac{f(t,l)+d\sqrt{f(l)p_{H}(t|l)(1-p_{H}(t|l))}}{f(l)}\]
Approximating $1-p_{L}(t|l)$ and $1-p_{H}(t|l)$ to 1 on the grounds that the
probabilities are small, and solving the equations gives
\[p_{L}(t|l) = \frac{(\sqrt{d^{2}+4f(t,l)}-d)^{2}}{4f(l)} \]
\[p_{H}(t|l) = \frac{(\sqrt{d^{2}+4f(t,l)}+d)^{2}}{4f(l)} \]
The calculation requires marginally more computational effort than the first
case, and in practice we use it for all but very large values of $f(t,l)$,
where the approximation of $1-p_{L}(t|l)$ and $1-p_{H}(t|l)$ to 1 would break
down.

For very small values of $f(t,l)$, say less than 10, the normal
approximation is not good enough, and we calculate the probabilities by
reference to the binomial equation for the probability of $m$ (=$f(t,l)$)
successes in $n$ (= $f(l)$) trials:
\[p(m) = \frac{p^{m}(1-p)^{n-m}n!}{m!(n-m)!} \]
$p$ is the underlying probability of the distribution, and this is what we are
after. By choosing values for $p(m)$ and solving to find $p$ we can
obtain a given confidence range. To obtain a 95\% interval, $p(m)$ is set to
0.025, 0.5 and 0.975, yielding $p_{L}(t|l)$, $p_{B}(t|l)$, and $p_{H}(t|l)$,
respectively. In fact, this is not exactly how the probability ranges for low
frequency items should be calculated: instead the cumulative probability
density function should be calculated and the range estimated from
it\footnote{Thanks to one of the referees for pointing this out.}. For the
present purposes, the low frequency items do not make much of a contribution
to the overall success rate, and so the approximation is unimportant. However,
if similar techniques were applied to problems with sparser data, then the
procedure here would have to be revised.

Finally, we need a probability for tokens which were not seen in the training
data, called the zero probability, for which we set $m=0$ in the above
equation giving
\[p(0|l) = 1-\sqrt[n]{p(m)}\]
It is not clear what it means to have a confidence measure here, and so we use
a single value for base, low and high probabilities, obtained by setting $p(m)$
to 0.95.

Similar calculations using $f(t)$ in place of $f(t,l)$ and $F$ in place of
$f(l)$ give the a priori token probabilities $p(t)$. As already noted, base,
low and high value could have been calculated in this case, but as a minor
simplification, we use only the base probability.

\section{Evaluation}

To evaluate the technique, a test was run using similar data to
Sibun and Reynar. Corpora for eighteen languages from the European Corpus
Initiative CDROM 1 were extracted and split into
non-overlapping files, one containing 2000 tokens\footnote{Sibun and Spitz,
and Sibun and Reynar, present their results in terms of lines of input, with
1-5 lines corresponding roughly to a sentence, and 10-20 lines to a
paragraph. Estimating a line as 10 words, we are therefore working with
significantly smaller data sets.}, one containing 200 tokens, and 25 files
each of 1, 5, 10 and 20 tokens. The 2000 and 200 token files
were used as training data, and the remainder for test data. Wherever possible
the texts were taken from newspaper corpora, and failing that from novels or
literature. The identification algorithm was run on each test file and the
results placed in one of four categories:
\begin{itemize}
\item Definitive, correct decision made.
\item No decision made by the end of the input, but highest scoring language
was correct.
\item No decision, highest scoring language incorrect.
\item Definitive, incorrect decision made.
\end{itemize}
The sum of the first two figures divided by the total number of tests gives a
measure of accuracy; the sum of the first and last divided by the total gives
a measure of decisiveness, expressed as the proportion of the time a
definitive decision was made. The tests were executed using word shape tokens
on the same coding scheme as Sibun and Reynar, and using the words as they
appeared in the corpus. No adjustments were made for punctuation, case,
etc. Various activation thresholds were tried: raising the threshold increases
accuracy by requiring more information before a decision is made, but reduces
decisiveness. With shapes and 2000 tokens of training data, at a
threshold of 14 or more, all the 20 token files gave 100\% accuracy. For words
themselves, the threshold was set to 22. The results of these tests appear in
table~\ref{tab1}. The figures for the activation threshold were determined by
experimenting on the data. An interesting area for further work would be to
put this aspect of the procedure on a sounder theoretical basis, perhaps by
using the {\em a priori} probabilities of the individual languages.

\begin{table*}[htb]
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|}\hline
Test and & \multicolumn{5}{c|}{Accuracy (\%)} &
\multicolumn{5}{c|}{Decisiveness (\%)} \\ \cline{2-6} \cline{7-11}
threshold  &  \multicolumn{5}{c|}{Tokens of test data} &
\multicolumn{5}{c|}{Tokens of test data} \\ \cline{2-6} \cline{7-11}
 & 1 & 5 & 10 & 20 & {\bf All} & 1 & 5 & 10 & 20 & {\bf All} \\ \hline
Tokens (0) & 71.6 & 72.7 & 69.6 & 72.0 & {\bf 71.4} & 88.0 & 99.3 & 100 & 99.8 & {\bf 96.8}\\
Tokens (10) & 92.9 & 98.4 & 98.4 & 98.2 & {\bf 97.0} & 66.0 & 98.9 & 99.6 & 99.8 & {\bf 91.1}\\
Tokens (14) & 94.2 & 99.6 & 99.1 & 100 & {\bf 98.2} & 49.8 & 98.9 & 99.6 & 99.8 & {\bf 87.0}\\
\hline
Words (0) & 78.4 & 80.4 & 77.1 & 78.7 & {\bf 78.7} & 97.3 & 100 & 100 & 100 & {\bf 99.3}\\
Words (10) & 95.8 & 97.6 & 97.1 & 98.0 & {\bf 97.1} & 76.9 & 99.8 & 100 & 99.8 & {\bf 94.1}\\
Words (22) & 96.9 & 99.8 & 99.8 & 100 & {\bf 99.1} & 29.3 & 98.9 & 99.8 & 99.8 & {\bf 81.9}\\
\hline
\end{tabular}
\caption{Performance with 2000 tokens of training data}
\label{tab1}
\end{center}
\end{table*}
The accuracy figures are generally similar to or better than those of Sibun
and Reynar. The corresponding figures for 200 tokens of training data appear in
table~\ref{tab1.200}, for the token identification task only.
\begin{table*}[htb]
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|}\hline
Threshold &  \multicolumn{5}{c|}{Accuracy (\%)} &
\multicolumn{5}{c|}{Decisiveness (\%)} \\ \cline{2-6} \cline{7-11}
          &  \multicolumn{5}{c|}{Tokens of test data} &
\multicolumn{5}{c|}{Tokens of test data} \\ \cline{2-6} \cline{7-11}
          & 1 & 5 & 10 & 20 & {\bf All} & 1 & 5 & 10 & 20 & {\bf All} \\ \hline
0 & 63.3 & 72.4 & 48.9 & 47.1 & {\bf 57.9} & 72.2 & 89.6 & 96.7 & 96.2 & {\bf 88.7}\\
5 & 82.2 & 88.9 & 75.6 & 75.8 & {\bf 80.6} & 58.0 & 86.4 & 93.3 & 93.8 & {\bf 82.9}\\
10 & 86.0 & 94.0 & 88.0 & 87.6 & {\bf 88.9} & 45.6 & 85.1 & 91.1 & 92.4 & {\bf 78.6}\\
\hline
\end{tabular}
\caption{Performance with 200 tokens of training data (word shape tokens only)}
\label{tab1.200}
\end{center}
\end{table*}

One of the strengths of the algorithm is that it makes a decision as soon as
one can be made reliably. Table~\ref{tab2} shows the average number of tokens
which have to be read before a decision can be made, for the cases where the
decision was correct and incorrect, and for both cases together. Again, the
results are for word shape tokens, and for words alone.
\begin{table*}[htb]
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|}\hline
Threshold & \multicolumn{3}{c|}{Shape tokens} & \multicolumn{3}{c|}{Words} \\
\cline{2-4} \cline{5-7}
   & Correct & Incorrect & All & Correct & Incorrect & All \\ \hline
0  & 3.22 & 1.23 & 2.65 & 1.81 & 1.07 & 1.66 \\
10 & 7.33 & 4.55 & 7.28 & 5.31 & 3.88 & 5.28 \\
14 & 9.35 & 6.50 & 9.33 &      &      &      \\
22 &      &      &      & 10.6 & 8.00 & 10.6 \\
\hline
\end{tabular}
\caption{Average number of tokens read before convergence}
\label{tab2}
\end{center}
\end{table*}
The figures show that convergence usually happens within about 10 words, with
a long tailing off to the results. The longest time to convergence was 153
shape tokens.

A manual inspection of one run (2000 lines of training data, tokens,
threshold=14) shows that errors are somtimes clustered, although quite
weakly. For example, Serbian, Croatian and Slovenian show several confusions
between them, as in Sibun and Reynar's results. There are two observations to
be made here. Firstly, there are about as many other errors between these
language and languages which are unrelated to them, such as Italian, German
and Norwegian, and so the errors may be due to poor quality data rather than a
lack of discrimination in the algorithm. For example, Croatian is incorrectly
recognised as Serbian 3 times and as Slovenian once, while the languages which
are misrecognised as Croatian are German and Norwegian (once each). Secondly,
even where there are errors, the range of possibilities has been substantially
reduced, so that a more powerful process (such as full-scale OCR followed by
identification on words rather than shape tokens, or a raising of the
threshold and adding more data) could be brought in to finish the job
off. That is, the confidence limits have provided a benefit in reducing the
search space. The confusion matrix for this case appears in an appendix.

\subsection{Broader applicability}

Although the algorithm was developed with language identification in mind, it
is interesting to explore other classification problems with it. A simple and
rather crude experiment in ``genre'' identification was carried out, using the
Brown corpus. Each section of the corpus (labelled A, B, C ... R in the
original) was taken as a genre, and files of similar distribution to the
previous experiment were extracted. Because this is a more unconstrained
problem, the training set and tests sets were about 10 times the size of the
language identification task. A 20000 word file was used as training data, and
the remaining files as test data. Accuracy and decisiveness results appear in
table~\ref{tab3}.
\begin{table*}[htb]
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|}\hline
Threshold &  \multicolumn{5}{c|}{Accuracy (\%)} &
\multicolumn{5}{c|}{Decisiveness (\%)} \\ \cline{2-6}\cline{7-11}
          &  \multicolumn{5}{c|}{Words of test data} &  \multicolumn{5}{c|}{Words of test data} \\\cline{2-6}\cline{7-11}
          & 10 & 50 & 100 & 200 & {\bf All} & 10 & 50 & 100 & 200 & {\bf All} \\ \hline
0 & 47.7 & 76.0 & 83.7 & 80.8 & {\bf 72.1} & 36.3 & 38.7 & 39.5 & 42.9 & {\bf 39.3}\\
10 & 50.9 & 86.9 & 96.8 & 99.5 & {\bf 83.5} & 2.13 & 15.5 & 16.5 & 18.1 & {\bf 13.1}\\
12 & 50.9 & 86.9 & 96.8 & 99.7 & {\bf 83.6} & 1.07 & 14.1 & 14.9 & 16.0 & {\bf 11.5}\\
\hline
\end{tabular}
\caption{Performance on genre identification}
\label{tab3}
\end{center}
\end{table*}
Beyond the activation threshold of 12, there is no significant improvement in
accuracy. The technique seems to give good accuracy when there is sufficient
input (100 words or more), but at the cost of very low decisiveness. Excluding
a fixed list of common words such as function words might increase the
decisiveness. These results should be taken with a pinch of salt, as the
notion of genre is not very well-defined, and it is not clear that sections of
the Brown corpus really represent coherent categories, but they may provide a
starting point for further investigation.

\subsection{On decisiveness}

Decisiveness represents the degree to which a unique decision has been made
with a high degree of confidence. In cases where no unique decision has been
made, the range of possibilities will often have been reduced: a category is
only still possible at any stage if its high accumulator value is greater than
the low accumulator value of the best rated category. To illustrate this, the
number of categories which are still possible when all the input was exhausted
was examined. The results appear in tables~\ref{tab4} and ~\ref{tab5}, for the
tests of language identification from word shape tokens with an activation
threshold of 14 and a training set of 2000 tokens, and for genre
identification with a threshold of 12 and a training set of 20000
tokens. Results are shown for the cases of a correct decision, an incorrect
one, and all cases.
\begin{table}[htb]
\begin{center}
\begin{tabular}{|c|r|r|r|}\hline
Languages & \multicolumn{3}{c|}{Number of tests}\\ \cline{2-4}
remaining  & Correct & Incorrect & All \\ \hline
1 & 1560 & 6 & 1566 \\
2 & 128 & 7 & 135 \\
3 & 37 & 9 & 46 \\
4 & 18 & 2 & 20 \\
5 & 5 & 1 & 6 \\
6 & 5 & 1 & 6 \\
7 & 2 & 0 & 2 \\
8 & 2 & 0 & 2 \\
9 & 1 & 0 & 1 \\
10 & 3 & 0 & 3 \\
11 & 1 & 0 & 1 \\
12 & 2 & 0 & 2 \\
13 & 2 & 0 & 2 \\
17 & 1 & 0 & 1 \\
18 & 7 & 0 & 7 \\
\hline
\end{tabular}
\caption{Categories remaining at end of input (language identification from
word shape tokens)}
\label{tab4}
\end{center}
\end{table}
\begin{table}[htb]
\begin{center}
\begin{tabular}{|c|r|r|r|}\hline
Genres & \multicolumn{3}{c|}{Number of tests}\\ \cline{2-4}
remaining  & Correct & Incorrect & All \\ \hline
1 & 173 & 0 & 173 \\
2 & 22 & 0 & 22 \\
3 & 31 & 1 & 32 \\
4 & 45 & 3 & 48 \\
5 & 34 & 4 & 38 \\
6 & 65 & 8 & 73 \\
7 & 73 & 4 & 77 \\
8 & 83 & 3 & 86 \\
9 & 84 & 3 & 87 \\
10 & 89 & 2 & 91 \\
11 & 84 & 0 & 84 \\
12 & 128 & 1 & 129 \\
13 & 131 & 0 & 131 \\
14 & 175 & 1 & 176 \\
15 & 253 & 0 & 253 \\
\hline
\end{tabular}
\caption{Categories remaining at end of input (genre identification)}
\label{tab5}
\end{center}
\end{table}
The average number of possibilities remaining is 1.3 out of 18 for the
language identification test, and 9.7 out of 15 for the genre test, showing
that we are generally near to convergence in the former case, but have only
achieved a small reduction in the possibilities in the latter, in keeping with
the generally low decisiveness.

\subsection{A further comparison}

The classification algorithm described above was originally developed in
response to Sibun and Spitz's work. There is another approach to language
identification, which has a certain amount in common with ours, described in a
patent by Martino and Paulsen~(1996). Their approach is to build tables of the
most frequent words in each language, and assign them a normalised score,
based on the frequency of occurrence of the word in one language compared to
the total across all the languages. Only the most frequent words for each
language are used. The algorithm works by accumulating scores, until a preset
number of words has been read or a minimum score has been reached. They also
apply the technique to genre identification. Since there is a clear
similarity, it is perhaps worth highlighting the differences. In terms of the
algorithm, the most important difference is that no confidence measures are
included. The complexities of splitting the data into different frequency
bands for calculating probabilities are thus avoided, but no
test analogous to overlapping confidence intervals can be applied. Martino and
Paulsen say they obtain a high degree of confidence in the decision after
about 100 words, without saying what the actual success rate is; we can
compare this with around 10 words (or tokens) for convergence here.

\section{Conclusions}

We have examined a simple technique for classifying a stream of input tokens
in which confidence measures are used to determine when a correct decision can
be made. The results in table~\ref{tab1} show that there is a tradeoff between
accuracy and the degree to which the algorithm selects a single language. Not
surprisingly, the amount of training data also affects the performance, with
2000 tokens being adequate for accuracy close to 100\%, and convergence
typically being reached in the first 10 tokens. On a more unconstrained
problem, such as genre identification from words alone, the algorithm performs
less well in both accuracy and decisiveness even with significantly more
training data, and is probably not adequate except as a preprocessor to some
more knowledge intensive technique.

In a sense, language identification is not a very interesting problem. As we
have noted, there are plenty of techniques which work well, each with its own
characteristics and suitability for different application areas. What is
perhaps more important is the way the statistical information has been used
here. When we take a statistical or data-led approach to NLP, there are two
things which can help us trust that the technique is accurate. The first is a
belief that the statistical technique is an adequate model of the underlying
process which ``generates'' the data, using theoretical considerations or some
external source of knowledge to inform this belief. The second is quantitative
evaluation on test data which has been characterised by an outside source (for
example, in the case of part of speech tagging, a corpus which has been
manually annotated, or at least automatically tagged and manually
corrected). The problem with quantitative evaluation is that we do not know
whether it will generalise, so that if we train on one data set, we have only
the theoretical model to reassure that the same model will work on a different
data set. The idea I have been presenting here is to get the statistical
process itself to provide feedback about itself, through the use of confidence
limits which are themselves based in the statistical model. In doing so, we
hope to avoid presenting a result for which we lack adequate evidence.
\newline
\newline
{\bf Acknowledgement}
\newline
Thanks to Robert Keiller of Canon Research Centre Europe for his advice on
computing accurate statistics.


\section*{Acknowledgments}
K. C and K. W. H.  are  supported in part by KOSEF
grant 98-0201-004-2, KRF Distinguished Scholar Exchange Program
and KRF  BSRI Program 1998-015-D00071.
J. H. wishes to thank Prof. Kamionkowski for informing us his work
during 19th Texas Symposium in Paris.