\documentclass[12pt]{article}
\usepackage[english]{babel}
\usepackage{amsmath, amssymb, bm, float}

\newcommand{\cvec}{\textbf{c}}
\newcommand{\evec}{\textbf{e}}
\newcommand{\fvec}{\textbf{f}}
\newcommand{\gvec}{\textbf{g}}
\newcommand{\honevec}{\textbf{h}_\textbf{1}}
\newcommand{\htwovec}{\textbf{h}_\textbf{2}}
\newcommand{\hvec}{\textbf{h}}
\newcommand{\lvec}{\textbf{l}}
\newcommand{\pvec}{\textbf{p}}
\newcommand{\qvec}{\textbf{q}}
\newcommand{\rvec}{\textbf{r}}
\newcommand{\svec}{\textbf{s}}
\newcommand{\vvec}{\textbf{v}}
\newcommand{\wvec}{\textbf{w}}
\newcommand{\xvec}{\textbf{x}}
\newcommand{\xnotvec}{\textbf{x}_\textbf{0}}
\newcommand{\xonevec}{\textbf{x}_\textbf{1}}
\newcommand{\xtwovec}{\textbf{x}_\textbf{2}}
\newcommand{\yvec}{\textbf{y}}
\newcommand{\zvec}{\textbf{z}}
%
\newcommand{\zerovec}{\textbf{0}}
\newcommand{\onevec}{\textbf{1}}
\newcommand{\gammavec}{\textbf{\gamma}}
\newcommand{\phivec}{\textbf{\phi}}
\newcommand{\psivec}{\textbf{\psi}}
\newcommand{\rhovec}{\textbf{\rho}}
\newcommand{\thetavec}{\textbf{\theta}}
\newcommand{\sensai}{\textsc{sensai}}
%
\newcommand{\Amat}{\textbf{A}}
\newcommand{\Bmat}{\textbf{B}}
\newcommand{\Cmat}{\textbf{C}}
\newcommand{\Fmat}{\textbf{F}}
\newcommand{\Imat}{\textbf{I}}
\newcommand{\Smat}{\textbf{S}}
\newcommand{\Tmat}{\textbf{T}}
\newcommand{\Vmat}{\textbf{V}}
\newcommand{\Wmat}{\textbf{W}}
\newcommand{\Xmat}{\textbf{X}}
\newcommand{\zeromat}{\textbf{0}}
\newcommand{\bzero}{\textbf{0}}
%
\newcommand{\reals}{\mathbb R}
\newcommand{\DD}{\displaystyle}
\newcommand{\norm}[1]{\left\lVert#1 \right\rVert_{2}}
%
%\newcommand{\xdim}{{\rm xdim}}
%\newcommand{\kdim}{{\rm kdim}}
\newcommand{\xdim}{{\rm M}}
\newcommand{\kdim}{{\rm K}}
\newcommand{\xonedim}{{\rm M_1}}
\newcommand{\xtwodim}{{\rm M_2}}


\begin{document} \title{SENSAI: README file}
\author{}
\date{}
\maketitle

%\pagestyle{empty}

\section{Software requirements}

\begin{enumerate}
    \item M{\sc atlab} 7.10.0 (R2010a)
    \item {\sf Maple} 14
\end{enumerate}

NOTE: This version of {\sensai} is no longer being maintained. It requires a M{\sc atlab} and {\sf Maple} interface, which is no longer maintained as of M{\sc atlab} version 7.12 (R2011a). If you have a version of M{\sc atlab} at or above 7.12, you must use one of the other versions of {\sensai}, which uses MuPAD. Unfortunately at this point, some of the examples do not work in the newer version, as MuPAD is not as effective as {\sf Maple} just yet.
  

\section{A Practical Guide to Using {\sc sensai}}

There are two basic ways to use {\sensai}, directly from the GUI or indirectly using a {\sf Maple} input file.  If the model contains up to four equations with six parameters, the user may run {\sensai} directly through the GUI (see section \ref{sec:direct}).  Alternatively, the user may construct and save his or her model in a separate {\sf Maple} worksheet and run it indirectly through the GUI (see section \ref{sec:indirect}). A glossary of {\sensai} terms appears in section \ref{sec:glossary}.

\newpage
\subsection{Using SENSAI directly}\label{sec:direct}

\begin{enumerate}
    \item Open {\sc Matlab}.
    \item Within {\sc Matlab}, change the directory to the location in which the {\sensai} software has been downloaded. We will call this {\tt SENSAI\_DIRECTORY}, e.g. C:/SENSAI/.
    \item Open the {\sc sensai} GUI by typing {\tt sensai} in the {\sc Matlab} command window.
    \item Select check box ``Input from GUI?".
    \item Select check box ``Iterated nonlinear map?" if the model is in the form of a map.
    \item Select check box ``Compute solutions only?" to compute solutions but not sensitivities and elasticities.
    \item Enter a job name to create a subfolder in which solutions will be stored.
    \item Enter the model, using [$\cdot$] to denote vector elements, i.e., {\sf Maple} syntax. (An example of a simple SIR model is provided.)
    \item Select ``Create {\sc Matlab} files using {\sf Maple}" which creates {\sc Matlab} files {\tt gvec.m}, {\tt dgvec\_dxvec.m}, {\tt dgvec\_dparam.m}, {\tt qoi.m} and {\tt dcp\_dparam.m} in {\tt SENSAI\_DIRECTORY}.  Wait until a popup box appears that says ``{\sc Matlab} files successfully created" before continuing.
    \item Enter parameter values.
    \item Select the desired plots in ``Plotting Information'' box.
    \item Select ``Execute {\sc Matlab} file created by {\sf Maple}".
    \item {\sensai} creates a directory {\tt SENSAI\_DIRECTORY/JOB\_NAME} in which it stores the {\sc Matlab} files {\tt gvec.m}, {\tt dgvec\_dxvec.m}, {\tt dgvec\_dparam.m}, {\tt qoi.m} and {\tt dcp\_dparam.m}, the input data, the plots and a binary file {\tt output.mat} containing {\tt xdim}, {\tt kdim}, {\tt tfinal}, {\tt t}, {\tt x}, {\tt p}, {\tt q}, {\tt dxdp}, {\tt dqdparam}, {\tt elxp} and {\tt elqp}.
    \item Results from a new run can be saved into another folder in {\tt SENSAI\_DIRECTORY} by changing the Job Name in the second box in the upper right corner of the {\sensai} GUI.  (E.g. ``run2").
%    \item Back to step 8 if the model equations or quantity of interest is new.  Back to step 9 if the only changes are in plotting data, parameter values, or initial conditions.
 \end{enumerate}

%\newpage
\subsection{Using SENSAI indirectly}\label{sec:indirect}

\begin{enumerate}
    \item Open {\sc Matlab}
    \item Within {\sc Matlab}, change the directory to the location in which the {\sensai} software has been downloaded. We will call this {\tt SENSAI\_DIRECTORY}, e.g. C:/SENSAI/.
    \item Open the {\sc sensai} GUI by typing {\tt sensai} in the {\sc Matlab} command window.
    \item Make sure the check box ``Input from GUI?" is \textit{unselected}.

    \item Create {\sf Maple} worksheet {\sf Model\_Name.mw} in what will be called {\tt WORKING\_DIRECTORY} by copying and modifying templates e.g.,
       \begin{description}
           \item[] {\sf Examples/ODE\_examples/SIR/SIR.mw}, or
           \item[] {\sf Examples/MAP\_examples/Caswell08/Caswell08.mw}.
       \end{description}
       It is important that the right-hand side is a vector named $g$ and the variables are in a vector named $x$, as {\sensai} will look for these specific variable names.

   \item Copy the path of {\tt WORKING\_DIRECTORY} into the box in the upper right hand within the GUI.
       \begin{description}
           \item[]  (e.g. {\sf C:/SENSAI/Examples/ODE\_examples/SIR/})
       \end{description}
   Note: \textit{The active directory within {\sc Matlab} must be the same one that contains the sensai.m program, i.e., {\tt SENSAI\_DIRECTORY}. This is especially important for step 8.}

    \item Execute {\sf Maple} worksheet {\sf Model\_Name.mw} which creates {\sc Matlab} files {\tt user\_equations.m}, {\tt user\_inputs.m}, {\tt user\_plotdata.m}, {\tt user\_qoi.m} and {\tt user\_parameters.m} in  {\tt WORKING\_DIRECTORY}.

    \item In the {\sensai}  GUI, select ``Create {\sc Matlab}  files using {\sf Maple}" which creates the files {\tt gvec.m}, {\tt dgvec dxvec.m}, {\tt dgvec dparam.m}, {\tt qoi.m}, and {\tt dcp dparam.m} in {\tt SENSAI\_DIRECTORY}.  Wait until a popup box appears that says ``{\sc Matlab} files successfully created" before continuing.

    \item Within {\sc Matlab}, control of the program is through the files {\tt user\_inputs.m} and {\tt user\_plotdata.m}, in {\tt WORKING\_DIRECTORY}.
        \begin{enumerate}
            \item Via {\tt user\_inputs.m}, you control parameter values, initial conditions, and the name of the subfolder of {\tt WORKING\_DIRECTORY} in which you wish to save your work (using ``JOB").
            \item Via {\tt user\_plotdata.m}, you control which solutions ($x$-values) to output and plot (using {\tt ilist}), and which parameters to have their sensitivities tested (using {\tt klist}).
        \end{enumerate}

    \item Within {\sc Matlab} in the GUI, select ``Execute {\sc Matlab} file created by {\sf Maple}". \\

          All plots of the solutions, their sensitivities and elasticities are saved in the subfolder {\tt WORKING\_DIRECTORY/JOB}. All solutions including sensitivities and elasticities are stored in the binary file \\ {\tt WORKING\_DIRECTORY/JOB/output.mat}.

%          \begin{enumerate}
%              \item To work with the numerical solutions, sensitivity values, and elasticities for manipulation or plotting, either work within {\sc Matlab} on the data by loading the file output.mat, or \dots
%              \item Use the automatically generated .eps files as plots for the run specified above.
%          \end{enumerate}

    \item Before performing another run using the {\sc sensai} GUI, return to the {\sc sensai} directory and enter the {\sc Matlab} commands\\ {\tt $>>$ close all, clear all} \\ to clear both plots and active memory before continuing. (Optional, but probably a good idea). Results from a new run can be saved into another folder in the {\tt WORKING\_DIRECTORY} by changing the name of ``JOB" in {\tt user\_inputs.m}, e.g. JOB = `run2'. (This can also be achieved by changing this line in the {\sf Maple} file \textit{but then the {\sf Maple} file must be re-run and you will return to step 7.})
        %{\tt $>>$ JOB\_NAME:= "run2";  //  Sets the folder name for the output.}

\end{enumerate}

\subsection{Glossary}\label{sec:glossary}
%
\subsubsection{.m-files used in {\sc sensai} }
\begin{table}[H]
\begin{center}
\begin{tabular}{||l||p{10cm}||}
\hline\hline
Name & Description \\
\hline\hline
{\tt sensai.m} & Main {\sc Matlab} file that controls the {\sensai} GUI\\
\hline
{\tt mm\_interface.m} & Main {\sc Matlab} file that interfaces {\sc Matlab} with {\sf Maple}\\
\hline
{\tt gtype.m} & Main {\sc Matlab} function that directs the solving, plotting, and storing of data \\
\hline
{\tt gvec.m} & Automatically generated {\sc Matlab} file that creates vector of right-hand side functions\\
\hline
{\tt dgvec\_dxvec.m} & Automatically generated {\sc Matlab} file that computes derivatives of right-hand side with respect to variables \\
\hline
{\tt dgvec\_dparam.m} & Automatically generated {\sc Matlab} file that computes derivatives of right-hand side with respect to parameters\\
\hline
{\tt qoi.m} & Automatically generated {\sc Matlab} file that constructs quantity of interest and its sensitivity information\\
\hline
{\tt r0\_matrix.m} & Automatically generated {\sc Matlab} file that computes $R_0$ and its sensitivity information\\
\hline
{\tt cp.m} & Automatically generated {\sc Matlab} file that constructs user-defined parameter and its sensitivity information\\
\hline\hline
\end{tabular}
\end{center}
\caption{Names and descriptions of .m-files used in {\sc sensai} }
\label{tab:glossary_mfile}
\end{table}


\subsubsection{Terms used in {\sc sensai}}

\begin{table}[H]
\begin{center}
\begin{tabular}{||l||p{10cm}||}
\hline\hline
Name & Description \\
\hline\hline
{\tt xdim} & Dimension of variables $\xvec$ \\
\hline
{\tt kdim} & Dimension of parameters $\pvec$ \\
\hline
{\tt tfinal} & Final time step\\
\hline
{\tt imap} & Logical that defines whether the system is a map or ode\\
\hline
{\tt solution\_only} & Logical that defines whether only the solutions are to be computed\\
\hline
{\tt x0} & Vector of initial conditions\\
\hline
{\tt p} & Vector of parameters\\
\hline
{\tt ilist} & List of a subset of variables to be included in plots\\
\hline
{\tt klist} & List of a subset of parameters to be included in plots\\
\hline
{\tt iplot} & Structure that determines which figures to plot \\
\hline
{\tt NextGen} & List of a subset of variables that describe infective classes\\
\hline
{\tt R0\_only} & Logical that defines whether $R_0$ is to be computed without its sensitivities\\
\hline
{\tt DIR} & {\tt WORKING\_DIRECTORY}\\
\hline
{\tt JOB} & Folder name where data from a run is stored\\
\hline\hline
\end{tabular}
\end{center}
\caption{Names and descriptions of terms used in {\sc sensai} }
\label{tab:glossary_terms}
\end{table}

\subsubsection{Variables in \textit{output.mat}}
\begin{table}[H]
\begin{center}
\begin{tabular}{||l|c|p{8cm}||}
\hline\hline
Name & Dimensions & Description \\
\hline\hline
{\tt xdim} & 1 & Number of variables\\
\hline
{\tt kdim} & 1 & Number of parameters\\
\hline
{\tt tfinal} & 1 & Final time step\\
\hline
{\tt p} & kdim$\times 1$ & Parameters\\
\hline
{\tt t} & $T\times 1$ & Vector of times at which the solutions and sensitivities are determined \\
\hline
{\tt x} & xdim $\times T$ & Matrix of solutions at all times $t(\cdot)$\\
\hline
{\tt q} & $T\times 1$ & Quantity of interest\\
\hline
{\tt r0} & 1 & The basic reproduction number $R_0$\\
\hline\hline
{\tt dxdp} & xdim $\times$ kdim $\times T$ & Sensitivity of variables with respect to parameters\\
\hline
{\tt dqdparam} & kdim$\times T$ & Sensitivity of the quantity of interest with respect to parameters\\
\hline
{\tt dR0dp} & kdim $\times 1$ & Sensitivity of $R_0$ with respect to parameters\\
\hline
{\tt elxp} & xdim $\times$ kdim $\times T$ & Elasticity of variables with respect to parameters\\
\hline
{\tt elqp} & kdim $\times T$ & Elasticity of quantity of interest with respect to parameters\\
\hline
{\tt elR0p} & kdim $\times 1$ & Elasticity of $R_0$ with respect to parameters\\
\hline\hline
{\tt cparam} & 1 & User defined parameter value\\
\hline
{\tt dxdc} & xdim $\times T$ & Sensitivity of variables with respect to user defined parameter\\
\hline
{\tt dqdc} & $T\times 1$ & Sensitivity of the quantity of interest with respect to user defined parameter\\
\hline
{\tt dR0dc} & 1 & Sensitivity of $R_0$ with respect to user defined parameter\\
\hline
{\tt elxc} & xdim $\times T$ & Elasticity of variables with respect to user defined parameter\\
\hline
{\tt elqc} & $T\times 1$ & Elasticity of quantity of interest with respect to user defined parameter\\
\hline
{\tt elR0c} & 1 & Elasticity of $R_0$ with respect to user defined parameter\\
\hline\hline
\end{tabular}
\end{center}
\caption{Names, dimensions, and descriptions of variables in \textit{output.mat}. }
\label{tab:glossary_output}
\end{table}
For notation purposes, $T$ is the dimension of time steps {\sc Matlab} uses to solve the system, not to be confused with {\tt tfinal}.  The value $T$ is \textit{not} a variable in {\sensai}.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Defining $R_0$ in {\sc sensai}}

For appropriate epidemiological models, {\sc sensai} is capable of automatically calculating the basic reproduction ratio, $R_0$, as defined by the Next Generation method.  However, since {\sc sensai} is not limited to infection modeling, a specific syntax is required so that {\sc sensai} recognizes if a
model is compatible with the definition of $R_0$.  The following guide instructs the user how to edit existing {\sf Maple} templates so that {\sc sensai} will
produce $R_0$ and its sensitivities with respect to parameters and initial conditions.

\begin{enumerate} \item Edit the {\sf Maple} templates to define your model equations.

    \item Indicate which equations describe the dynamics of infected classes.  Store these indices in the variable $NextGen$.  Note that $NextGen$ must be written in \textit{matrix} syntax.  Consider the following examples.
        \begin{enumerate}
            \item If the model includes three states, $S$, $I$, and $R$, in that order, $NextGen := matrix([2]);$.
            \item If the model has more than one equation describing an infected class, list them in the order they appear.  For example, if the model describes $S_1, I_1, R_1, S_2, I_2, R_2$ in that order, $NextGen := matrix([2,5]);$.
            \item \textit{If you do not wish to calculate $R_0$ for the model, define $NextGen = matrix([0])$, or let the first entry of $NextGen$ be 0.}
        \end{enumerate}

    \item If the model has four or more infected classes, you may want to consider computing $R_0$ without its sensitivities.  $R_0$ will be a very
    lengthy expression for such models, and the analytic derivatives will require a lot of time to compute.  If this is the case, define ``R0\_only"
    to be 1.  If you wish to calculate the sensitivities anyway, define ``R0\_only" = 0. \\

    (While running the {\sc sensai} GUI, you may encounter large delays in ``Create {\sc Matlab} files using {\sf Maple}" if R0\_only = 0.  If your patience has run thin, you must terminate the program through the task manager.  The emergency stop in M{\sc{atlab}} of CTRL+c in the command window will not work, as the computation of $R_0$ is done externally in a {\sf Maple} procedure call.)

    \item If the analytic expression for $R_0$ is already known, it may be faster (and more accurate if {\sf Maple} can not solve $R_0$) to use this expression of $R_0$ as the quantity of interest ($qoi$) instead of re-deriving the expression during the ``Create {\sc Matlab} files using {\sf Maple}" phase.

\end{enumerate}

\subsection{Automatic construction of $\mathcal{F,T,V}$}

\begin{enumerate}
    \item For ODEs, the Next Generation definition of $R_0 = \rho(FV^{-1})$ where $F = \dfrac{\partial \mathcal{F}_i}{\partial x_j}(x^{\star}) \quad 1 \leq i,j \leq m$ describes new infections and $V = \dfrac{\partial \mathcal{V}_i}{\partial x_j}(x^{\star}) \quad 1 \leq i,j \leq m$ describes transfer of existing infections, $x^{\star}$ is the disease-free equilibrium, the infected classes are $1,\dots,m$, and $\rho(\cdot)$ denotes the spectral radius operator.

    \item For maps, the Next Generation definition of $R_0 = \rho(F(I-T)^{-1})$, where $I$ is the $m\times m$ identity and $F$ and $-T$ are defined the same as $F$ and $V$ for ODEs, respectively.

\end{enumerate}

The following criteria are used by {\sc sensai} to determine the placement of each term.  If the terms of the model will
not be placed in the biologically correct vectors, {\sensai} fails to compute the Next Generation $R_0$.

\begin{enumerate}

\item If the term $X$ in an equation describing an infective class involves a state variable from a noninfectious class, $X \in \mathcal{F}$, unless
    the occurrence of the noninfectious state variable is part of a sum of all state variables (that is, the term is scaled by the total
    population).\\

\textbf{Important Note for MuPAD}:  Within the current release, MuPAD is not as effective as {\sf Maple} at analytical computations.  There may be
examples (like the Hantavirus model) where MuPAD can not account for the scaling by the total population.  This problem is due to limitations of the {\tt subs()} command and will hopefully be addressed in a future versions of MuPAD.\\

\item If the term $X$ in an equation describing an infective class does not involve any state variables and is only a parameter, product of
    parameters, or quotient of parameters, $X \in \mathcal{F}$.  However, if terms like these exist, the disease-free subspace will not be invariant, and the model will not have a valid Next Generation $R_0$.

\item Every other term $X$ that does not satisfy the above will be placed in $\mathcal{V}$ for ODEs, or $\mathcal{T}$ for maps.

\end{enumerate}

\subsection{Potential problems with $R_0$}

There are some examples in which the Next Generation construction of $R_0$ is not valid, or is not compatible with {\sc sensai}.  %The following are
possible problems the user might encounter when trying to define $R_0$.


\subsubsection{Potential problems with ODE models.}  For $R_0$ to be valid, the model must satisfy the conditions of Theorem 2 of \cite{vanden2002}, specifically:

 \begin{description}
  \item[(A1)] The fecundity matrix $F$ is nonnegative.  %This is part of assumption (A1) of \cite{vanden2002}.
  \item The transition matrix $V$ is nonsingular.  %{\tt SJT: Huh?? This can occur if an equation is in the model as a placeholder, but the right-hand side is identically 0.  This state must be removed from the system for $R_0$ to be valid.}
  \item[(A4)] The disease-free subspace is invariant.  That is, infection can enter a disease-free population through a nonzero component in a state
      that is identified as disease-free.  This can occur in models with background infection rates, or in models where the infective classes are
      not identified properly.  %This is assumption (A4) of \cite{vanden2002}.
  \item[(A5)] The equilibrium is  asymptotically stable in the absence of disease.  That is, if $\mathcal{F} = 0$, all eigenvalues of the
      Jacobian of the full system evaluated at $x^{\star}$ have negative real part.  %This is assumption (A5) of \cite{vanden2002}.
 \end{description}

Notice that assumptions (A2) and (A3) of \cite{vanden2002} for ODE models are not automatically checked by {\sensai}.  These assumptions must be verified by the user, but are usually true.

\subsubsection{Potential problems with map models.}  For $R_0$ to be valid, the model must satisfy the conditions of Theorem 2.1 in \cite{Allen2008}, specifically:

 \begin{enumerate}
  \item The fecundity matrix $F$ is nonnegative.
  \item The transition matrix $T$ is nonnegative.
  \item The transition matrix $T$ is nonsingular.  %{\tt SJT: Huh?? Singular $T$ can occur if an equation is in the model as a placeholder, but the right-hand side is identically 0.  This state must be removed from the system for $R_0$ to be valid.}
  \item The transition matrix $T$ is  asymptotically stable.  That is, $\rho(T) < 1$.
  \item The equilibrium is asymptotically stable in the absence of disease.  That is, $\rho(C) < 1$ where $C$ is the Jacobian of the right-hand side of the noninfectious states.
 \end{enumerate}

For map models, the assumption of a unique DFE is not checked by {\sensai}, nor is the condition that $F+T$ is irreducible.  These should also be checked by the user to ensure a valid $R_0$.  It is difficult to check both of these conditions, but again, for most models, $F+T$ is irreducible based on the structure of $T$ having a nonzero main diagonal and a sub-diagonal and the structure of $F$ having a nonzero top row.

%There are a number of reasons for any of the problems in lists 1 and 2 to occur.  Perhaps the model does not have a valid Next Generation construction of
%$R_0$.  If this is the case, some alternative means to calculate $R_0$ should be sought, if desired.  Alternatively, {\sc sensai} may not be able to
%recognize which terms describe new infections and belong in $\mathcal{F}$ and which terms describe transfer of existing infections and belong in
%$\mathcal{V}$ or $\mathcal{T}$.



\section{Template Examples}

\subsection{MAP Examples}

\subsubsection{Caswell 08}

This model can be found in \cite{Caswell2008}.  It involves two stages, juveniles ($x_1$) and adults ($x_2$) and the map
%
\begin{equation*}
\xvec(t+1) = \begin{pmatrix} \sigma_1 (1-\gamma) & f \\ \sigma_1 \gamma & \sigma_2 \end{pmatrix} \xvec(t),
\end{equation*} %
where the juvenile survival
$\sigma_1(\xvec) = \tilde{\sigma} e^{-\textbf{e}^\textrm{T}\xvec}$, where $\textbf{e}$ is a vector of ones, $\sigma_2$ is the adult survival, $\gamma$ is
the maturation probability, and $f$ is the adult fertility.  The parameter values given are $(f, \gamma, \tilde{\sigma}, \sigma_2) = ( 0.25, 1/15, 0.98,
0.95)$.

The main purpose of this model is to verify that {\sensai} gives the same equilibrium and sensitivities to those mentioned in the paper.  This model also serves as a template for MAP examples.

\subsubsection{Hantavirus}

This model is from \cite{Allen2008}.  It involves susceptible and infected male and female rodents. The map is
%
\begin{equation*}\label{eq:hantavirus}
\left.
\begin{gathered}\begin{aligned}
S_m(t+1) &= \Big[ \frac{B}{2}+e^{-\beta_mI_m(t) - \beta_f I_f(t)}S_m(t)\Big]D(N)\\
I_m(t+1) &= \Big[
(1-e^{-\beta_mI_m(t)-\beta_f I_f(t)})S_m(t) + I_m(t)\Big]D(N)\\
S_f(t+1) &= \Big[ \frac{B}{2}+e^{-\beta_f I_m(t)-\beta_f I_f(t)}S_f(t)\Big] D(N)\\
I_f(t+1) &= \Big[ (1-e^{-\beta_f I_m(t)-\beta_f I_f(t)})S_f(t)+I_f(t)\Big]D(N)
\end{aligned}\end{gathered}
\right \}
\end{equation*}
%
where the logistic growth is scaled by
%
\begin{equation*}
D(N) = \frac{K}{K+(b/2)N},
\end{equation*}
%
where $K$ is the carrying capacity and $N$ is the total population. The birth function is the harmonic mean
%
\begin{equation*}
B(N_m,N_f) = \frac{2bN_mN_f}{N},
\end{equation*}
%
where $N_m = S_m+I_m$ is the total number of males, $N_f = S_f+I_f$ is the total number of females, and $b>0$ is the average litter size, $k$ is the number of contacts that result in an infection, and $\beta_m$ and $\beta_f$ are the infection rate constants of males and females, respectively.  Parameter values are not provided in \cite{Allen2008}, but some reasonable values are $K= 1000$, $\beta_f=0.09$, $\beta_m=0.9$, and $b=6$.

The main purpose of this model is to verify that {\sensai} gives the same value of $R_0$ as provided by \cite{Allen2008}.  This example exhibits two common practices in model formulation which {\sensai} performs extremely well \textit{in this version}.  First, every term is scaled by the total population $N$. Second, the infection rate is given by a probability of an infection occurring.  Even with these two difficulties, {\sensai} is able to identify which terms belong to $\mathcal{F}$ and $\mathcal{V}$ correctly to produce a valid $R_0$.  (Unfortunately, for the users with a later version, {\sensai} is not able to handle the first issue due to limitations in MuPAD's {\tt subs()} command.)

\subsubsection{Six Stage Genetics}

This model is based on \cite{Caswell2008}, adding partially dominant selection based on viability and two alleles.  The full model is explained in detail in \cite{Tavener2011}.

\subsection{ODE Examples}

\subsubsection{SIR}

This model is a typical SIR model with logistic growth.
%
\begin{equation*}\label{eq:SIR_log}
\left .
\begin{gathered} \begin{aligned}
\frac{dS}{dt} &= rN
\left( 1-\frac{N}{K}\right) - \beta SI - \delta S, \quad \\
\frac{dI}{dt} &= \beta SI - \gamma I - \mu I - \delta I, \\
\frac{dR}{dt} &= \gamma I - \delta R,
\end{aligned} \end{gathered}
\right \}
\end{equation*}
%
where $N = S+I+R$ is the total population at any time $t$, $r$ is the per capita growth rate, $K$ is the carrying capacity, $\beta$ is the infection rate, $\delta$ is the natural death rate of the species, $\gamma$ is the recovery rate, and $\mu$ is the disease specific death rate.  Some reasonable parameter values are $r = 0.5$, $K=1000$, $\beta = 0.1$, $\delta = 0.2$, $\gamma$ = 0.02, and $\mu = 0.1$.

The main purpose of this model is to demonstrate the calculation of $R_0$  for a standard ODE infection model which can be easily verified by hand.  This model also serves as a template for ODE examples.

\subsubsection{SI (Indirect Transmission)}

This model is an SI model that involves indirect transmission of the infection, and is
%
\begin{equation*}\label{eq:background_si}
\left.
\begin{gathered}\begin{aligned}
\frac{dS}{dt} &= rN\left(1-\frac{N}{K}\right)\\
\frac{dI}{dt} &= \beta -\gamma I
\end{aligned}\end{gathered}
\right\}
\end{equation*}
%
where $N(t) = S(t) +I(t)$ is the total population, $r$ is the per capita growth rate, $K$ is the carrying capacity, $\beta$ is the background transmission probability, and $\gamma$ is the recovery rate.

A detailed explanation of this model can be found in \cite{Mikucki2012}.  Parameter values may be chosen as  $r= 0.5$, $\beta = 0.8$, $\gamma = 0.02$, and $K = 1000$.

The main purpose of this model is to show that models with a background (indirect) transmission of the disease through the environment or some alternative source do not have a valid $R_0$.

\subsubsection{Plague}

This model is given by \cite{Buzby2008}.  The first three classes are the SIR classes of rats, $N$ is the average number of fleas living on a rat, and $F$ is the number of free infectious fleas that are searching for a new host. The system of ordinary differential equations is

\begin{equation*}
\left.
\begin{gathered}\begin{aligned}
\dot{S}_R &= r_RS_R\left(1-\frac{T_R}{K_R} \right) + r_RR_R(1-p)-d_RS_R -
\beta_R\frac{S_R}{T_R}F(1-e^{-aT_R}) \\
\dot{I}_R &= \beta_R\frac{S_R}{T_R}F(1-e^{-aT_R})-(d_R+m_R)I_R\\
\dot{R}_R &= r_RR_R\left(p-\frac{T_R}{K_R}\right) + m_Rg_RI_R -d_RR_R\\
\dot{N} &= r_FN\left(1-\frac{N}{K_F}\right) + \frac{d_F}{T_R}F(1-e^{-aT_R})\\
\dot{F} &= (d_R+m_R(1-g_R))I_RN-d_FF
\end{aligned}\end{gathered}
\right\}
\end{equation*}
%
where $T_R = S_R+I_R+R_R$ is the total size of the rat population, $r_R$ is the net rat reproduction rate, $K_R$ is the rat carrying capacity, $p$ is the proportion of offspring that inherit the disease, $d_R$ is the natural rat death rate, $\beta_R$ is the transmission rate from rats to fleas, $m_R$ is the rate that rats leave the infected class, $g_R$ is the fraction of rates that become resistant, $a$ is the searching efficiency of the fleas, $r_F$ is the net flea reproductive rate, $d_F$ is the natural flea death rate, and $K_F$ is the flea carrying capacity.

\subsubsection{Dengue}

This model is given by \cite{Garba200811}.  It is an SEIR model that describes the dynamics of dengue fever, an infection carried by a vector, mosquitos. The equations are
%
\begin{equation*} \label{eq:dengue}
\left. \begin{gathered} \begin{aligned}
\frac{dS_H}{dt} &= \Pi_H-\lambda_HS_H-\mu_HS_H\\
\frac{dE_H}{dt} &= \lambda_HS_H-(\sigma_H+\mu_H)E_H\\
\frac{dI_H}{dt} &= \sigma_HE_H-(\tau_H+\mu_H+\delta_H)I_H\\
\frac{dR_H}{dt} &= \tau_HI_H-\mu_HR_H\\
\frac{dS_V}{dt} &= \Pi_V-\lambda_VS_V-\mu_VS_V\\
\frac{dE_V}{dt} &= \lambda_VS_V-(\sigma_V+\mu_V)E_V\\
\frac{dI_V}{dt} &= \sigma_VE_V-(\mu_V+\delta_V)I_V
\end{aligned} \end{gathered}
\right \}
\end{equation*}
%
where $\lambda_H = \frac{C_{HV}}{N_H}(\eta_VE_V+I_V)$ is the human infection rate, $\lambda_V = \frac{C_{HV}}{N_H}(\eta_HE_H+I_H)$ is the vector infection rate, and $N_H = S_H+E_H+I_H+R_H$ is the total human population. The other parameter values and interpretations can be found in Table \ref{tab:dengue}.

\begin{table}[H]
\begin{center}
\begin{tabular}{|| l || c | p{8cm} ||}
\hline\hline Parameter & Numerical Value & Interpretation \\
\hline\hline $\mu_H$    & $0.0195$ & $1/\mu_H$ is average human lifespan \\
$\sigma_H$     & $0.5300$ & Rate of transfer of exposed to infected humans \\
$\Pi_H$        & $10$     & Human recruitment rate \\
$\delta_H$     & $0.9900$ & Disease specific human death rate\\
$\eta_H$       & $0.9900$ & Infectiousness factor of exposed to infected humans\\
$\tau_H$       & $0.2000$ & Human recovery rate \\
\hline $\mu_V$ & $0.0140$ & $1/\mu_V$ is average vector lifespan\\
$\sigma_V$     & $0.2000$ & Rate of transfer of exposed to infected vectors\\
$\Pi_V$        & $30$     & Vector birth rate\\
$\delta_V$     & $0.0057$ & Disease specific vector death rate\\
$\eta_V$       & $0.9800$ & Transmissibility factor of exposed to infectious vectors\\
\hline $C_{HV}$   & $0.038$  & Infection rate of mosquitoes\\
\hline\hline
\end{tabular}
\end{center}
\caption{Parameter Values for Dengue Model.}
\label{tab:dengue}
\end{table}

An interesting initial condition for the model is $\textbf{x}_\textbf{0} = (\frac{\Pi_H}{\mu_H}, 0, 0, 0, \frac{\Pi_V}{\mu_V}, 0, 200)$, which will show that even though $R_0 <1$, infection may still persist in the population.

\newpage
\subsubsection{Typhoid}

This model is given by \cite{Bailey1980}.  It is a 9-stage SIR type model with the variables as described in Table \ref{tab:typhoid_variables} and the model equations are

\begin{equation*}\label{eq:typhoid}
\left.
\begin{gathered}\begin{aligned}
\dot{x_1} &= -(\rho_{12}+\rho_{13})x_1y + \rho_{41}x_4+\rho_{51}x_5+\rho_{61}x_6+\rho_{81}x_8+\rho_{91}x_9-\mu x_1+\mu\\
\dot{x_2} &= \rho_{12}x_1y - (\rho_{23}+\rho_{24} +\rho_{25}+\mu)x_2 + \rho_{32}x_3\\
\dot{x_3} &= \rho_{13}x_1y - (\rho_{32}+\rho_{34}+\rho_{35}+\mu)x_3 + \rho_{23}x_2\\
\dot{x_4} &= \rho_{24}x_2 +\rho_{34}x_3 + \rho_{54}x_5-(\rho_{41}+\rho_{45}+\rho_{46}+\rho_{48}+\mu)x_4\\
\dot{x_5} &= \rho_{25}x_2 + \rho_{35}x_3+\rho_{45}x_4-(\rho_{51}+\rho_{54}+\rho_{58}+\mu)x_5 \\
\dot{x_6} &= \rho_{46}x_4-(\rho_{61}+\rho_{67}+\rho_{68}+\mu)x_6\\
\dot{x_7} &= \rho_{67}x_6 - \mu x_7\\
\dot{x_8} &= \rho_{48}x_4 + \rho_{58}x_5 + \rho_{68}x_6 - (\rho_{81}+\rho_{89}+\mu)x_8\\
\dot{x_9} &= \rho_{89}x_8 - (\rho_{91}+\mu)x_9.
\end{aligned}\end{gathered}
\right \}
\end{equation*}

Parameter values are given in Table \ref{tab:typhoid_variables}.
\begin{table}[H]
\begin{center}
\begin{tabular}{||l|l||}
\hline\hline
Name & Description \\
\hline\hline
$x_1$ & susceptibles \\
$x_2$ & incubating noninfectious \\
$x_3$ &  incubating infectious \\
$x_4$ &  sick infectious\\
$x_5$ &  sick noninfectious\\
$x_6$ & temporary carrier\\
$x_7$ & permanent carrier\\
$x_8$ & short resistance\\
$x_9$ & long resistance\\
\hline\hline
\end{tabular}
\end{center}
\caption{Variables for Typhoid Model.}
\label{tab:typhoid_variables}
\end{table}
%

This example demonstates {\sensai}'s ability to implement a large system and compute $R_0$ effectively.  In this model, states $x_2$ through $x_7$ are considered infective, so the next generation matrix is a $6\times 6$ matrix with analytical (not numerical) components.  Parameter values can be found in Table \ref{tab:typhoid}.

\begin{table}[H]
\begin{center}
\begin{tabular}{||l||c|l||}
\hline\hline Parameter & Numerical Value & Interpretation \\
\hline\hline $\rho_{12}$ & $8.43381\times10^{-3}$ & Infection rate \\
$\rho_{13}$ & $8.51900\times10^{-5}$  & Infection rate \\
\hline $\rho_{23}$ & $2.85720\times10^{-3}$ & Transfer rate from state of incubation\\
$\rho_{24}$ & $6.78585\times10^{-2}$ & Transfer rate from state of incubation\\
$\rho_{25}$ & $7.14300\times10^{-4}$ & Transfer rate from state of incubation\\
\hline $\rho_{32}$ & $7.14300\times10^{-4}$ & Transfer rate from state of incubation\\
$\rho_{34}$ & $6.42870\times10^{-2}$  & Transfer rate from state of incubation\\
$\rho_{35}$ & $6.42870\times10^{-3}$  & Transfer rate from state of incubation\\
\hline $\rho_{41}$ & $3.46000\times10^{-3}$ & Transfer rate from state of sickness\\
$\rho_{45}$ & $3.46000\times10^{-3}$ & Transfer rate from state of sickness\\
$\rho_{46}$ & $3.46000\times10^{-3}$ & Transfer rate from state of sickness\\
$\rho_{48}$ & $2.40124\times10^{-2}$ & Transfer rate from state of sickness\\
\hline $\rho_{51}$ & $3.46000\times10^{-3}$ & Transfer rate from state of sickness\\
$\rho_{54}$ & $6.92000\times10^{-3}$ & Transfer rate from state of sickness\\
$\rho_{58}$ & $2.40124\times10^{-2}$ & Transfer rate from state of sickness\\
\hline $\rho_{61}$ & $1.11100\times10^{-3}$   & Transfer rate from temporary carrier\\
$\rho_{67}$ & $3.33300\times10^{-3}$ & Transfer rate from temporary carrier\\
$\rho_{68}$ & $6.66600\times10^{-3}$ & Transfer rate from temporary carrier\\
\hline $\rho_{81}$ & $2.74000\times10^{-4}$  & Transfer rate from short resistance\\
$\rho_{89}$ & $2.46600\times10^{-3}$ & Transfer rate from short resistance\\
\hline $\rho_{91}$ & $2.74000\times10^{-4}$  & Transfer rate from long resistance\\
\hline $\mu$  & $5.48000\times10^{-5}$ & Overall birth and death rate\\
\hline\hline
\end{tabular}
\end{center}
\caption{Parameter Values for Typhoid Model according to Bailey, all units are $days^{-1}$.  Note that each $\rho_{ij}$
is a transfer rate from state $i$ to state $j$.}
\label{tab:typhoid}
\end{table}



\newpage

%\bibliographystyle{ecology}
\bibliographystyle{plain}
\bibliography{bibliography}


\end{document}
