%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%	chap2.tex
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\psdraft
\chapter{MODEL DEFINITIONS AND DEVELOPMENT}\label{ch:model}   

In this chapter we will illustrate the use of references, equations, equation
arrays, tables, and some figures. Below we have numerous references.

The first biophysical mechanisms of bursting related to pancreatic
$\beta$-cell electrical activity were developed by Atwater {\it et al}
(1980)~\cite{atwater:nob80}. Chay and Keizer (1983)~\cite{chay:mmm83} used this
work to create the first 'minimal' (biophysical) mathematical model, based on
the Hodgkin-Huxley model. Since then there have been a large number of
$\beta$-cell models~\cite{chay:ebp88},~\cite{himmel:tse87},~\cite{keizer:asp89},
\cite{sherman:eob88},~\cite{keizer:bea91},~\cite{chay:ecc90},~\cite{smolen:svi92}
and other cellular models exhibiting bursting behavior.
 
In addition to studying bursting models in the context of pancreatic
$\beta$-cells (or other endocrine cells), there has been great interest in
modeling nerve cells~\cite{deschenes:tbm82},~\cite{harris-warrick:mmb87},
\cite{wong:agh81},~\cite{wang:fbf99}, \cite{bose:ass00},~\cite{kepecs:acb00}.
Nerve cells, or neurons, are specialized for rapid
electrical signaling over long distances. Neurons, structurally consist of a
soma, axons and dendrites. Whereas, axons and dendrites are unique to neurons,
the soma resembles the cells in other organs. The basic mechanism of
communication between neurons is the transmission of action potentials along
axons. Many neurons exhibit bursting behavior. It has been speculated that
bursting is a process that makes communication across synapses (between neurons)
reliable~\cite{lisman:bun97},~\cite{izhikevich:bun03},~\cite{hoppensteadt:tci98},~\cite{wang:fbf99}.

I will prove the Contraction mapping principle of Banach-Cacciopoli on page~\pageref{thm:fixed_pt}.
% An example of a page reference


        \subsection{Equations}      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

An example of an equation,
\begin{verbatim}
\begin{equation}     \label{eq:HEQN3}
    h(u) = \beta (u-\alpha),
\end{equation}
\end{verbatim}
yields the equation shown below in equation \eqref{eq:HEQN3}.  
By the way we label equations, figures, tables, etc by $\backslash$label\{name\}
and reference the labels as: $\backslash$ref\{name\}, $\backslash$eqref\{name\}, or $\backslash$pageref\{name\}.

An equation array, {\bf eqnarray}, as in \eqref{eq:UEQN3}-\eqref{eq:ZEQN3}, is
used to align a number of equations, the alignment is on the \&=\&. Any symbol
may be used in place of the '=' sign here. Any number of equations may be used;
{\bf $\backslash$notag} or {\bf $\backslash$nonumber} omits one equation number,
{\bf eqnarray*} omits all.
\begin{verbatim}
\begin{eqnarray}
    LHS1 &=& RHS1
       \label{eq:1} \\
    LHS2 &=& RHS2
       \label{eq:2} \\
    ...
\end{eqnarray}
\end{verbatim}
The aforementioned Pernarowski polynomial model is:
\begin{eqnarray}
    \frac{du}{dt} & = & f(u) - w - z \ ,
       \label{eq:UEQN3} \\
    \frac{dw}{dt} & = & \frac{1}{\tau} ( g(u) - w) \ ,
       \label{eq:WEQN3} \\
    \frac{dz}{dt} & = & \varepsilon ( h(u) - z) \ ,
       \label{eq:ZEQN3}
\end{eqnarray}
where $0<\varepsilon \ll 1$ is a parameter, $f$ and $g$ are cubic
polynomials, and
\begin{equation}     \label{eq:HEQN3}
    h(u) = \beta (u-\alpha),
\end{equation}
where $\alpha, \beta$ are parameters. 


%% The symbol  $\reals$ is defined in mystyle.sty 
A simple array is shown and not numbered in the $P$ definition below. Use
'$\backslash$notag' to prevent numbering. Note, as already mentioned, you must
escape the special symbol '\{' as ``$\backslash$\{``.

A two dimensional manifold in $ \reals^4$ defined as follows: \label{'SL'}
\begin{equation}
    {\cal S}_{L} = \{ (u,w,x,y) : x+\gamma y = G(u), w=g(u), u<u_{-}\} ,
    \label{eq:SS_L}
\end{equation}
where $u_{-}$ is the value of $u$ at the limit point.
If ${P:\reals^4~\rightarrow~\reals^4}$ is the projection operator:
\begin{equation}
   P=\left( {\begin{array}{c} u \\ w\\x\\y \end{array}} \right) =
       \left( {\begin{array}{c} 0 \\ 0\\x\\y \end{array}} \right) ,
   \notag
\end{equation}
then we define the projection onto the $(x,y)$-plane, $P({\cal S}_L)=S_L$:
\begin{equation*}  %%% the {equation*} is equivalent to using \notag
    S_{L} = \{ (u,w,x,y) : x+\gamma y = G(u), u<u_{-},u=0,w=0 \} .
\end{equation*}

More examples,
\begin{equation}
  I_x = \bar g_x n^p_x (v-v_x) = g_x(v-v_x) .
  \notag   %%%%%% use this command when you do not want an equation number
\end{equation}
Then we may rewrite the current balance equation as the sum of the fast currents
$ I_f$ and the two slow currents mentioned above: 
\begin{equation}
  C_m \frac{dv}{dt}  =  - \sum_f I_f -\bar g_x x (v-v_x)  -\bar g_y y (v-v_x) .
     \notag
\end{equation}
On the fast $t$ time scale the dynamics of (FULL) %use \quad or qquad for large
				%spacing in math mode. 
is governed by the Fast Subsystem (FS) obtained by letting $\varepsilon = 0$,
\begin{eqnarray}
 \frac{du}{dt} & = &  f(u) - w -z \quad , \quad z \equiv x + \gamma y \>,
  \label{eq:FS1} \\
 \frac{dw}{dt} & = & g(u) - w \>,
     \label{eq:FS2} 
\end{eqnarray}


We may rewrite the orthogonality condition in matrix form as
\begin{eqnarray}
  \langle Q,b \rangle &=& \int^T_0 Q(t,\alpha)^{\top}  b(t,\alpha) dt \\
  &=& \int^T_0 Q(t,\alpha)^{\top} \left ( G(U(t,\alpha),t)-  D_\alpha U(t,\alpha) \frac{d\alpha}{d\tau} \right ) dt \\
  &=& \int^T_0 Q(t,\alpha)^{\top} G(U(t,\alpha),t) dt -  \frac{d\alpha}{d\tau}
  \int^T_0 Q(t,\alpha)^{\top} D_\alpha U(t,\alpha) dt 
  \label{eq:thm3.5b}\\
  &=& 0  \label{eq:thm3.5c}.
\end{eqnarray}

The Jacobian of (FULL)
\begin{equation}
   D\vec{\bf  F} =
  \left[ {\begin{array}{cccc} f'(u) & -1 &  -1 & - \gamma \\
      g'(u) &- 1& 0&0 \\
      \varepsilon \frac{h_1'(u)}{\tau_1}&0 & -\frac{ \varepsilon}{\tau_1}& 0 \\
      \varepsilon \frac{h_2'(u)}{\tau_2}&0&0&- \frac{ \varepsilon}{\tau_2}
  \end{array}} \right] ,
  \label{eq:jacobian}
\end{equation}
has the characteristic polynomial
\begin{equation*}
    P(\lambda) = det\left(D\vec{\bf F} - \lambda I \right) .
\end{equation*}

Example using {\bf align} environment instead of eqnarray and the
$\backslash$intertext command.
\begin{align}
   \varepsilon \frac{du}{d\tilde t} &= f(u) - w - x - \gamma y \ , \\
   \varepsilon \frac{dw}{d\tilde t} &= g(u) - w \ , \\ 
  \frac{dx}{d\tilde t} &= \frac{h_1(u) - x}{\tau_1} \ , \\ 
  \frac{dy}{d\tilde t} &= \frac{h_2(u) - y}{\tau_2} \ . \\
\intertext{Setting $\varepsilon = 0$ we get, to leading order,}
    x + \gamma y &= f(u) - g(u) \equiv G(u) \ , 
    \label{eq:ss1a}\\
   w &= g(u) \ , 
    \label{eq:ss2a}\\ 
  \frac{dx}{d\tilde t} &= \frac{h_1(u) - x}{\tau_1} \ , 
  \label{eq:ss3a}\\ 
  \frac{dy}{d\tilde t} &= \frac{h_2(u) - y}{\tau_2} \ . 
  \label{eq:ss4a}
\end{align}
For $u<u_-$, equations \eqref{eq:ss1a}-\eqref{eq:ss4a} define the (SS) of
(FULL). Solutions of \eqref{eq:ss1a}-\eqref{eq:ss4a} yield trajectories on
${\cal S}_L$ which are leading order approximations to the silent phase of
(FULL).

We allow the three parameters ($\lambda$, $x_f$, $T$) to
vary in this run with $x_0$ fixed until $\lambda = 1$:
\begin{eqnarray}
    \frac{dx}{d\tau} &=& T F_1(x,y) = T \left( \frac{\beta_1(u_{\LB}(x+\gamma y)-\alpha_1)-x}
	 {\tau_1} \right) ,\\ 
    \frac{dy}{d\tau} &=& T F_2(x,y) = T \left(\frac{\beta_2(u_{\LB}(x+\gamma
         y)-\alpha_2) - y}{\tau_2}\right) ,
\end{eqnarray}
where $u_{\LB}$, is the variable $u$ restricted to the lower branch.


        \subsection{Tables}      %%%%%%%%%%%%%%%%%%

% Example of a table:
% Table caption can be selected as paragraph style or centered style
% (for an inverted pyramid title). Use \oldstylecaptiontrue (paragraph)
% or \oldstylecaptionfalse (centered) to select the style.
%
We present several different tables in this section.
In Table~\ref{table:var} we outline the model variations we will consider.
\begin{table}[hb]
\begin{center}
\begin{tabular}{|c ||c|}
\multicolumn{2}{c}{\bf Activation/ Inactivation }\\
\hline
Variation & Slow variables are \\ 
\hline \hline
$(+,+)$ & Both activation \\ \hline
$(-,-)$ & Both inactivation  \\ \hline
$(+,-)$  & Mixed activation/inactivation \\ \hline
\multicolumn{2}{c}{\vspace{-6pt} }\\
\multicolumn{2}{c}{\bf Time Constants }\\
\hline
  & Degenerate case reduces\\
 $\tau_1 = \tau_2$ & to the One-Slow variable model.\\ \hline
 & Non-degenerate case \\
$\tau_1 \neq \tau_2$ &  here $(+,-)$ is different from $(-,+)$. \\ \hline \hline
\end{tabular}
\caption[Two slow variable model variations (overview)]{Possible variations in the form of the
  two-slow variable model. The activation/inactivation variation refers to the signs of
  the parameters  ($\beta_1,\beta_2$). }
\label{table:var}
\end{center}
\end{table}


I better check, but I don't think the grad office allows captions above tables
(read: you better check).
\begin{table}[ht]
\centering
\caption{\label{timing1} PDE solve times, $15^3+1$ equations\label{pde.tab1}}
\begin{tabular}{||l|l|l|l|l|l||}\hline
Precond. & Time & Nonlinear & Krylov
& Function & Precond. \\
 & & Iterations & Iterations & calls & solves \\ \hline
None & 1260.9u & 3 & 26 & 30 & 0  \\
 &(21:09) & & & &  \\ \hline
FFT  & 983.4u & 2  & 5  & 8  & 7 \\
&(16:31) & & & & \\ \hline
MILU & 629.7u & 3  & 11 & 15 & 14 \\
& (10:36) & & & & \\ \hline
\end{tabular}
\end{table}


A slightly more complicated interior grid is shown in Table~\ref{table:lambda1}.
\begin{table}[h]
\begin{center}
\begin{tabular}{|c|c|c|c|c|}
\hline
Case  & sign($q_2$) & sign($q_1$) & sign($q_0$) & Equilibria \\ \hline \hline
$(+,+)$ & +  & +   & +   & stable \\ \hline
\multirow{2}{11mm}{$(-,-)$} 
     & \multirow{2}{3mm}{+}  & +   & + & stable  \\ \cline{3-5}
        &    & $+/-$ & $-$ & saddle \\ \hline
$(+,-)$ & +  & $+/-$ & $+/-$ & HB possible \\ \hline
$(-,+)$ & +  & $+/-$ & $+/-$ & HB possible \\ \hline
\end{tabular}
\caption[Dependence of (SS) equilibria on parameters.]{Dependence of (SS)
  equilibria $\vec{\bf x}_e \in {\cal S}_L$ on the parameters $(\beta_1,
  \beta_2)$, where $(q_2, q_1, q_0)$ are the coefficients of the
  $O(\varepsilon^2)$ term in the characteristic polynomial of (FULL) with
  $\lambda_0=0$. In the table $+/-$ indicates that both signs are possible}
\label{table:lambda1}
\end{center}
\end{table}

\clearpage


	\section{Figures}   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




% Figures in LaTeX will go on the bottom of the same page, or the top of
% the next page, but never before the first reference. All figures must
% be referenced. The syntax is below.
%
%       \begin{figure}[x]    % x = b, t, h, p for bottom, top, here(if possible) and
%              page. Can also use [ht] or other combos.
%              May use H if you use the float package; put the float here period.
%       ...
%       \caption[short form]{My title - full caption}  % Captions are below the figure!
%       \end{figure}
%
%       The 'short form' is what will appear in the list of tables.
%       The portion in { } is what appears below the figure.
%
% We are using the package pstricks below. This allows us to utilize native TeX
% fonts, and gives us more control over the figures.
% Uncomment the \psgrid line to overlay a positioning grid

A range of parameters, for which \eqref{eq:UEQN3}-\eqref{eq:ZEQN3} yields square wave
bursts such as those in Figure~\ref{fig:burst}, will be discussed in a
subsequent section.

I will have many more examples of figures in Chapter~\ref{ch:numerics}.

\begin{figure}[h]
   \pspicture(0,0)(13.1,8.9)
   %\psgrid[subgriddiv=1,griddots=10,gridlabels=7pt](0,0)(13,9)
   \rput[bl](0,0){\includegraphics[angle=-90,width=14cm]{burst.ode.ps}}
   \rput[l](0.9,5.5) {$u$}
   \rput[l](8.3,0) {$t$}
   %\rput[b]{90}(8,5){\psframebox{STUFF}}   %% the '{90}' here rotates the text - But only shows up rotated in the .ps or .pdf
   %\rput[b]{-34}(8,1){Some comment}               %% here the '{-34}' rotates the text -34 degrees, as above NOT in the .dvi file.
   \endpspicture
   \caption[Square wave bursting]{Square wave bursting - solution of
     \eqref{eq:UEQN3}-\eqref{eq:ZEQN3}}
   \label{fig:burst}
\end{figure}


% These graphs were created by gnuplot. For simple graphs this is a
% great utility.  For example, if you want a sin curve in your thesis
% try the following:
%
% (terminal window): gnuplot
% (in gnuplot):
%                 set terminal latex
%                 set output "foo.tex"
%                 plot sin(x)
%                 quit
%
\begin{figure}[ht]       % Place it on the bottom of page
\centering              % Put \label{} into \caption.
\inputpicture{figs/sin.tex}
\caption{\label{fig:sincurve} A simple $\sin$ curve}
\end{figure}


Next we compare this approximation of the Hopf bifurcation
to one numerically computed using xppaut~\cite{ermentrout:91} on
(FULL) with 
$(\beta_1, \beta_2) = (5.035, -2.858)$ associated with the Hopf point
found above for $u=-1.16$.
\begin{figure}[H]
  \begin{center}
    \pspicture(0,0)(14,8.1)
%    \psgrid[gridcolor=cyan,subgriddiv=1,griddots=10,gridlabels=7pt](0,0)(14,9)
    \rput[bl](0.5,0){\includegraphics[angle=-90,width=12cm]{HB_run4.ps}}
    \rput[l](7.3,0.1){ $x$}
    \rput[l](1.3,4.3){ $y$}
    \rput[l]{-27}(5.5,2.6){\footnotesize $x+\gamma y = z_{-}$}
    \psline[linestyle=dotted,arrowsize=0.2]{->}(6.65,5.55)(6.6,5.6)
     \psline[linecolor=red,linewidth=.5pt](2.2,4.85)(3.3,4.22)
    \endpspicture
    \caption[Slow flow near a Hopf point.]{Slow flow near the numerically
    approximated Hopf point. Trajectories of (FULL) with parameter values
    $(\beta_1, \alpha_1, \beta_2 \alpha_2,\tau_1, \tau_2) = (5.0348,-1,-2.8583, -0.5, 1, 0.3)$}
    \label{fig:HB_run4}
  \end{center}
\end{figure}
 


\clearpage % dump figures where they below



        \subsection{Theorems, Lemmas, etc}      %%%%%%%%%%%%%%%%%%{Still waiting...}





\begin{thm}    %%%%%%%%%%%%%%% THEOREM %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\label{Malkin thm}
Consider a T-periodic dynamical system of the form
\begin{equation}
   \dot X= F(X,t) + \varepsilon G(X,t), \quad X \in \reals^n,
   \label{eq:thm3.1}
\end{equation}
and suppose that the unperturbed system, $\dot X= F(X,t)$, has a k-parameter
family of T-periodic solutions
\begin{equation}
   X(t) = U(t,\alpha),
   \label{eq:thm3.1X}
\end{equation}
where $\alpha=(\alpha_1, ... ,\alpha_k)^{\top} \in \reals^k$ is a vector of
independent parameters, which implies that the rank of the $n \times k$ matrix
$D_\alpha U$ is k. Suppose the adjoint linear problem
\begin{equation}
   \dot Q_i = - \{D F(U(t,\alpha))\}^{\top}Q_i
\end{equation}
has exactly k independent T-periodic solutions $Q_1(t,\alpha),..., Q_k(t,\alpha)
\in \reals^m$. Let Q be the matrix whose columns are these solutions such that
\begin{equation}
   Q^{\top} D_\alpha U = I
   \label{eq:thm3.1i}
\end{equation}
where I is the identity $k \times k$ matrix. Then the perturbed system
\eqref{eq:thm3.1} has a solution of the form
\begin{equation}
  X(t) = U(t,\alpha(\varepsilon t)) + O(\varepsilon),
   \label{eq:thm_expansion}
\end{equation}
where
\begin{equation}
   \frac{d\alpha}{d\tau} = \frac{1}{T} \int^T_0 Q(t,\alpha)^{\top}  G(U(t,\alpha),t) dt ,
   \label{eq:thm3.2}
\end{equation}
where $\tau = \varepsilon t$ is slow time.
\end{thm}

\pf Let ......            %%% add text to a reference:
The Fredholm alternative (see Hale~\cite[pg 146]{hale:ode80} for a proof of this
variation) implies that the linear $T$-periodic equation~(\#) has a
unique solution if and only if the orthogonality condition.....

.... as was to be shown. {\scriptsize $_\blacksquare$} %need amssymb package
%%%%%%%%%%% END PROOF



%%%%%%%%%%%%%%% THEOREM %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{thm}[Contraction mapping principle of Banach-Cacciopoli]
\label{thm:fixed_pt}
Let $\phi:D(\phi)\rightarrow \reals $ and assume that there exists $M\subset
D(\phi)$ which is closed and nonempty such that
\begin{itemize}
  \item[(i)] $\phi:M \rightarrow M$ ; and
  \item[(ii)] $\phi$ is k-contractive, ie,  there exists $k \in [0,1)$ such that\\
  $\|\phi(x)-\phi(y)\|\leq k \|x-y\| \quad$ for all $x,y \in M$.
\end{itemize}
Then the following are true:
\begin{itemize}
     \item[(a)] Existence and uniqueness. The function $\phi$ has exactly one
     fixed point $\bar x \in M$.
     \item[(b)] Convergence of the iteration method. For each $x_0 \in M$, the
     sequence $\{x_n\}^{\infty}_{n=1}$ constructed by $x_{n+1}=\phi(x_n)$ converges to the unique
     fixed point  $\bar x$.
%      \item[(c)] Error estimates. {\bf should I mention this??}.
%      \item[(d)] Rate of convergence. {\bf should I mention this???}.
\end{itemize}
\end{thm}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% COMMENTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% I don't ever prove that k=1.
% But since U is a 2 dimensional system and since Q are the T-periodic solns to
% the linearized system AND since the limit cycles are attractive there must be
% one eigen direction that flows to the limit cycle. As in we can have both
% directions T-periodic




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{comment}
words here
\begin{equation}

\end{equation}


\begin{eqnarray}

\end{eqnarray}

\notag
%%%Combination to do page references
 \label{pg:compress}
\pageref{pg:compress}

TO rotate things in a figure (Won't show up in the .dvi file):
\rput[b]{90}(8,5){\psframebox{STUFF}}
\rput[b]{-34}(8,1){Sanoe}

\end{comment}