%what I did this spring 2005 lateset draft
\documentclass[10pt]{article}

\usepackage{graphicx, verbatim, ulem, epsfig}

\oddsidemargin 0.0in
\evensidemargin 1.0in
\textwidth 6.5in
%\headheight 0.5in
%\topmargin 0.25in
\textheight 9.0in
%\footheight 1.0in

\begin{document}
\author{Fayette Shaw \\ \texttt{fayshaw@u.washington.edu}}
\date{\today}
\title{\LARGE Reduced-order Model for DNA Simulations \\
{\it Self-Organizing Systems Group} \\ Spring 2005 Status Report}
\maketitle

\tableofcontents
%\listoffigures
\begin{abstract}
As simulation in the DNA Machines project develops, research becomes motivated by two hypothesis:
(1) nucleotides can be modeled as rigid bodies and (2) this model will be comparable to current state 
of the art techniques.  As proposed in the winter status report, the first step of identifying
rigid body simplifications has taken place for Cytosine and Guanine bonded along a backbone.
Comparisons to other modeling techniques
must take place as well as proposing new experiments.
\end{abstract}


\section{Objective}
Current methods for computing energies and reaction pathways for large molecules are computational demanding.
State of the art computation techniques trade efficiency for accuracy, illusrating a need for
an algorithm with both capabilities.  We propose a new model for DNA folding by reprsenting nucleotides
as rigid bides, which may reduce computational demands without sacrificing accuracy,

Included in this report is the work done in AA549 Estimation and System Identification in taking the 
first steps toward identifying parameters in a reduced order model.  While this project did not have
satisfactory results, the procedure will be revisited.  Additionally, this method of identifying
spring constants brings up several questions.  How might this approach be extrapolated for larger
systems?  As this problem only involved five bodies, there were only ten spring constants to 
identify.  Increasing each ith additional body introduces i-1 parameters.  It may be possible to implement
a cutoff region beyond which forces are neglibly small.  Also, the covariance matrix needs to be calculated 
for various combinations of A,T,C,G, along a backbone as well as for the A-T and C-G hydrogen bond pairings.

\section{Representation}
The following representation is taken from Yakushevich\cite{Yakushevich}.  There are many different
ways to represent the various scales encompassed by molecular dynamics simulations with
respect to chemical reactions.

\begin{tabular}[h]{|c|c|c|c|c|c|}
\hline
time scale  & picosecond     & nanosecond   & microsecond   & milliseond     & second \\
\hline
interaction & oscillation    & bending      & bending       & dissociation   & amplitudes \\
\hline
energy      & 0.6 kcal / mol & 2/5 kcal/mol & 5/20 kcal/mol & 10/50 kcal/mol & -  \\
\hline
experimental & NMR, IRS & NMR, EPR & NMR, EPR & NMR, hydrogen exchange & bulk  \\
\hline
theoretical  & harmonic & harmonic  & helix-coil & topological & - \\
\hline
\end{tabular}




%\input{kinetics}

\section{Approach}

The bond elongation parameters can be calculated by applying an extended Kalman filter to a trajectory
generated in simulation.

The procedure for generating the covariance data is as follows: 
\begin{itemize}
\item Generate DNA strand \cite{Aksimentiev04}
\item Using NAMD \cite{NAMD} input: ``CG'' , set up configuration file to run for 10E10 timesteps
\item Equilibrate the system (\texttt{minimization on}) 
\item Using Matlab vercorrblock.m script,  calculate covariance of atoms over time and colorize matrix. 
(rearrange atoms to be next to their neighbors
in the covariance matrix)
\end{itemize}

\vspace{0.5in}

In order to perform parameter identification, equations were taken from from Crassidis \cite{Crassidis} Table 5.  Modeling
choices are explained.

$\dot{x}(t) = f(x(t)) $Initially there was no input to the system, and therefore no u(t), however, upon introduction of the 
estimation of $k_{ij}$ the state was augmented to be %\begin{eqnarray}
$Z = \left( \begin{array}{c}\eta \\
                           \dot{\eta}\\ k_{ij} \end{array}\right)$ \\
%\end{eqnarray}
$\tilde{y}_k = h(x_k) $ Measured y, $ \tilde{y}_k$ is the calculated centers of mass from the simulation. 
 The h matrix is the identity, as we have full state knowledge of the positions and velocities. 
 However, when we augment the parameters to be estimated $(k_{ij})$the system no longer has full state feedback. \\
$\hat{x}(t_0) = \hat{x}$ The system is initialized to the starting positions and velocities after energy is minimized\\
$P_0 = E{\hat{x}(t_0)\hat{x}^T(t_0)}}$ \\ This quantity was implemented as the expected value added to a randomly generated
number to prevent the matrix from being singular.
$K_k = P_k^-H_k^T[H_kP_k^-H_k^T + R_k]^{-1}$  various values of R were tried, usually about 0.01 \\
$\hat{x}^+ = \hat{x}_k^- + K_k[\tilde{y}_k - h]\\
P_k^+ = [I - K_kH_k]P^- \\
\dot{\hat{x}}(t) = f(\hat{x}) \\
\dot{P}(t) = F(\hat{x},t)P(t) + P(t)F^T(\hat{x},t) + G(t)Q(t)G^T(t)$ \\  F is defined below in function {\bf evalF.m}

\section{Results}
To view the matlab script see {\bf AA549.m} in the Appendix.  The script ran well and estimated all spring constants to be 
1.35 N/m, given equal masses.  Upon entering masses, which are on the order of $10^{-27}$, the problem became singular and
badly scaled.  This problem can be alleviated with careful choice of inital conditions and will be revisited for
cytosine and guanine connected by a hydrogen bond, adenine and thymine for a hydrogen bond, and various other
combinations of the nucleic acids for comparison.  These results can provide insight to the possibility
of a reduced order model for DNA kinetics.  DNA is ``programmable;'' the motions are predictable
(unlike that of protein) and the sequences can be chosen so that desired behaviors may arise.  Making computations more
tractable is one step toward characterizing the kinetics of DNA.  Currently there is no single widely accepted
representation.

%The potential energy is modeled the sum of spring potentials between each body. 

\newpage
\section{Appendix}\label{App}

%\subsection{Protein Data Bank File}
%\input{CGpdb}
%\subsection{Protein Structure File}
%\input{CGpsf}

%\subsection{CHARMM Force Field Parameterization}
%The CHARMM (Chemistry at HARvard Molecular Mechanics) force field parameterization is the most common way of expressing
%molecular mechanics.
%\input{charmm}
%\newpagey
%\input{centerofmass}

%\input{Amat}
%\subsection{Matlab routines}

%{\bf AA549.m} routine estimates the spring constant parameters.
%\input{function}

%{\bf evalF.m} routine calculates the linearized F matrix for each data point.
%\input{evalF}

%{\bf CytGuatraj.m} routine generates the trajectory of the centers of mass corresponding to the reduced order model.
%\input{CytGuatraj}

%{\bf getCOMdcd.m} is a subroutine is called by CytGuatraj.m to retrieve the center of mass from a dcd file.
%\input{getCOMdcd}

%{\bf velcorrblock.m} was used to produce the colorized covariance matrix.
%\input{velcorrblock}

%\subsection{Matrices}

%\input{Vmat}
%\input{matrixA}

\newpage
\bibliographystyle{hplain}
\bibliography{springreport05}

\end{document}