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


\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{Simon Tavener, Michael Mikucki, Stewart Field and Michael Antolin} 
\date{} 
\maketitle

%\pagestyle{empty}

\section{Software requirements}

\begin{enumerate}
    \item M{\sc atlab} 7.13.0 (R2010a)
\end{enumerate}

NOTE:  This version of {\sensai} uses the new M{\sc atlab} software called MuPAD.  If you have an earlier version of M{\sc atlab}, you must use one of the other (un-maintained) version of {\sensai}.  If you have M{\sc atlab} 7.12 (R2011a), you must use the {\sensai} created specifically for this version.  If you have M{\sc atlab} 7.10 (R2010a) or 7.11 (R2010b), you must use the {\sensai} which also requires {\sf Maple}.  Unfortunately at this point, some of
the examples do not work in the newer versions (M{\sc atlab} 7.12 and higher), 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 MuPAD input file.  If the model contains up to four equations with six parameters and does not need to be saved, 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 MuPAD worksheet and run it indirectly through the GUI (see section \ref{sec:indirect}).

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

\begin{enumerate}
    \item Open M{\sc atlab}
    \item Within M{\sc atlab}, change the directory to the location in which the {\sensai} software has been downloaded. We will call this the {\tt SENSAI\_DIRECTORY}, e.g. C:/SENSAI/.
    \item Open the {\sc sensai} GUI by typing {\tt sensai} in the M{\sc atlab} 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 data using [$\cdot$] to denote vector elements, i.e., MuPAD syntax. (An example of a simple SIR model is provided.)
    \item Select ``Create M{\sc atlab} files using MuPAD" which creates {\tt gvec.m}, {\tt dgvec\_dxvec.m}, {\tt dgvec\_dparam.m}, {\tt qoi.m} and {\tt dcp\_dparam.m}.  Wait until a popup box appears that says ``M{\sc atlab} files successfully created" before continuing.
    \item Select the desired plots in ``Plotting Information'' box.
    \item Select ``Execute M{\sc atlab} file created by MuPAD".
    \item {\sensai} creates a directory {\tt SENSAI\_DIRECTORY/JOB\_NAME} in which it stores the M{\sc atlab} 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, kdim, tfinal, t, x, p, dxdp, q, dqdparam, elxp} and {\tt elqp}.
    \item Results from a new run can be saved into another folder in the {\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}

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

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

    \item Create MuPAD worksheet {\sf user\_input.mu} in what will be called the {\tt WORKING\_DIRECTORY} by copying and modifying templates e.g.,
       \begin{description}
           \item[] {\sf Examples/ODE\_examples/SIR/SIR.mu}, or 
           \item[] {\sf Examples/MAP\_examples/Caswell08/Caswell08.mu}.
       \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 the {\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 M{\sc atlab} must be the same one that contains the sensai.m program, i.e. the {\tt SENSAI\_DIRECTORY}.
   This is especially important for step 8.}

    \item Execute MuPAD worksheet {\sf user\_input.mu} which creates M{\sc atlab} files {\tt user\_equations.m}, {\tt user\_inputs.m}, {\tt
        user\_plotdata.m}, {\tt user\_qoi.m} and {\tt user\_parameters.m} in the {\tt WORKING\_DIRECTORY}.

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

    \item Within M{\sc atlab}, control of the program is through the files {\tt user\_inputs.m} and {\tt user\_plotdata.m}, in the {\tt
        WORKING\_DIRECTORY}.
        \begin{enumerate}
            \item Via {\tt user\_inputs.m}, you control parameter values, initial conditions, and the name of the folder 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 ``ilist"), and which parameters to             have their sensitivities tested (using ``klist").
        \end{enumerate}

    \item Within M{\sc atlab} in the GUI, select ``Execute M{\sc atlab} file created by MuPAD". \\
    
          All of the plots of the solutions, sensitivities, and elasticities specified in the run of the model, and a file with all of the outputs from the model (output.mat) will be saved in the {\tt WORKING\_DIRECTORY} in a subfolder named by the variable string ``JOB."
                
          \begin{enumerate}
              \item To work with the numerical solutions, sensitivity values, and elasticities for manipulation or plotting, either work within M{\sc atlab} 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, within , return to the {\sc sensai} directory and enter the M{\sc atlab} 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 MuPAD file \textit{(but then the MuPAD 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}

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

\newpage

\begin{table}[htbp] 
\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 JOB} & Folder name where data from a run is stored\\
\hline
{\tt DIR} & {\tt WORKING\_DIRECTORY}\\
\hline\hline
\end{tabular}
\end{center}
\caption{Names and descriptions of terms used in {\sc sensai} }
\label{tab:glossary_terms}
\end{table}

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




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
\section{How to Define $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 MuPAD templates so that {\sc sensai} will
produce $R_0$ and its sensitivities with respect to parameters and initial conditions.

\begin{enumerate} \item Edit the MuPAD 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 M{\sc atlab} files using MuPAD" 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 MuPAD procedure call.)

    \item If the analytic expression for $R_0$ is already known, it may be faster (and more accurate if MuPAD 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 M{\sc atlab} files using MuPAD" phase.

\end{enumerate}

\subsection{Possible 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$.

\begin{enumerate}

\item \textbf{ Problems with ODE models.}  For $R_0$ to be valid, the model must satisfy the conditions of Theorem 2 of \cite{vanden2002}.  Recall 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.

 \begin{enumerate}

  \item The fecundity matrix $F$ is not nonnegative.  This is part of assumption (A1) of \cite{vanden2002}.
  \item The transition matrix $V$ is singular.  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 The disease-free subspace is not 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 The equilibrium is not asymptotically stable in the absence of disease.  That is, if $\mathcal{F} = 0$, there is an eigenvalue of the
      Jacobian of the full system evaluated at $x^{\star}$ that has a positive real part.  This is assumption (A5) of \cite{vanden2002}.

 \end{enumerate}

\item \textbf{Problems with map models.}  For $R_0$ to be valid, the model must satisfy the conditions of Theorem 2.1 in \cite{Allen2008}.  Recall 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.

 \begin{enumerate}

  \item The fecundity matrix $F$ is not nonnegative.
  \item The transition matrix $T$ is not nonnegative.
  \item The transition matrix $T$ is singular.  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 The transition matrix $T$ is not asymptotically stable.  That is, $\rho(T) \geq 1$.
  \item The equilibrium is not asymptotically stable in the absence of disease.  That is, $\rho(C) \geq 1$ where $C$ is the Jacobian of the
      right-hand side of the noninfectious states.

 \end{enumerate}

\end{enumerate}

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.  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}$.  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}$.  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}

\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$). 
% 
\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. 
% 
\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 the Maple version}.  First, every term is scaled by the total population $N$. Second, the infection rate is given by a probability of a non-infection not occurring.  Unfortunately, while {\sensai} handles the second issue, it is not able to handle the first.  This is due to limitations in MuPAD's {\tt subs()} command.

As an illustration, consider the following term.  Suppose there are two states $x_1$ and $x_2$ and the term in consideration is $p_1x_2 / (p_1 +
x_1+x_2)$.  If $x_2$ is the infectious state, then this term belongs in $\mathcal{V}$, since the only appearance of a non-infective state is in the
denominator, which is just a scaling by the total population.  This illustrates what occurs for the Hantavirus model with the scaling $D(N)$, with $p_1$
acting as $K$.  While the {\sf Maple} command {\tt subs(term, $x_1 + x_2 = N$)} converts this term to $p_1x_2 / (p_1 + N)$, the identical MuPAD command
does not.  If $p_1$ were not in the denominator, the MuPAD command would work as desired.  Unfortunately, the computation of $R_0$ for models with the
scaling of the total population is not available in the MuPAD version of {\sensai}, although it is available for the {\sf Maple} version.

\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 that it is a standard ODE infection model with a simple $R_0$ that 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. 
% 
\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.

\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}[htbp] 
\begin{center} 
\begin{tabular}{||l||c|l||}
\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{table}[htbp] 
\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}
% 
\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*}

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}.

\newpage

\begin{table}[htbp] 
\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}
