\section{Introduction}\label{mainobj}
\subsection{Stochastic dynamics for experiments on bacterial genetic evolution }


We focus our theoretical study on stochastic models for the genetic evolution of bacterial populations. Such models can be
recast as discrete-time Markov chains in large-dimensional space of histograms for bacteria with different genotypes.
These models are also called the ``locked-box'' models; they have been developed to describe daily evolution 
of finite-size populations with multiple genotypes \cite{azencoop,rice,heger}.

Our analysis is motivated by the long-term laboratory experiments on genetic evolution of \emph{Escherichia coli}. In these experiments (see \cite{cooper1,heger,len94a,Fox2015,Barrick2013,Barrick2010,Cooper2001,Deatherage2017,Woods2011,Levy2015,Plank1979}), on day $n$, the current cell population $pop_n$ has a roughly fixed large size $N$ and grows freely until nutrient exhaustion. One then extracts (by dilution) a random sample of approximately $N$ cells, which constitutes the next day population $pop_{n+1}$. Widely used genetic evolution models for such experiments implement a succession of  ``daily'' evolutions comprising of three steps - (i) growth phase of a fixed duration, (ii) Poisson distributed mutations, and 
(iii) random selection of a sub-sample of fixed size $N$.

We focus our attention on rare events where the frequency of some intermediate-strength genotype can become unusually large.
Such events are called fixations, and it is important from the biological point of view to understand the evolutionary paths which
can lead to emergence of a large-sub-population of mutants with a non-dominant fitness.

For the stochastic dynamics described above, computationally usable formulas for fixation probabilities (see \cite{heger,rice}) have been either limited to $g =2$ genotypes, have assumed very small ``selective advantages", or have been restricted to extremely small mutation rates to make sure that at most only one new genotype emerges before fixation. Estimation techniques based on intensive simulations have been implemented (see for instance \cite{azencoop,heger,Gordo2013,Illinworth2012,Simon2013,Ross1996}) to evaluate mutation rates and selective advantages from experimental data.



\subsection{Rare events and genetic evolution}
For large bacterial populations (e.g. $N \approx 10^6$), many potential genotypes fixations become \emph{rare events} with probabilities far too small to be correctly evaluated by simulations. Large Deviations approach is a natural tool to study rare events for such stochastic dynamics.
Large Deviations results have been obtained for the trajectories of a wide range of vector-valued stochastic processes (see, for example, \cite{dembo,varadhan,ventcell,azen77,azen80,azenLDT2012}), including Markov chains, Gaussian processes, stochastic differential equations with a small diffusion coefficient, etc. Here we extend the Large Deviations approach to discrete-time Markov chains
modeling genetic evolution of bacterial populations.

The numerical applicability of rare events analysis has not often been exploited in concrete models of cell populations experiments. Previous large deviations results for stochastic population evolution have involved theoretical asymptotic studies such as \cite{champ01,champ02,champ6,champ8}. These papers have studied very general asexual population evolution in the space of phenotypic traits vectors. In such models, Darwinian evolution of asexual populations is driven by birth and death rates, which are themselves dependent on phenotypic traits. These traits are approximately transmitted to offspring with rare but important variations due to genes mutations. Competition for limited resources forces a permanent or roughly periodic selection. Mutations are assumed to be rare enough so that most of the time, only one currently dominant trait vector can coexist with the traits vectors of emerging mutants. 

Evolutionary models considered here are quite different because they combine two random steps - mutations and dilution. Therefore, a combination of these two effects can lead to observed rare events and it is important to understand contributions of these two events to the experiments in evolutionary dynamics \cite{Gordo2012,Desai2007,Barrick2013,Korona1994,Peng2017}.


\subsection{Scope of our large deviations analysis}
Here, we focus our large deviations study on developing rigorous but computationally implementable algorithms to quantify rare genetic events for a discrete Markov chain of  genetic histograms $H_n$ in the following context.

We consider sequences $pop_n$ of large bacterial populations submitted  daily to growth, mutations, and random selection of $N$ cells constituting $pop_{n+1}$. We characterize $pop_n$ by the histogram $H_n \in \mbb{R}^g$ of genotype frequencies in $pop_n$. For large $N$, and any fixed time horizon $T$, we develop a Large Deviations framework for the space $\Omega(T)$ of all histograms trajectories $\mbold{H}(1,T) = [H_1\,H_2 \ldots H_T ]$.

For large $N$, the random histogram trajectory $\mbold{H}(1,T) = [H_1 \,H_2\ldots H_T ]$ starting at an initial (given) histogram $H_1$ has a high probability of being close to the deterministic ``mean trajectory'' $\mbold{Y}(1,T) = [Y_1\, Y_2 \ldots Y_T]$ recursively defined by $Y_{n+1} = \mbb{E}[H_{n+1} \; | \; H_n = Y_n ]$ with $Y_1 = H_1$.
For any histogram trajectory $\mbold{w}$, we derive an \emph{explicit} large deviations rate functional $\lambda(\mbold{w}) \geq 0$ such that $Prob(\mbold{w})$ is roughly equivalent to $ \exp(- N \lambda(\mbold{w}))$. For any subset $A$ of $\Omega(T)$, the rate functional 
$$\Lambda(A) = \inf_{\mbold{w} \in A} \; \lambda(\mbold{w}) \geq 0$$
then verifies, under weak conditions on $A$,
$$
\lim _{N \to \infty} \frac{1}{N} \log P(A) = - \Lambda(A)
$$
When $\Lambda(A) > 0$, the event $ ( \mbold{w} \in A )$ is then a ``rare event'' with probability vanishing exponentially fast as $N \to \infty.$
For any given initial and terminal histograms $H$ and $G$, we characterize for large $N$ the \emph{most likely evolutionary paths} from $H$ to $G$, by minimizing $\lambda(\mbold{w})$ over all $\mbold{w}$ starting at $H$ and ending at $G$ at some finite time $T$. We derive a new, \emph{explicit} second-order reverse recurrence equation satisfied by these optimal paths. This essentially solves in $\mbb{R}^g$ the discretized Hamilton-Jacobi-Bellmann PDEs verified by the extremal histogram trajectories. 




\section{Stochastic model for bacterial evolution experiments}\label{stochmodel}
To model the main features of random bacterial evolutions, we focus on a class of Markov chains often used in this context \cite{azencoop,rice,heger}. 
The finite set of distinct genotypes is of size $g$ and denoted $\{1, 2, \dots, g \}$. Cells of genotype $j$ are called \emph{j-cells} here and have exponential growth rate $ f_j > 0$ called the \emph{fitness} of genotype $j$. We detail below the three phases (growth, mutation, random selection) comprising a single evolutionary cycle and present the space of histograms that quantifies the concentration of $j$-cells for each $j \in \vset{1,\dots,g}.$  Then, we will describe the transition probability associated to a single cycle.

\subsection{The space of genetic histograms} \label{histograms}
For any population $pop$ involving at most $g$ specific genotypes, the genetic state of $pop$ will be described by the vector $H = \left( H(1), H(2), \ldots, H(g) \right)$ of the $g$ genotype frequencies within $pop$, and $H$ will be called the \emph{genetic histogram} of this population.  The set $\mathcal{H} \subset \mbb{R}^g$ of all potential histograms is the compact, convex set of all vectors $H$ with coordinates $0 \leq H(j) \leq1$, and $\sum_{j=1}^{g} H(j) = 1$.

The genetic histogram of a population of size $N$ actually belongs to the finite grid $\mathcal{H}_N$ of \emph{N-rational histograms} $H$ such that all coordinates of $N H$ are integers. Note that $card(\mathcal{H}_N) \leq (N+1)^g$.\\
The set $\mathcal{H}$ will be endowed with the \emph{$L_{\infty}$-distance}
$$
\vnorm{H - G} = \max_j\,{ \abs{H(j) - G(j)}} \quad \text{for all } \; H,G\in\mathcal{H}.
$$ 
Since our large deviations rate functions are strongly sensitive to the smallest \emph{positive} coordinate of histograms, we make the following definition.
\begin{mydef}
For any $H \in \mathcal{H}$, define its \emph{support} $spt(H)$ and its \emph{essential minimum} $b(H) > 0$ by 
\begin{equation}\label{spt(H)b(H)}
spt(H) = \vset{j \, | \,H(j)>0} \quad \text{and} \quad b(H) = \min_{j \in spt(H)} \; H(j) 
\end{equation}

\end{mydef}
Note that $b(H) \geq 1/N$ for all $H \in \mathcal{H}_N$.  Stochastic genetic evolution after $T$ evolutionary cycles is then described by the (random) histogram trajectory $\mbold{H} \equiv \mbold{H}(1,T) = [H_1\, H_2 \ldots H_T]$ of arbitrary duration $T$ where $H_n \in \mcal{H}_N$ is the histogram on day $n$ for $1 \leq n \leq T$. We now present formally the three phases that occur during a ``daily'' cycle.  We note that we consider a day in this context a single evolutionary cycle, which often takes $8-12$ hours in real time.  For the sake of clarity in notation when referencing histograms and histogram trajectories, we will use bold letters to describe histogram trajectories where the time duration is often suppressed.  Non-bolded letters will typically denote histograms unless otherwise stated.
\subsection{Path space of histograms sequences}
Genetic evolution will be modeled as a sequence of  daily cycles indexed by day $n$. To simplify our evolution model, we have artificially split each daily population growth into two successive steps: deterministic growth with no mutations, followed by random mutations assumed to occur simulataneously. This rough simplification is only introduced here to facilitate the presentation of our large deviations analysis and of its numerical implementation. A companion paper studies the more realistic stochastic dynamics where random mutations can occur at any time during daily growth. In that more sophisticated model, the large deviations theory is similar but more technical, and will be outlined in another paper.

At the beginning of day $n$, the current population $pop_n$ will always have the same large size $N$ and will be identified by its genetic histogram $H_n \in \mbb{R}^g$. The day $n$ cycle involves three successive phases to generate $pop_{n+1}$, and its genetic histogram $H_{n+1}:$ \smallskip \\
{\bf Phase 1:}  pure deterministic growth where the $j$-cells have multiplicative growth factor $F_j > 1$ . \smallskip \\
{\bf Phase 2:}  random mutations occurring simultaneously at fixed very small mean rates. \smallskip \\
{\bf Phase 3:}  random sub-sample of fixed size $N$, which constitutes $pop_{n+1}$.\smallskip \\
\subsection{Phase 1 - Deterministic growth with no mutations} Call $H_n $ the genetic histogram of $pop_n$. Let $D$ be the \emph{fixed} duration of this pure growth phase during which no mutation is allowed to occur. The impact of Phase 1 on $pop_n$ is to multiply the initial size $N H_n(j)$ of each $j$-cell colony by the \emph{growth factor} $F_j =\exp ( D f_j ) > 1$, where $f_j$ is the “fitness” coefficient of $j$-cells. Let $F=( F_1, \ldots, F_g )$, and denote the scalar (dot) product in $\mbb{R}^g$ as $\langle \cdot , \cdot \rangle$ . We always order the $F_j$ by \emph{strictly increasing} fitness $F_1 < F_2 < \ldots < F_g$. Genotype $1$ is called the \emph{ancestor}, and genotype $g$ is the fittest. The difference $f_j - f_1 >0$ is often called the \emph{selective advantage} of genotype $j$ over the ancestor genotype $1$.

In actual experiments on bacterial evolution (see\cite{cooper1,heger}), the daily multiplicative growth factors $F_j$ for ``observable'' genotypes $1\leq j \leq g$ typically range from 20 to 300, and "detectable" selective advantages over the ancestor are typically larger than $0.01$.
Based on the description of the deterministic growth phase, we can now calculate the histogram at the end of Phase 1.  To this end, for any $u \geq 0,$ denote $ \ceil{u}$ the smallest integer greater than or equal to $u$ (the ceiling function). At the end of Phase 1, each $j$-cell colony reaches a size $N F_j H_n(j)$. The terminal population size given $H_n = H$ is then 
$
N_{ter} = \ceil{N \langle F, H \rangle} = \ceil{ \sum_j N F_j H(j) }.
$
So after Phase 1 with $H_n = H$, the population genetic histogram $\varphi(H)$ is given by 
$
\varphi(H) = \ceil{ N F_j H(j) }  /  \ceil{ N \langle F, H \rangle}. 
$
One can naturally approximate $\varphi(H)$ by $\Phi(H)$ where $\Phi : \mathcal{H} \to \mathcal{H} $ is defined for all genotypes $j$ by
\begin{equation} \label{phi}
\Phi_j(H) = F_j H(j) / \langle F, H \rangle.
\end{equation} 
Precise but elementary computations indeed show that for $N>20$ and all $H \in \mathcal{H}$ that
\begin{equation} \label{intervalphi}
\vnorm{\Phi(H) - \varphi(H)} \leq  \frac{4 F_g }{N F_1}.
\end{equation}

This completes Phase 1.  Now, we will focus on Phase 2, the mutation phase.  This phase is a bit more intricate since the feasibility of mutations will naturally be constrained by the population size parameter $N$ and the histogram $\Phi(H_n)$ after Phase 1.  Namely, this will force the mutation matrices (described below) to live in a set of $N$-rational matrices.  As was done above in approximating the histogram after deterministic growth using \eqref{phi}, we would like to make approximations so that the $N$-rationality condition may be dropped.  Therefore, we will perform these approximations in detail, which will greatly aid large deviations analysis further on.



\subsection{Phase 2 - Random mutations} 
At the end of the growth phase on day $n$, Phase 2 allows all random mutations to occur \emph{simultaneously}. Since we only have a finite amount $j$-cells for each genotype $j$ after Phase 1 and some mutations may be impossible, there are natural constraints that limit the possible number of cells that can mutate from one genotype $j$ to another genotype $k$, which we describe below.  Then, we describe the relevant probability distributions associated to Phase 2.
\subsubsection{Mutation matrices and mutation rates} \label{lincons} 
Let $R_n(j,k)$ be the random number of $j$-cells mutating into $k$-cells for $k \neq j$, and set $R_n(j,j) = 0$ for all $j$. The random matrix $R_n $ belongs to the set $\mathcal{Z}$ of $g \times g$ matrices $Z$ with integer coefficients $Z_{j,k} \geq 0$ and all $Z_{j,j}= 0$.

The $g \times g$ matrix $M$ of \emph{mean mutation rates} with entries $M_{j,k} \geq 0$ will be fixed and of the form $ M = m Q$, where $m>0$ is a \emph{very small, fixed mutation scale} process parameter shared by all cell types. The transfer matrix $Q$ is also fixed with entries $Q_{j,k} \geq 0$ and $Q_{j,j}= 0$ for all $j$. Within the $j$-cell colony, the fixed global mean emergence rate for mutants will be 
$
m_j= m \; \sum_k Q_{j,k}= \sum_k M_{j,k}.
$ 
For bacterial populations, these $m_j$ typically range from $10^{-9}$ to $10^{-6}.$  Thus, in order to simplify formulas further on, we will systematically assume that $m$ and $Q$ verify $m \leq 10^{-6}$ and $m_j \leq 10^{-6}$ for all $j$.

Given $H_n= H$, the total number of mutants $\sum_k R_n(j,k)$ emerging from $j$-cells must be inferior to the number $N F_j H(j)$ of $j$-cells after growth. As specified below, this imposes $g$ \emph{linear constraints} on $R_n/N$, and forces $R_n/N$ to live in a convex set of $N$-rational matrices $K_N(H)$ which we present below.
\subsubsection{Mutations constraint set $K(H)$} 
We begin with a few definitions describing the space where the mutation matrices reside along with the linear constraint space.
\begin{mydef}[$N$-Rational Matrices]
\label{def:mut-mat}
A matrix $A$ will be called \emph{N-rational} if $N A$ has non-negative integer coefficients. The set $\mathcal{Z}(N)$ of $N$-rational $g \times g$ matrices includes all \emph{standardized} random mutation matrices $R_n / N$.  
\end{mydef}

\begin{mydef}[Constraint Sets]
\label{KNconstr}
For each histogram $H$ and each $j$, define $K(j, H)$ as the set of all $g \times g$ matrices $r$ with non-negative coefficients verifying for all $k$
\begin{equation} \label{Kj}
\begin{cases}
\sum_{k=1}^g  r_{j,k} <  Fj H(j), & \text{when } \quad H(j) > 0;\\
r_{j,k}  =  0, & \text{when }\quad H(j) = 0;\\
r_{j,k}  =0 , & \text{when }\quad Q_{j,k} = 0.
\end{cases}
\end{equation}
Then, define two convex sets of matrices by $K(H)  = \bigcap_{j=1}^g \; K(j, H)$ and $K_N(H) = \mathcal{Z}(N) \cap K(H).$
\end{mydef}
Note that $K_N(H) \subset K(H)$ by definition.  Due to \eqref{Kj}, we have $ r_{j,k} \leq F_g$ and $card( K_N(H) )\leq [ (1+F_g)(N+1) ]^{g^2}$ for all $H\in \mathcal{H}$, $r \in K(H)$, $j,k$, and $N$. We state and prove a quick lemma describing the density of $K_N(H)$ in $K(H).$

\begin{lemma} \label{density}
Fix $a>0$. For any $N> \frac{g^2}{a F_1}$, any $H \in \mathcal{H}$ with $b(H) \geq a$, and any $r \in K(H)$, there is an $N$-rational matrix $s = \; [ s_{j,k} ]$ such that 
\begin{equation}\label{ss}
s \in K_N(H) \subset K(H); \qquad spt(s) = spt(r); \qquad|| s - r || \leq g/N 
\end{equation}
\end{lemma}

\begin{proof}
In this proof only, denote $\floor{z}$ as the largest integer less than or equal to $z$ (the floor function).  Fix any $j \in spt(H)$, and select any $k^*= k^*(j)$ such that $r_{j, k^*} = \max_ k \; r_{j,k}.$  For any $k$, if  $r_{j,k*} < g/N ,$ define $s_{j,k}$ by
\begin{equation} \label{defsjk}
s_{j,k} =
\begin{cases}
0, & r_{j,k} = 0; \\
1 / N, & 0 < r_{j,k} < 1/N; \\
 \floor{N r_{j,k}} / N, & 1/N \leq r_{j,k}.
\end{cases} 
\end{equation}
If $g/N \leq r_{j,k*}$, define 
\begin{equation} \label{defsjk*}
s_{j,k*} = [ N \; r_{j,k*} ] / N - (g-1 )/N.
\end{equation}
Definitions \eqref{defsjk} and \eqref{defsjk*} imply $\vnorm{s - r} \leq g/N$ and $spt(s) = spt(r)$.
Let $S(j) = \sum_k \, s_{j,k} $. To prove $s \in K(j, H)$, we only need to show that $S(j) < F_j H(j)$ for all $j \in spt(H)$, which involves two cases.  Thus, let $j \in spt(H).$\smallskip\\
{\bf Case 1:} Suppose $r_{j,k*} < g/N.$  Since $F_j H(j) > F_1 b(H) > F_1 a$, impose $N > g^2/a F_1$ to force $F_j H(j) > g^2/N $. Definition \eqref{defsjk} gives $ s_{j,k} < g/N$ for all $j, k$, and hence $S(j) < g^2/N < F_j H(j)$.\smallskip\\
{\bf Case 2:} Suppose $g/N \leq r_{j,k*}$.  Let $u(j) \leq (g-1)$ be the number of indices $k$ such that $0 < r_{j,k} < 1/N$. Definitions \eqref{defsjk*} and \eqref{defsjk} imply $s_{j,k*}\leq r_{j,k*} - (g-1)/N$ and $\sum_{k \neq k*}  s_{j,k} \leq u(j)/N + \sum_{k \neq k*}  r_{j,k}.$  This yields $S(j) \leq \sum_k  r_{j,k} < Fj H(j)$ since $r \in K(H)$.
\end{proof}
We now build our stochastic mutations model to force 
$
P( R_n/N \in K(H) \; | \; H_n = H ) = P( R_n/N \in K_N(H) \; | \; H_n = H ) = 1
$
We also want our stochastic mutations model to ensure that for large $N$ and given $H_n$, the distributions of $R_n(j,k)$ are approximately independent Poisson distributions with means $M_{j,k} N F_j H_n(j)$ proportional to the size $N F_j H_n(j)$ of the $j$-cell colony at the end of Phase 1. To this end, we introduce companion matrices $Z_n$.
\subsubsection{The Poissonian companion matrices $Z_n$} 
\label{compMat}
Given $H_n = H$, let $Z_n$ be a matrix of independent random variables $Z_n(j,k)$ having Poisson distributions with respective means 
$
\mbb{E}[ Z_n(j,k) \; | \; H_n= H ] = M_{j,k} N F_j H(j) = m Q_{j,k} N F_j H(j).
$ 
The probability $P( Z_n/N \in K(H) \; | \; H_n = H ) $ is \emph{smaller than 1}, but we will show that this probability tends very fast to 1 for $N$ large (see Theorem \ref{companionTH}). Thus, for all histograms $H$ and matrices $z$ with non-negative coefficients, define the conditional distribution of $R_n$ given $H_n=H$ by 
$
P( R_n/N = z \; | \; H_n = H )  = P(Z_n/N = z \; | \; Z_n/N \in K(H)). 
$
Since $Z_n/N$ is $N$-rational, this forces $R_n/N$ to be $N$-rational with
$$
P( R_n/N \in K_N(H) \; | \; H_n = H ) = P( R_n/N \in K(H) \; | \; H_n = H ) = 1.
$$
We can now complete the analysis of random mutations by describing the population histogram after Phase 2.

\subsubsection{Population histogram $J_n$ after mutations}
\label{JnMut}
After accouting for all the random mutations $R_n(j,k)$ at the end of Phase 2, each $j$-cell colony has lost $\sum_k R_n(j,k)$ outgoing mutants and gained $\sum_k R_n(k,j)$ incoming mutants, but the total population size is still $N_{ter} = \ceil{N \langle F, H_n \rangle}$. Therefore, the number $J_n(j)$ of $j$-cells at the end of Phase 2 is given by
\begin{equation} \label{Wn}
J_n(j) = \frac{1}{\ceil{ N \langle F, H_n \rangle}}  \left(\ceil{N F_j H_n(j)} - \sum_k R_n(j,k) + \sum_k R_n(k,j)\right)
\end{equation}
where we note that $R_n(j,k) = 0$ for all $j \notin spt(H)$.

Hence, given $H_n= H \in \mathcal{H}_N$ and $R_n/N = r \in K_N(H)$, the population histogram $J_n$ at the end of Phase 2 is a deterministic function $J_n= J(H_n, R_n/N, N)$ of $H_n, R_n/N$ and $N$. The histogram-valued function $(H, r, N) \to J(H,r,N)$ does not depend on $n$ and is actually defined for all $H \in \mathcal{H}$ and all $r \in K(H)$ by the formula (derived from \eqref{Wn}), 
\begin{equation} \label{WWn}
J_j (H,r,N) =\frac{ \ceil{ N F_j H(j)} }{ \ceil{ N \langle F, H \rangle }} + \; \frac{N}{\ceil{ N \langle F, H \rangle } } \left[- \sum_k  r_{j,k} + \sum_k  r_{k,j} \right].
\end{equation}
One has $\left|\frac{\ceil{ N F_j H(j)}}{\ceil{ N \langle F, H \rangle }} - \frac{F_j H(j)}{\langle F, H \rangle} \right| \leq \frac{4 F_g}{N F_1}$ due to \eqref{intervalphi}. Consequently, given $H_n= H \in \mathcal{H}_N$ and $R_n/N = r \in K_N(H)$, we expect $J_n(j) $ to be well-approximated by
\begin{equation} \label{psi}
\Psi_j(H, r) = \frac{1}{\langle F, H \rangle} \left(F_j H(j) - \sum_k  r_{j,k} + \sum_k  r_{k,j}\right) .
\end{equation}
Then, all $\Psi_j(H, r)$ remain well-defined for all $H \in \mathcal{H}$ and $r \in K(H)$. Furthermore, all $\Psi_j(H,r)$ verify $\Psi_j(H, r) \geq 0$ and $\sum_j  \Psi_j(H, r) =1$. Hence, they determine a histogram $\Psi(H, r)$.

We now prove the above claim. Precise but elementary accuracy computations show that for $N > 20,$ one has 
\begin{equation} \label{JJn}
\vnorm{J(H, r, N) - \Psi(H, r)} \leq \frac{13 F_g}{N F_1}
\end{equation}
for all $H \in \mathcal{H}$ and  $r \in K(H).$ In particular, one has with probability 1
\begin{equation} \label{Jn}
\vnorm{J_n - \Psi(H_n , R_n/N)} \leq \frac{13 F_g}{N F_1}.
\end{equation}
Since $J_n$ is an $ N_{ter}$-rational histogram, its essential minimum verifies 
\begin{equation} \label{bJn}
b(J_n) \geq \frac{1}{N_{ter}} \geq \frac{1}{N \langle F, H_n \rangle} \geq \frac{1}{N F_g}. 
\end{equation}
For $H \in \mathcal{H}$ and $r \in K(H)$, equations \eqref{Kj} and \eqref{psi} yield for each $j$ 
\begin{equation} \label{zeroPsi}
\Psi_j(H, r) = 0 \quad \text{if and only if} \; \; H(j) = r_{j,k} = r_{k,j}= 0 \;\; \text{for all} \; k. \\
\end{equation}
The relations \eqref{Wn} and $R_n/N \in K_N(H_n)$ imply $J_n(j) = 0$ if and only if $H_n(j) = R_n(k,j) = 0$ for all $k.$ Similarly, \eqref{WWn} shows for all $N$-rational $H$ and $r \in K(H)$
that $J_j(H, r, N ) = 0$ if and only if $H(j) = r_{k,j} = 0$ for all $k.$ Hence, one has $spt(J_n) = spt( \Psi(H_n , R_n/N) )$ for all $H_n$ and $R_n/N \in K(H_n).$  In addition, $spt(J(H, r, N))  =  spt( \Psi(H, r) )$ for all $N$-rational $H$ and $r \in K(H)$.
Equation \eqref{zeroPsi} proves that for $H, H' \in \mathcal{H}$, $r \in K(H),$ and $r' \in K(H')$,
\begin{equation} \label{sptpsiHh}
spt(\Psi(H,r) ) = spt(\Psi(H', r') \; \; \text{whenever} \; spt(H) = spt(H') \; \; \text{and} \; spt(r) = spt(r').
\end{equation} 
For all $H$ and $r \in K(H)$, the partial derivatives of $\Psi(H,r)$ are given by 
\begin{equation} \label{difPsi}
\begin{cases}
\frac{\partial}{\partial H(i)} \Psi_j  =  \frac{- F_i F_j H(j)}{\langle F, H \rangle^2} + \frac{1_{\{ i=j \}} F_j}{\langle F, H\rangle}, \\
\frac{\partial}{\partial r_{j,k}} \Psi_j  =  \frac{- F_j H(j)}{\langle F, H \rangle} \\
 \frac{\partial}{\partial r_{k,j}} \Psi_j  =  \frac{F_j H(j)}{\langle F, H \rangle}.
\end{cases}
\end{equation}
where $1_{\{i=j\}}$ denotes the indicator function that is $1$ when $i=j$ and $0$ otherwise.  Since $F_j H(j) \leq  \langle F, H \rangle$ and $\langle F, H\rangle  \geq  F_1$, we have that
\begin{equation} \label{lipPsi}
\vnorm{\Psi(H',r') - \Psi(H,r)} \leq \frac{3 g F_g}{F_1} \left(\vnorm{ H' - H } + \vnorm{ r' - r }  \right)
\end{equation}
for all $H', H \in \mathcal{H}$, $r \in K(H)$, and $r' \in K(H').$

\subsection{Phase 3 - Random selection:} 
At the end of Phase 2 on day $n$, the population has size $N_{ter} \simeq N \langle F,H \rangle$ and histogram $J_n \simeq \Psi(H_n , R_n/N)$ with $\Psi$ given by \eqref{psi}. From this population, Phase 3 extracts a random sample of size $N$, which becomes the new initial population $pop_{n+1}$ on day $n+1$ and has genetic histogram $H_{n+1}$. Phase 3 is thus a very simplified emulation of natural selection. 

The multinomial distribution $\mu_{N,J}(V)$ parameterized by $N$ and any histogram $J$ is defined for all $V \in \mbb{R}^g$ with integer coordinates $V(j) \geq 0$ and $\sum_j V(j) = N$ by 
\begin{equation} \label{multinomial}
\mu_{N, J}(V) = N! \; \prod_{j \in spt(G)} \; J(j)^{V(j)} / V(j)!. 
\end{equation} 
When $J_n(j) =0$, no mutant of genotype $j$ is present before selection so that $H_{n+1}(j) = 0$. Hence $spt(H_{n+1}) \subset spt(J_n)$ with probability 1.
Therefore, for $G \in \mathcal{H}_N$, all coordinates of $V = N G$ are integers, and one has
\begin{equation}\label{multinomial2}
P( H_{n+1} = G \; | \; H_n ; R_n ) = P(H_{n+1} = G \; | \; J_n) = 
\begin{cases}
\mu_{N, J_n}(N G),&  spt(G) \subset spt(J_n), \\
 0 & \text{otherwise.}
\end{cases}
\end{equation}
The multinomial distribution $\mu_{N, J}$ has mean $N J$. Hence, equation \eqref{multinomial2} gives 
$
\mbb{E}[ H_{n+1} \, | \, H_n , R_n ]  =  \mbb{E}[ H_{n+1} \, | \, J_n] = J_n. 
$
\subsection{Markov chain dynamics in the space of histograms} \label{TransProb}
At the completion of the three phases during a daily cycle, the population size returns to the large but \emph{fixed size $N$}. The cycle on day $n$ thus induces a stochastic transition $H_n \to H_{n+1} $ in the space of genetic histograms $\mathcal{H}\subset \mbb{R}^g$. The succession of daily cycles comprised of the three phases described above generates a \emph{time-homogeneous Markov chain} $n \to H_n$ on the state space $\mathcal{H} $ of all histograms. However, for each fixed population size $N$, the actual state space of this Markov chain is the finite set $\mathcal{H}_N$ of $N$-rational histograms where $card(\mathcal{H}_N) \leq (N+1)^{g}$.

The \emph{transition kernel} $\mathcal{Q}_N (H, G) = P(H_{n+1} = G \; | \; H_n = H)$ for $H, G \in \mathcal{H}_N$ is given by the finite sum 
\begin{equation}\label{kernel}
\mathcal{Q}_N (H, G) = \sum_{r \; \in \; K_N(H) } \; P( R_n/N = r \; | \; H_n = H ) \; P( H_{n+1}= G \; | \; R_n/N = r, H_n=H ).
\end{equation}
Recall that $card(K_N(H)) \leq [ (1+F_g)(N+1) ]^{g^2}$. 

This Markov chain generates stochastic histograms trajectories $ \mbold{H} = [H_1\, H_2 \ldots H_T]$ of arbitrary duration $T$ belonging to a \emph{path space} $\Omega_T$ studied further on (see Section \ref{pathspace}), which will heavily depend on the various parameters introduced in this section.  We will only consider realistic parameter sets, though we do note that most results presented in later sections hold generically.
\subsection{The set $\mathcal{P}$ of process parameters} \label{PAR}
The Markov chain of population histograms $H_n$ is defined by the very small scale factor $m$ affecting all mutation rates and a fixed finite parameter set $\mcal{P}$, namely,
\begin{itemize}
\item the number $g \geq 2$ of genotypes, 
\item the multiplicative daily growth factors $F_j > 1$ for genotype $j$ where the $F_j$ are \emph{distinct} and \emph{ordered} by increasing fitness from $F_1$ to $F_g$,
\item the $g \times g$ \emph{transfer matrix} $Q_{j,k}$ which  after scaling by $m$ defines the matrix of mutation rates $M = m \; Q$ 
\end{itemize}
Recall that $m_j = m \left(\sum_k Q_{j,k}\right)$. Our theoretical large deviations results hold uniformly when $\log(1/m)$ and all $\log(1/m_j)$ are larger than $1 + \delta$ for some $\delta >0$. Our numerical implementations use realistic parameter values with $m$ and all $m_j$ inferior to $10^{-6}$ so that $\delta \geq 12.8$.
 
During the random evolution of the genetic histograms $H_n$, the size $N$ of populations $pop_n$ at beginning of day $n$ remains fixed for all times $n$ due to daily random selection. Since realistic bacterial evolutions involve very large values of $N$ ranging from $10^5$ to $10^8$, our study naturally focuses on accurate asymptotics of random genetic evolution as $N \to \infty$.  The initial genetic histogram $H_1$ and the population size $N$ will never belong to the parameter set $\mcal{P}$.










\section{Large Deviations for Daily Cycles}
\label{LDcycles}
Naturally, the distribution of histogram trajectories is composed of sequences of daily transitions whose distributions are given by \eqref{kernel}.  Consequently, a large deviations analysis of random histogram trajectories necessarily starts with a large deviations analysis of daily transitions.  These daily transitions involve two probabilistic steps: random mutations and selection.  Therefore, we will carry out a full large deviations analysis below for random mutations and selection, culminating in an explicit large deviations framework for daily transitions.  We note that since mutations approximately follow a Poisson distribution and random selection follows a multinomial distribution, we will make repeated use of Stirling's Formula.  In addition, we will seek to relax the $N$-rationality conditions for our histograms as done in the previous section. 

\begin{mydef}
\label{balls}
For any histogram $H$ define the ball $V_N(H)$ and the $N$-rational ball $B_N(H) \subset V_N(H)$ by $V_N(H)  =   \vset{ H' \in \mathcal{H} \;: \; \vnorm{H' - H} \leq \frac{2}{3 N}} $ and
$B_N(H)  =   \mcal{H}_N \cap V_N(H).$
\end{mydef} 
The following lemma gives a sense of the density of the $N$-rational balls in the space of histograms $\mcal{H},$ which will greatly aid in relaxing $N$-rationality.
\begin{lemma} \label{basicBN}
For all histograms $H$, one has $card(B_N(H)) \leq 2^g.$  Furthermore, if $N > \frac{2}{b(H)}$, one has $spt(H') = spt(H)$ and  $b(H') > b(H)/2$ for all $H' \in B_N(H).$
\end{lemma}
\begin{proof} 
The bound on $card(B_N(H))$ is obvious. Take $H' \in B_N(H)$. For $k \in spt(H')$, one has $H'(k) \geq 1/N$ since $H'$ is $N$-rational. This forces $H(k) \geq H'(k) - \frac{2}{3 N} >\frac{1}{3 N} >0 $. Hence $k \in spt(H)$ and $spt(H') \subset spt(H)$.
Since $\vnorm{ H' - H} \leq \frac{2}{3 N}$, one has for $N > \frac{2}{b(H)}$ and any $j \in spt(H)$, 
$$
H'(j) \geq H(j) - \frac{2}{3 N} \geq b(H) - \frac{2}{3 N} > 2 b(H)/3
$$
Thus $j \in spt(H')$ and $ spt(H) \subset spt(H')$ so that $ spt(H) =spt(H')$. The last inequality then yields $b(H') \geq 2 b(H)/3$, which proves the lemma.
\end{proof}
\subsection{Accuracy of the Poissonian mutation approximations}
We will explicitly quantify the accuracy of the Poissonian approximations made in Section \ref{compMat}. The well-known proof of Stirling approximation for factorials can be easily modified to provide the following \emph{universal} inequality verified by all integers $N \geq 1$
\begin{equation}\label{stirling}
\abs{\log(N!) - N \log(N/e)} \leq 2 \log(N). 
\end{equation}
We will need a few classic large deviations lemma and a definition.
\begin{lemma}[Poisson Large Deviations] 
\label{poissonLD}
Fix any $u >0,$ and let $X$ be a random variable having a Poisson distribution with mean $N u$. For any $v > u > w > 0$ and for any integer $N,$ the following inequalities hold:
\begin{align} 
P( X \geq N v )  &\leq  \exp (- N \left[ u + v \log (v / u) - v) \right]) \label{cramerN1}\\
P( X \leq N w )  &\leq  \exp (- N \left[ u + w \log (w / u) - w) \right]) \label{cramerN2}
\end{align}
Moreover, for any $v \geq 0$ such that $N v$ is an integer, one has 
\begin{equation}\label{XNv}
\frac{1}{N}\log{P( X= N v )} = - (u + v \log(v/u) - v) + o(N) 
\end{equation}
with $| o(N) | \leq 2 \log(N)/N $.
\end{lemma}
\begin{proof}
Write $X = \sum_{s=1}^N  X_s$ where $X_1, \ldots, X_N$ are independent random variables having identical Poisson distribution $\theta$ with mean $u$. The large deviations rate $\lambda$ of $\theta$, also called Cramer transform of $\theta$, is classically defined (see \cite{azenLDT2012} \cite{dembo}) as the Legendre dual function of the Laplace transform of $\theta$, which is given by the well known formula $\lambda(v) = u + v \log (v / u) - v$ for all $v \geq 0.$  The empirical mean $X/N$ must verify the large deviations inequality
$
P ( X \geq N v) \leq \exp(- N \lambda(v))
$
for all $N$ and $v > u$ (see \cite{azenLDT2012,dembo}).  This proves equation \eqref{cramerN1}.  An analogous line of reasoning proves equation \eqref{cramerN2}.

For any $v \geq 0$ such that $N v$ is an integer, one has $P(X= N v) = e^{- N u} (N u)^{N v} / (N v) !$ so that
$$
\frac{1}{N}\log{P(X= N v)} = - u + v \log(N u) - \frac{1}{N} \log ((N v) !).
$$
For $v > 0$, the uniform Stirling's formula \eqref{stirling} gives
$$
\frac{1}{N} \log ((N v)!)= v \log(N) + v \log(v) -v + o(N)
$$
with $\abs{o(N)} \leq 2\log{(N)} / N.$ The last two equations prove \eqref{XNv} for $v >0$. Finally, \eqref{XNv} is trivially true for $v =0$ with the convention $0 \log(0) = 0$. This proves the lemma.
\end{proof}
\begin{lemma} \label{sumexp} 
Fix a positive sequence $\epsilon(N)$ such that $\epsilon(N) \to 0$ as $N \to \infty$. Fix any set $E$ and any ``rate function'' $\lambda(s) \geq 0$ defined for all $s \in E$. Fix $c > 0,\beta > 0,$ and let $E_N$ be a finite subset of $E$ with $card(E_N) \leq c N^\beta$. Consider fast-vanishing exponentials $p_N(s) > 0$ indexed by $s \in E_N$ such that 
$
\frac{1}{N} \log{p_N(s)} = - \lambda(s) + o_s(N)
$ 
where $\abs{o_s(N)} \leq \varepsilon(N).$  Define $\Lambda(E_N) = \inf_{s \in E_N} \lambda(s) $.  The sum $p(E_N) = \sum_{ s \in E_N } p_N(s) $ satisfies
$$
\frac{1}{N}\log {p(E_N)} = - \Lambda(E_N) + o(N) 
$$
for all $N$ with $|o(N)| \leq \epsilon(N) + \beta \log(N) / N + \log(c) /N.$ 
\end{lemma}
\begin{proof}
Select $s(N) \in E$ such that $\Lambda(E_N) = \lambda(s(N))$. This yields the lower bound 
\begin{equation} \label{low}
\frac{1}{N}\log {p(E_N)} \geq \frac{1}{N}\log {p(s(N)} = - \lambda(s(N)) + o_{s(N)}(N) \geq \; - \Lambda(E_N) - \epsilon(N).
\end{equation}
We have $ p_N(s) \leq \exp (  - N [ \Lambda(E_N) + \epsilon(N) ]  )$ for all $s \in E_N$ by definition of $\Lambda(E_N) $.  This gives $ p(E_N) \leq \exp (  - N [ \Lambda(E_N) + \epsilon(N) ] ),$ which implies
$$
\frac{1}{N}\log {p(E_N)} \leq - \Lambda(E_N) + \epsilon(N) + \log(card(E_N)) / N
$$
Combining this upper bound with \eqref{low} concludes the proof.
\end{proof}
\begin{mydef} 
\label{d(H)}
We define the \emph{decay coefficient} $d(H)$ of any histogram $H$ by $d(H) = \exp( - \delta F_1  b(H) )$ where $\delta = \log(1/m) -1 > 0.$
\end{mydef}
Recall that $b(H) >0$ is the essential minimum of $H$ given by \eqref{spt(H)b(H)}, which implies $0 < d(H) <1 $. Moreover $N b(H) \geq 1$ for $H \in \mathcal{H}_N$ so that $d(H)^N \leq e^{- \delta \; F_1} < 1$ for all $H \in \mcal{H}_N.$ For $m \leq 10^{-6},$ one has $ \delta > 12.8.$ In actual laboratory experiments \cite{len94a,azencoop,Gordo2012,Desai2007} , one has $F_1 >100$ so that $d(H)^N < e^{-12800} < 10^{- 5500}$ is practically zero.
Given the histogram $H_n = H$ on day $n$, recall that $Z_n = [ Z_n(j,k) ]$ is a matrix of independent Poissonian random variables with respective means $\mbb{E}[Z_n(j,k)]=  N m Q_{j,k} F_j H(j) = N M_{j,k} F_j H(j)$. We now prove that for $N$ large, the random mutations matrix $R_n$ and its Poissonian companion $Z_n$ have nearly identical distributions conditioned on $H_n = H$.
\begin{theorem} \label{companionTH}
Let $R_n$ be the random mutations matrix with Poissonian companion matrix denoted $Z_n$. Fix any histogram $H^*$ with decay coefficient $d(H^*)$ given in Definition \ref{d(H)}. Let $c = 2 +\log(g) / F_1 $. 
Then, provided $N > \frac{c}{b(H^*)} $, one has for all $n$, $H \in B_N(H^*)$, and matrices $z \in K(H),$
\begin{equation} \label{control}
1 \leq  \frac{P( R_n/N = z \; | \; H_n = H )} {P( Z_n/N = z \; | \; H_n = H )} \leq 1 + 2 g\,  d(H^*)^{N/2}.
\end{equation}
Denote $f_N(z, H) =\frac{P( R_n/N = z \; | \; H_n = H )}{P( Z_n/N = z \; | \; H_n = H )},$ the conditional density of $R_n/N$ with respect to $Z_n/N$.  Fix any $0< a < 1$. Then, $f_N(z, H) \to 1$ at exponential speed as $N \to \infty$ where convergence is uniform for $b(H^*) \geq a$, $ H \in B_N(H^*)$, and $z \in K(H)$.  Therefore, the conditional joint distribution of the mutations matrix $[R_n(j,k)]$ given $H_n \in B_N(H^*)$ becomes, for large $N,$ practically equal to the product of the Poisson distributions with respective means $mQ_{j,k} N F_j H_n(j) = M_{j,k} N F_j H_n(j)$.
\end{theorem}

\begin{proof}
Consider any $N$-rational histogram $H$. Given $H_n = H$, let $Z_n$ be a matrix of independent Poissonian random variables $Z_n(j,k)$ with means $M_{j,k}N F_j H(j)$. The sums $S(j) = \sum_{k= 1}^g  Z_n(j,k) $ then have Poisson distributions with respective means $\mbb{E}[S(j)] = N m_j F_j H(j)$.  First, we will first show that for each $j$
\begin{equation} \label{controlCj}
P (Z_n \in N  K(j, H) \; | \; H_n = H ) \geq 1- d(H)^N 
\end{equation}
holds for all $N$ and $H$.  To do so, we consider the cases when $m_j H(j) > 0$ and $m_j H(j) = 0$ separately.\smallskip\\
{\bf Case 1:} Suppose $m_j H(j) > 0.$  For any $s > m_j$, apply \eqref{cramerN1} to $X = S(j)$ with $v = s F_j H(j)$ and $u = m_j F_j H(j)$ to obtain
\begin{equation} \label{tail.j}
P( S(j) \geq s N F_j H(j) \; | \; H_n = H) \leq \exp \left( - s N F_j H(j)[\log(1/m_j) -1] \right).
\end{equation}
Since $ \log(1/m_j) \geq \delta + 1,$ \eqref{tail.j} implies 
\begin{equation} \label {tail.j2}
P( S(j) \geq s N F_j H(j) \; | \; H_n = H) \leq \exp \left(- \delta s  N  F_j H(j) \right) \leq d(H)^{s N}
\end{equation}
for all $N, H$ and $s > m_j $.  By definition \eqref{Kj} of $K(j, H),$ 
$$P (Z_n / N \in K(j, H) \; | \; H_n = H ) = 1 - P(S(j) \geq N F_j H(j) \; | \; H_n = H).$$
Equation \eqref{tail.j2} with $s=1$ implies \eqref{controlCj} for all $N$ and $H$. \smallskip\\
{\bf Case 2:} Suppose $m_j H(j) = 0.$ Both situations $m_j = 0$ or $H(j) = 0$ imply $P(Z_n(j,k)=0 \:| \:H_n= H ) = 1$ for all $k.$ Therefore, $S(j) = 0$ so that $P ( Z_n \in N K(j, H_n) \; | \; H_n = H ) =1$, which trivially satisfies \eqref{controlCj}
\smallskip\\
Now, since $K(H) = \bigcap_{j = 1}^g  K(j, H)$ and \eqref{controlCj} holds for all $N$ and $H,$
\begin{equation} \label{control2}
P (Z_n / N \in K(H_n) \; | \; H_n = H) \geq 1- g \,d(H)^N
\end{equation}
also holds for all $N$ and $H$. The constraint $ N > N(H) = \log(2 g) / (  12  F_1 b(H) ) $ forces $g \,d(H)^N < 1 / 2$ so that $1/(1 - g \, d(H)^N) \leq 1+2 g \, d(H)^N $. So for $N > N(H)$, equation \eqref{control2} implies
\begin{equation} \label{controlC}
1 \leq \frac{1}{P( Z_n / N \in K(H) \; | \; H_n = H)} \leq 1 + 2 g \, d(H_n)^N.
\end{equation}
For all matrices $z \in K(H)$, we have
$$
P (R_n/N = z \; | \; H_n = H) = \frac{P ( Z_n / N = z \; | \; H_n =H )} {P ( Z_n/N \in  K (H) \; | \; H_n = H)}.
$$
For $H_n= H$, $z \in K(H)$, and $N > N(H)$, equation \eqref{controlC} yields
\begin{equation} \label{compar}
P ( Z_n/N = z \; | \; H_n =H) \leq P ( R_n/N = z \; | \; H_n = H) \leq \left(1+ 2 g \, d(H)^N\right)  P ( Z_n/N = z \; | \; H_n = H ). 
\end{equation}
Fix any $H^* \in \mathcal{H}$. Impose $H_n = H \in V_N(H^*)$, which is equivalent to $H \in B_N(H^*)$ since $H_n$ is $N$-rational. By Lemma \ref{basicBN}, for $N > \frac{2}{b(H^*)}$, we must have $spt(H) = spt(H^*) $ and $b(H) \geq b(H^*)/2$ so that $d(H)< d(H^*)^{1/2}$.  Let $N_1 = \frac{c}{b(H^*)}$ with $c = 2 + \frac{\log(g)}{F_1}$. Then $N > N_1$ will force $N > N(H)$ provided $H \in B_N(H^*)$. Equation \eqref{compar} then implies, for all $N$-rational $H \in B_N(H^*)$, $z \in K_N (H)$, and $N > N_1$,
\begin{equation} \label{N3}
1 \leq  \frac{P( R_n/N = z \; | \; H_n = H )} {P( Z_n/N = z \; | \; H_n = H )} \leq 1 + 2 g \, d(H^*)^{N/2}, 
\end{equation}
which proves \eqref{control}. The uniformity in the statement of the theorem is an easy consequence of \eqref{control}.
\end{proof}
\subsection{Large deviations for random mutations}
With the Poissonian approximation complete and explicitly verified, we can now begin with the large deviations analysis of random mutations, which will yield a \emph{rate function} for random mutations.  Key in this analysis will be the topology of the mutation matrices under consideration and the regularity of this rate function. 

We begin by endowing the set of $g \times g$ matrices $r$ with the $L_{\infty}$ norm $\vnorm{ r } = \max_{j,k} \abs{ r_{j,k} } $. Most results will hinge on the feasibility of mutations and will also explicitly depend on the smallest nonzero mutation rate, which motivates the following definition.
\begin{mydef}
For any $g \times g$ matrix $r$ with all $r_{j,k} \geq 0$ and all $r_{j,j} = 0$, define the \emph{support} $spt(r)$ and the \emph{essential minimum $b(r)$ of $r$,} respectively, by
$$
spt(r) = \{ (j,k) \; | \; r_{j,k} >0 \} \qquad  \qquad b(r) = \min_{(j,k) \in spt(r)} \; r_{j,k}. 
$$
\end{mydef}
Given $H_n = H$, we have seen that the mutation matrix $ R_n / N$ belongs almost surely to the convex set $ K_N(H) \subset K(H)$ given in Definition \ref{KNconstr}. For any matrix $r$ in $K(H)$ and pairs of indices $(j,k) \in spt(r),$ one must have $H(j) >0$ and $M_{j,k} = m Q_{j,k} > 0$.   With this in mind, we are now ready to define the rate function for random mutations.  For $r \in K(H),$ let $V_N(r)$ be the $L_{\infty}$ ball of center $r$ and radius $\frac{2}{3 N}$. Let $B_N(r) \subset V_N(r)$ be the set of $N$-rational matrices in $V_N(r)$. Note that $card( B_N(r) ) \leq 2^{g^2}$.
\begin{mydef}[Rate function for mutations]
 For any histogram $H\in \mcal{H}$ and $r \in K(H),$ define the $g \times g$ matrix $L(r,H)$ of Poissonian rate functions by
\begin{equation} \label{Ljk}
L_{j,k}(r,H)  = 
\begin{cases}
M_{j,k} F_j H(j) + r_{j,k} \log\left( \frac{r_{j,k}} { e M_{j,k} F_j H(j) }\right), & H(j) M(j,k) >0; \\
 0,& H(j) M(j,k) = 0.
\end{cases} 
\end{equation}
The large deviations \emph{rate function for mutations} $mut(r,H)$ is given by
\begin{equation}\label{mutcramer}
mut(r,H) = \sum_{j,k} \; L_{j,k}(r,H) = \sum_{(j,k) \; | \; M_{j,k} H(j) > 0} 
M_{j,k} F_j H(j) + r_{j,k} \log r_{j,k} - r_{j,k} \log(e M_{j,k} F_j H(j) ) 
\end{equation}
with the convention $0 \log(0) = 0.$
\end{mydef}
For $r \in K(H)$ and $H(j) M_{j,k} >0,$ the function $ L_{j,k} (r,H) \geq 0$ is the classical rate function from large deviations theory of the Poisson distribution with mean $M_{j,k} F_j H(j)$. Hence, $ L_{j,k} $ is a \emph{strictly convex} function of $r_{j,k}$ for $r_{j,k} >0$, and one has $L_{j,k}(r,H)=0$ if and only if $r_{j,k}= M_{j,k} F_j H(j)$.  Consequently, $mut(r,H)$ is a finite, non-negative, \emph{continuous, strictly convex} function of $r \in K(H)$ since each $L_{j,k}$ is strictly convex in $r_{j,k}$. Moreover, $mut(r,H)$ is a continuous, convex function of $H$ for fixed $r$.
Note the basic inequalities
\begin{align*} 
 x \abs{ \log x } \leq \mcal{D}(v) &:= 1/e + \log(1+v), \qquad 0 \leq x \leq v; \\
 \abs{ \log x } \leq \mcal{E}(u,v) &:= \log(1+1/u) + \log(1+v), \qquad  0 < u \leq x \leq v.
\end{align*}
For $b(H) >a$, $r \in K(H),$ and $M_{j,k} H(j) > 0$, one has $r_{j,k} \leq F_j H(j) \leq F_g $ and 
$$
0 < a \, b(M) \leq M_{j,k} F_j H(j)) \leq F_g \vnorm{M}.
$$ 
Then, we get $\abs{ r_{j,k} \log(r_{j,k}) }\leq\mcal{D}(F_g)$ and 
$$
\abs{ r_{j,k} \log(e M_{j,k} F_j H(j)) }\leq F_g \abs{\log(e M_{j,k}) }+ F_j H(j) \abs{ \log(F_j H(j)) } \leq \mcal{E}\left(a \, b(M) , F_g || M || \right) + \mcal{D}(F_g). 
$$
Hence, for all $H$ with $b(H) \geq a >0$ and $r \in K(H),$ we have the uniform bound
\begin{equation} \label{mutbound}
mut(r,H) \leq g^2 \left[F_g \vnorm{M} + \mcal{E}\left(a\, b(M), F_g\vnorm{ M }\right) + 2 \mcal{D}(F_g)\right] = c
\end{equation} 
where the constant $c$ depends only on $a$ and the parameter set $\mcal{P}$.  Note also that for $r \in K(H)$, 
\begin{equation}\label{zeromutcramer}
mut(r,H) = 0 \iff L_{j,k}(r,H) = 0  \iff r_{j,k} = m Q_{j,k} F_j H(j).
\end{equation}
Our continuity and asymptotic results presented below will be uniform over a certain compact set of histograms, defined below.  These results will frequently make use of the above analysis as well.
\begin{mydef}
For a fixed constant $0 < a < 1,$ and define the compact set of histograms $\mcal{H}(a) \subset \mcal{H}$ by
\begin{equation} \label{HIST(a)}
\mcal{H}(a) = \vset{H \in \mathcal{H} \; | \; b(H) \geq a} 
\end{equation}
\end{mydef}
\begin{proposition} \label{propholdermut}
Fix $0 < a < 1$ and any H\"{o}lder exponent $0 < \alpha < 1$. Set $c = \frac{ 14 g^2 F_g \log(F_g) }{a (1- \alpha)}$. For all histograms $H', H \in \mcal{H}(a)$ with $spt(H') = spt(H)$ and matrices $r$ and $r'$ with $r \in K(H)$ and $r' \in K(H'),$ the mutations rate function verifies 
\begin{equation}\label{holdermut}
\abs{ mut(r', H') - mut(r, H)} \leq c \left( \vnorm{r' - r}^{\alpha} + \vnorm{ H'- H}\right).
\end{equation} 
\end{proposition}
\begin{proof}
Fix $0 < \alpha <1$. For all $x,y\in [0,1/2],$ we claim that 
\begin{equation} \label{xlogxholder}
\abs{y \log y - x \log x } \leq \frac{4}{1-\alpha} \abs{x - y}^{\alpha}.
\end{equation}
Indeed, for $0 \leq x \leq 1$, one has
\begin{equation}\label{xlogx}
| x \log(x) | = x \log (1/x) \leq (x/(1-\alpha)) \log (1/x^{1-\alpha}) \leq x^{\alpha} / (1-\alpha).
\end{equation}
Consider $\abs{y \log y - x \log x }$ with $ \abs{y-x} < w \leq 1$ . Then, either $x, y \leq 2 w$ or $x,y \geq w$. When $x, y \leq 2 w$, equation \eqref{xlogx} gives 
$$
\abs{y \log y - x \log x } \leq \abs{x \log x} + \abs{y \log y} \leq \frac{4}{(1-\alpha) w^{\alpha}}. 
$$
When both $x,y \geq w$, Taylor's formula and \eqref{xlogx} yield 
$$
\abs{y \log y - x \log x } \leq w ( 1 + \log(1/w) ) \leq \frac{4}{1-\alpha} w^{\alpha}. 
$$
This proves \eqref{xlogxholder}.

Continuing further, for $x,y \in [0,A]$ with $A >1$, apply \eqref{xlogxholder} to $x/A$ and $y/A$ to easily get 
\begin{equation}\label{basic.holder}
\abs{x \log x  - y \log y } \leq \left(\frac{5 (1 + A ) \abs{\log(A)} }{1-\alpha}\right) \abs{x - y}^\alpha .
\end{equation}
Take $H, H' \in \mcal{H}(a)$ with $spt(H') = spt(H)$, $r' \in K(H'),$ and $r \in K(H)$. Fix $(j,k)$ with $H(j) M_{j,k} > 0$, which implies $H'(j) M_{j,k}> 0$. Then, $ x := r'_{j,k} $ and $y := r_{j,k}$ verify $\abs{ x- y } \leq \vnorm{ r' - r }$ and are bounded by $\max\vset{ F_j H'(j),  F_j H(j)} \leq F_g $.  By definition, 
\begin{align*}
L_{j,k}(r,H) &= y \log y - y \log(e F_j H(j)) + F_j H(j), \\
L_{j,k}(r',H') &= x \log x - x \log( e F_j H'(j)) + F_j H'(j).
\end{align*}
Set $\hat{L} := \abs{ L_{j,k}(r, H) - L_{j,k}(r', H') },$ which satisfies the inequality 
\begin{equation} \label{dL}
\hat{L} \leq \abs{ y \log y  - x \log x } + \abs{y-x}  \abs{ \log(e F_j H(j)) } + \abs{x}  \abs{ \log ( H'(j)/H(j) ) } + F_j  \abs{ H'(j)- H(j) }. 
\end{equation}
We decompose \eqref{dL} by writing $\hat{L} \leq U_1 + U_2 + U_3 + U_4$ where the term $U_i$ corresponds to the $i$th term of the sum in \eqref{dL}.
Clearly $U_3 + U_4 \leq 2 (F_g / a) \vnorm{ H' - H }$.  Then, \eqref{basic.holder} gives $U_1 \leq \frac{c_1}{1-\alpha}\vnorm{ r' - r }^\alpha$ with $c_1 =5 (1+F_g) \log(F_g) $. Since $F_g \geq F_j H(j) \geq a F_1$, one has $\abs{\log(e F_j H(j) )} \leq c_2 /a$ with $c_2 = 2 + \log(F_g/F_1)$, and hence, $ U_2 \leq c_2 \vnorm{r'- r}$. These bounds yield for $H(j) M_{j,k} > 0$
$$
\hat{L} \leq  \frac{\hat{c}}{a (1-\alpha) }\left( \vnorm{ r' - r }^{\alpha} + \vnorm{ H' - H }\right) 
$$
with $\hat{c} = c_1+ c_2 + 2 F_g \leq 14 F_g \log(F_g) $. This result still holds when $H(j) M_{j,k} = 0$ since $\hat{L} = 0$.  Hence, the mutations rate function $mut(r,H)= \sum_{j,k} L_{j,k}(r,H)$ verifies for all $r \in K(H), r' \in K(H') $ and $H,H' \in \mcal{H}(a),$
$$
\abs{ mut(r', H') - mut(r, H) }\leq \frac{\hat{c} g^2 }{a (1-\alpha)}\left(  \vnorm{ r' - r }^{\alpha} + \vnorm{ H' - H}\right) \leq c\left(  \vnorm{ r' - r }^{\alpha} + \vnorm{ H' - H}\right),
$$
where $c$ is the constant stated in this proposition. This concludes the proof.
\end{proof}
We now conclude the large deviations analysis for random mutations with a key asymptotic result that will be used for the large deviations analysis of daily transitions later on in this section.
\begin{proposition} \label{LDmut}
Let $R_n$ be the random matrix of mutations on day $n$. Let $mut(r, H)$ be the mutations rate function defined by \eqref{mutcramer}. Fix $a >0$ and the parameters $\mcal{P}$. There is a constant $N_0 = N_0(a, \mcal{P})$ such that for $N > N_0$, the large deviations formula
\begin{equation} \label{LDmut1}
\frac{1}{N} \log P( R_n / N = r  \; | \;  H_n = H )  =  - mut(r,H) + o(N)
\end{equation}
with  $\abs{ o(N) }  \leq  4 g^2  \log N / N $ holds uniformly for all $H$ and $r$ where $H\in \mcal{H}(a) \cap \mcal{H}_N$ and $r \in K_N(H).$ 
\end{proposition}
\begin{proof}
Take $H$ and $r$ as stated above. Given $H_n = H $, the coefficients $Z_n(j,k)$ of the companion matrix $Z_n$ are independent and have Poisson distributions with respective means $N M_{j,k} F_j H(j)$. For $j \in spt(H)$ and any $k$, set $u = M_{j,k} F_j H(j)$ and $v = r_{j,k}$. Since $N v$ is an integer, apply \eqref{XNv} to $X=Z_n(j,k)$ to obtain 
\begin{equation}\label{L0}
\frac{1}{N}\log P(Z_n(j,k) / N  =  r_{j,k} \; | \; H_n = H) = - L_{j,k}(r,H) + o_1(N)
\end{equation}
 with $\abs{ o_1(N) } \leq 2 \log N /N.$ This equation remains true for $j \notin spt(H)$ and all $k$ since $H(j) = 0$ and $Z_n(j,k) = 0$.

Given $H_n = H$, the coefficients of $Z_n$ are independent so that
$$ 
P( Z_n/N = r \: | \; H_n = H)  =  \prod_{j , k} \; \; P ( Z_n(j,k)/N = r_{j,k} \; | \; H_n = H). 
$$
Equation \eqref{L0} implies
\begin{equation}\label{L1}
\frac{1}{N} \log P( Z_n / N = r \; | \; H_n = H )  =  - \sum_{j, k} L_{j,k}(r, H) + o_2(N)  =  - mut(r, H) + o_2(N)
\end{equation}
with $\abs{ o_2(N) }  \leq g^2 \abs{ o_1(N) } \leq 2 g^2  \log N /N$.

Set $c = 2 + \log(g)$, and impose $N >c/a \geq \frac{c}{b(H)}$. Applying \eqref{control} yields
\begin{equation} \label{L3}
\left|\frac{1}{N} \log P (  R_n/N = r \; | \; H_n = H ) - \frac{1}{N}\log P ( Z_n/N = r \; | \; H_n = H )\right| \leq  \frac{1}{N} \log\left( 1 + 2 g \, d(H)^{N/2} \right). 
\end{equation}
The right-hand side of \eqref{L3} is bounded above by $\frac{1}{N} \log(1 + 2 g)$ so that equation \eqref{L1} implies
\begin{equation}\label{L4}
\frac{1}{N} \log P ( R_n/N = r \; | \; H_n = H )  =  - mut(r, H) + o(N) 
\end{equation}
with $\abs{ o(N) }  \leq  \frac{1}{N} \left( \log(1 + 2 g) + 2 g^2 \log N  \right) \leq 4 g^2 \log(N)/N $.
\end{proof}
\subsection{Large deviations for random selection} \label{LDsel}
The final phase of random daily cycles is multinomial selection, so we now present the explicit large deviations analysis for random selection.  The classical rate function for a multinomial distribution is linked to the well-known Kullback-Leibler divergence.  Thus, we will detail the resulting rate function for multinomial sampling in this bacterial context and discuss regularity of the rate function, which will culminate in an asymptotic result for random selection.
\begin{mydef}[Rate Function for Random Selection]
The classical \emph{Kullback-Leibler divergence} between two histograms $G$ and $J$ is defined by 
\begin{equation} \label{KL}
KL(G,J)  = 
\begin{cases}
\sum_{j \in spt(G)}  G(j) \log \frac{G(j)}{J(j)}, & spt(G) \subset spt(J); \\
\infty, & \text{otherwise}. 
\end{cases}
\end{equation}
\end{mydef}
Recall that $KL(G,J)\geq 0$ for all $G$ and $J$, and $KL(G,J) = 0$ if and only if $G = J.$ 
For all pairs of histograms $(G,J)$ such that $ spt(G) \subset spt(J), $ the function $KL(G,J)$ is finite, continuous, and has partial derivatives with respect to $G(j)$ and $J(j)$ for all $j \in spt(G)$, which are given by
\begin{align}
\frac{\partial}{\partial G(j)}  KL(G,J)  &=  1 + \log G(j) - \log J(j) \label{partialKLG},\\
\frac{\partial}{\partial J(j)}  KL(G,J)  &=  - \frac{G(j)}{J(j)} \label{partialKLJ},\\
\frac{\partial^2}{\partial G(i) \partial G(j)}  KL(G,J)  &= \left(\frac{1_{\vset{ i  =  j }}}{G(j)}\right)
\frac{\partial^2}{\partial J(i)\partial J(j)} KL(G,J)  = 1_{\vset{ i  =  j }} \frac{G(j)}{J(j)^2}.
\end{align}
For $G$ fixed, $KL(G,J)$ is a strictly convex, differentiable function of $J$ on the convex set of all $J$ such that $spt(G) \subset spt(J)$. For $J$ fixed, $KL(G,J)$ is also a strictly convex, differentiable function of $G$ on the convex set of all $G$ such that $spt(G) \subset spt(J)$. We now evaluate uniform continuity moduli for $KL(G,J)$.
\begin{lemma} \label{lemma.boundsKL}
For any histograms $G$ and $J$ with $spt(G) \subset spt(J)$, or equivalently $KL(G, J)< \infty $, define
\begin{equation} \label{betadef}
\beta(G,J) = \max_{k \in spt(G)} { 1/J(k) }.
\end{equation}
One has then 
\begin{equation} \label{betabound}
KL(G,J) \leq \log \beta(G,J) \leq \frac{KL(G,J) + \log(g)}{b(G)}
\end{equation}
\end{lemma}
\begin{proof}
The entropy $\mscr{E}(G)$ of $G$ verifies $0 \leq \mscr{E}(G) = - \sum_{j \in spt(G)} G(j) \log(G(j)) \leq \log(g)$.  By definition of $KL(G,J)$, one has 
\begin{equation} \label{beta2}
KL(G, J) + \mscr{E}(G) = \sum_{j \in spt(G)}  G(j) \log(1/J(j)).
\end{equation}
For each $k \in spt(G)$, this yields $b(G) \log{(1/J(k))} \leq KL(G, J) + \mscr{E}(G)$.  By definition \eqref{betadef}, we thus get $\log{\beta(G,J)} \leq [ KL(G, J) + \log(g) ] / b(G)$.
Since $G$ is a histogram, equation \eqref{beta2} forces 
$$ 
KL(G,J) \leq KL(G, J) + \mscr{E}(G) \leq \max_{j \in spt(G)}  \log(1/J(j)) = \log(\beta(G,J)).
$$
\end{proof}

\begin{proposition} \label{lipKL}
\hfill
\begin{enumerate}[label = (\roman*)]
\item For $G, J, I\in \mathcal{H}$ with $spt(G) \subset spt(J) \cap spt(I)$, one has 
\begin{equation}\label{lipKLJI}
\abs{KL(G, J) - KL(G, I)}  \leq  g^{1 + 1/b(G)}  e^{\kappa / b(G)}  \vnorm{ J - I } 
\end{equation}
with $\kappa = \min \vset{ KL(G, J),  KL(G, I) }$.
\item Fix $0 < a \leq 1/2$. For $G, G', J \in \mathcal{H}$ verifying $ spt(G) \cup spt(G') \subset spt(J), b(G) \geq a ,$ and  $b(G') \geq a,$  we have
\begin{equation} \label{holderKLG}
\abs{KL(G',J) - KL(G,J)} \leq c \vnorm{ G' - G } ^{1/2}
\end{equation} 
with $c = \frac{g}{a} \left(2 + \log g  + \max \vset{ KL(G,J), KL(G',J)}\right)< \infty.$ 
\item For $G, G', J$ verifying $sp(G) = spt(G'), spt(G) \subset spt(J), b(G) \geq a,$ and $b(G') \geq a,$ one has the Lipschitz continuity 
\begin{equation} \label{lipKLG}
\abs{KL(G',J) - KL(G,J)} \leq c_1 \vnorm{ G' - G } 
\end{equation} 
with $c_1 = \frac{g}{a} \left(2 + \log g + \min\vset{KL(G,J) , KL(G',J) }\right)$.
\end{enumerate}
\end{proposition}
\begin{proof}
\hfill
\begin{enumerate}[label = (\roman*)]
\item Take $G, J, I \in \mcal{H}$ as stated. By symmetry in $J$ and $I$, we may assume $KL(G, J) \leq KL(G, I)$.  Formula \eqref{partialKLJ} implies $\abs{\frac{\partial}{\partial J(k)} KL(G, J)}\leq \beta(G,J)$ for all $k \in spt(G)$ where $\beta(G,J)$ is given by \eqref{betadef}. By Taylor's formula and \eqref{betabound}, we get 
\begin{equation}
\abs{KL(G,J) - KL(G,I)}  \leq  g  \beta(G,J) \vnorm{ J - I } \leq c  \vnorm{ J - I }
\end{equation}
with $c = g  \exp ( [ K(G,J) + \log(g) ] / b(G) )$, proving item (i).
\item Now take $a, G, G',$ and $J$ as stated. For $j \in spt(G) \cup spt(G')$, let $u(j) = [G'(j) - G(j) ]\log J(j)$ and $v(j) = G'(j) \log G'(j) - G(j) \log G(j).$ By definition of the Kullback-Leibler divergence, we have
\begin{equation} \label{sumKL}
 \abs{KL(G',J) - KL(G,J)} \leq \sum_{ j \in spt(G) \cup spt(G') } ( \abs{u(j)} + \abs{v(j)}).
\end{equation} 
Equation \eqref{xlogxholder} with $\alpha = 1/2$ implies
\begin{equation} \label{generic.v(j)}
| v(j) | \leq \; 8 \; ||\; G' - G \;|| ^{1/2}, \qquad j \in spt(G) \cup spt(G'). 
\end{equation}
By definition \eqref{betadef} and the bound \eqref{betabound}, one has for $ j \in spt(G) \cup spt(G') \subset spt(J)$, 
\begin{align}
\abs{ u(j) } &\leq \vnorm{ G'- G }  \max\vset{\log  \beta(G,J)  , \log \beta(G,J)}\nonumber \\ 
 &\leq \vnorm{ G'- G }\left( \frac{1}{a}[ \log g + \max\vset{KL(G,J), KL(G',J)}]\right).  \label{generic.u(j)}
\end{align}
Combine \eqref{generic.u(j)}, \eqref{generic.v(j)}, \eqref{sumKL} to get 
\begin{equation}
\abs{KL(G',J) - KL(G,J)} \leq c_1 \vnorm{ G' - G } ^{1/2} 
\end{equation} 
with $c_1 = \frac{g}{a} \left(8 + \log g + \max\vset{ KL(G,J) , KL(G',J)}\right),$ proving item (ii).
\item Finally, take $G, G',$ and $J$ as stated in (iii). From \eqref{betadef} and \eqref{betabound}, one gets for $j \in spt(G) = spt(G') $
\begin{equation*} 
\abs{ u(j) } \leq \frac{1}{a} \vnorm{ G'- G }\left(\log g + \min \vset{KL(G,J), KL(G',J)}\right).
\end{equation*}
For $j \in spt(G) = spt(G'),$ Taylor's formula implies, 
\begin{equation} 
\abs{v(j)} \leq \vnorm{ G' - G} \left[ 1 + \log (1/a)\right] 
\end{equation}
These last two bounds along with \eqref{sumKL} prove (iii).
\end{enumerate}
\end{proof}
We now state the key asymptotic results for multinomial sampling, which coupled with Proposition \ref{LDmut}, will contribute to an asymptotic result for daily transitions.
\begin{proposition} \label{LDmultinomial}
Consider $J, G\in \mcal{H}$ with $spt(G) \subset spt(J)$. When $G$ is $N$-rational, the multinomial distribution $\mu_{N,J}$ defined by \eqref{multinomial} verifies
\begin{equation} \label{KLlogmu}
\frac{1}{N} \log(\mu_{N,J}(N G) )= - KL(G,J) + o(N) 
\end{equation}
with \emph{uniform} remainder $\abs{ o(N) } \leq 2 (g+1) \log N /N. $
\end{proposition}
\begin{proof}
The coordinates of $V= N G$ are non-negative integers with sum $N$ so that \eqref{multinomial2} gives 
\begin{equation}\label{logmu}
\frac{1}{N} \log(\mu_{N,J}(N G)) = \frac{1}{N} \log N! - \sum_{j \in spt(G)} \frac{1}{N} \log V(j)! + \sum_{j \in spt(G)}  \frac{V(j)}{N} \log J(j). 
\end{equation}
For $j \in spt(G) \subset spt(J)$, apply Stirling's formula \eqref{stirling} to $V(j)!= [ N G(j)]!$ to get, 
\begin{equation}\label{logV}
\frac{1}{N}\log V(j)! = G(j)\log N + G(j) \log G(j) - G(j) + o_j(N)
\end{equation}
with $\abs{ o_j(N) } \leq 2\abs{\log(N G(j)) / N} $.  Since $G \in \mathcal{H}_N$, the integer $V(j) = N G(j)$ is positive for all $j \in spt(G)$ so that $1 \leq N G(j) \leq N$.  Therefore, we have $0 \leq \log( N G(j) ) \leq \log(N)$ and $\abs{o_j(N)}  \leq  2  \log N /N .$  Using \eqref{stirling} for $\log(N!),$ along with\eqref{logV}, for $spt(G) \subset spt(J)$, equation \eqref{logmu} yields, 
\begin{equation} \label{KL1} 
\frac{1}{N} \log(\mu_{N,J}(N G)) = \sum_{j \in spt(G)} \left[- G(j) \log {G(j)} + G(j) \log J(j) \right] + o(N)
\end{equation}
with uniform remainder $\abs{o(N)} \leq 2 (g+1) \log N /N $. Notice that the sum in \eqref{KL1} is equal to $- KL(G,J)$ where $KL(G,J)\geq 0$ is the Kullback-Leibler divergence between $G$ and $J$ given by equation \eqref{KL}, which concludes the proof. 
\end{proof}
\begin{proposition} \label{LDKL}
Fix $0< a <1$, and set $N(a) = 2 / a + \log(F_g)$. On day $n$, let $J_n $ be the population histogram at the end of the mutations phase (Phase 2). For any $G$ with $b(G) > a$, $G' \in B_N(G)$, and any $N > N(a)$, one has the large deviations estimate
\begin{equation} \label{LDselect1}
\frac{1}{N} \log P ( H_{n+1} = G' \; | \;  H_n  ,  R_n  ) = - KL( G, J_n ) + o(N) 
\end{equation}
with uniform remainder $\abs{ o(N) } \leq (5g + 1) \log N / N$.
\end{proposition} 
\begin{proof}
Take $G,G',$ and $N > N(a)$ as stated above. This forces $N > 2/b(G) $ so that $spt(G') = spt(G)$ by Lemma \ref{basicBN}.  Suppose $KL(G,J_n)$ is finite so that $ spt(G') =spt(G) \subset spt(J_n)$. Recall that by construction of the Markov chain $H_n$, we have
$
P ( H_{n+1} = G' \; | \; H_n ,  R_n/ N )= P ( H_{n+1} = G' \; | \; J_n ).
$
Given $J_n$, the conditional distribution of $N H_{n+1}$ is the multinomial $ \mu_{N, J_n}$ given by \eqref{multinomial2}. Since $G'$ is $N$-rational, $P (  H_{n+1} = G'  \; | \; J_n) = \mu_{N, J_n}(N G').$ Then, \eqref{KLlogmu} yields
\begin{equation}\label{h1}
\frac{1}{N} \log P(  H_{n+1} = G' \; | \; J_n ) = - KL(G', J_n) + o(N)
\end{equation}
with $\abs{ o(N) }  \leq  2 (g+1) \log N/N$.  From \eqref{partialKLG}, we get 
\begin{equation}\label{h2}
\abs{KL(G', J_n) - KL(G, J_n)}  \leq  \frac{2 g}{3 N} \left( 1 + \abs{ \log b(G) } + \abs{ \log b(J_n) }\right). 
\end{equation}
We have $b(J_n) \geq 1/(N F_g) $ by \eqref{bJn}. Since $b(G) \geq 2/ N$, the right-hand side of \eqref{h2} is bounded above by $2 g \log N / N$ provided $N > N(a)$. Equation \eqref{h1} then implies 
\begin{equation}\label{h4}
\frac{1}{N} \log P(H_{n+1} = G' \; | \; J_n ) =  - KL(G, J_n) + o_1(N)
\end{equation}
with $\abs{ o_1(N) } \leq 5 g \log N /N$.  Hence, \eqref{LDselect1} is proved when $KL(G, J_n)$ is finite. \smallskip \\
When $KL(G, J_n) = + \infty$, we have $spt(G)\not \in spt(J_n)$. For $G' \in B_N(G)$, one has $spt(G') = spt(G) \not \in spt(J_n)$ so that the transition from $J_n$ to $H_{n+1} = G'$ is impossible during Phase 3. Thus, both sides of \eqref{LDselect1} are equal to $- \infty$.
\end{proof}
\subsection{Large deviations asymptotics for the one-step transition kernel}
With the large deviations analyses completed for each random phase during a daily cycle, we can now complete the full large deviations analysis for the one-step transition kernel.  As population size $N \to \infty$, the Markov transition kernel on the state space of histograms converges to a deterministic kernel at exponentially-fast speed and verifies a precise large deviations principle, as will be seen in this section.  This will be quantified by a composite transition rate, which combines the rate functionals for random mutations and random selection.

Equations \eqref{mutcramer} and \eqref{LDselect1} provide two explicitly computed ``partial'' rate functions, namely, the function $mut(r,H) \geq 0$ controlling large deviations for $P(R_n/N = r \; | \;H_n = H)$ and the Kullback-Leibler divergence $KL(G, J) \geq 0$ controlling large deviations for $P(H_{n+1} = G \; | \; J_n = J).$ We have also computed a deterministic linear function $\Psi(H,r)$ given by \eqref{psi} such that $J_n \approx \Psi(H_n,R_n/N)$ with accuracy approximately $1/N$. We thus expect $KL(G, \Psi(H,r))$ to control large deviations for the conditional probability $P(H_{n+1} = G \; | \; H_n = H, R_n/N = r )$.
\begin{mydef}[Composite Rate Function]
To control large deviations for the composite transition $(H_n = H) \to (R_n/N = r) \to (H_{n+1} = G),$ we introduce the \emph{composite transition rate} $\tau(H,r,G) \geq 0$ defined by
\begin{equation} \label{tauHrG}
\tau(H,r,G) = mut(r,H) + KL(G, \Psi(H,r)) 
\end{equation}
for $H, G \in \mcal{H}$ and $r \in K(H).$ 
\end{mydef}
As discussed previously, $mut(r,H)$ is a finite, continuous, convex function of $r \in K(H)$, and $KL(G,I)$ is a continuous, convex function of $I \in \mathcal{H}$ with strict convexity whenever $KL(G,I)$ is finite. However, $I = \Psi(H,r)) $ is an affine function of $r$ due to \eqref{psi}. Therefore, $\tau(H,r,G) \geq 0$ is continuous and convex in $r$ on the compact, convex closure of $K(H)$ with strict convexity in $r$ whenever $\tau(H,r,G)$ is finite.

When $\tau(H,r,G) < \infty,$ this composite transition rate inherits from $mut(r,H)$ the convexity in $H$ when $(r,G)$ are fixed and convexity in $G$ from $KL(G,J)$ when $(H,r)$ are fixed. Since $mut(r,H)$ is finite, $\tau(H,r,G) = + \infty$ iff $ KL(G, \Psi(H,r)) = + \infty$, which is equivalent to the existence of a genotype $j$ such that $G(j) > 0$ and $ \Psi_j(H,r) = 0$. Due to \eqref{zeroPsi}, we see that for $r \in K(H),$ one has $\tau(H,r,G) = + \infty$ if and only if there is a genotype $j$ such that $G(j) > 0$ and $H(j) = r_{k,j} = r_{j,k}= 0$ for all $k$.  These properties will be key in quantifying large deviations asymptotics for a daily transition, which frequently requires solving an important minimization problem.  This motivates the following definition.
\begin{mydef}[Feasible Transitions and One-Step Cost]
By construction, for any $H, G \in \mcal{H}$, the transition kernel $P(H_{n+1} = G \; | \; H_n = H)$ is strictly positive if and only if for any $j$ such that $G(j) > 0$ and $H(j) =0$, one can find a $k$ such that $H(k) M_{k,j} >0$. We will then say that $(H \to G)$ is a \emph{feasible transition}.  Let $\tau(H,r,G) \geq 0$ be the composite transition rate defined by \eqref{tauHrG}. For any $H,G \in \mathcal(H)$, define the \emph{one-step cost function} $C(H,G) \geq 0$ by 
\begin{equation} \label{costHG}
C(H,G) = \min_{r \in K(H)} \tau(H,r,G) = \min_{r \in K(H)} \left[ mut(r,H) + KL(G, \Psi(H,r) ) \right].
\end{equation}
\end{mydef}
The explicit expression of the one-step cost above is not necessary to obtain a large deviations result for daily transitions, so we will delay the presentation until later where the formula is given in Theorem \ref{thCHG}.  The following lemma and proposition will aid in proving regularity results of the one-step cost function.
\begin{lemma} \label{finitecost}
Let $H, G\in \mcal{H}$.
\begin{enumerate}[label = (\roman*)]
\item The one-step cost $C(H,G)$ is finite if and only if $(H \to G)$ is a feasible transition. 
\item When $C(H,G)$ is finite, there is a unique $\rho$ in the closure of $K(H)$ such that $C(H, G) = \tau(H, \rho, G)$.
\item When $C(H,G)$ is finite, the one-step cost $C(H,G)$ is convex in $H$ for $G$ fixed and convex in $G$ for $H$ fixed.
\end{enumerate}
\end{lemma}
\begin{proof} 
\hfill
\begin{enumerate}[label = (\roman*)]
\item By definition, $C(H,G)$ is finite if and only if there is at least one $r \in K(H)$ with $\tau(H,r,G) $ finite.  As discussed above, this occurs if and only if for each $j$ such that $G(j) >0$ and $H(j) = 0,$ there is a $k$ with $r_{k,j} >0$.  For any $(j,k)$ verifying such a condition, one must also have $M_{k,j} >0$ due to \eqref{Kj} so that $(H \to G)$ is a feasible transition, proving item (i).
\item Since $\tau(H,\rho,G)$ is continuous in $\rho$ on $\overline{K(H)}$, there exists at least one $\rho \in \overline{K(H)}$ such that $C(H, G) = \tau(H, \rho, G)$. When $C(H,G)$ is finite, we have uniqueness of $\rho$ in the compact, convex set $\overline{K(H)}$ since $\tau(H, \rho, G)$ is strictly convex in $\rho$ whenever $\tau(H, \rho, G)$ is finite, proving item (ii).
\item Item (iii) follows from the discussion preceding this lemma.
\end{enumerate}
\end{proof}
\begin{proposition} \label{holdertau}
Fix $A >0, a >0,$ and the parameter set $\mcal{P}$. Define the constant $\eta =\eta(A, a, \mcal{P})$ by
\begin{equation} \label{eta}
\eta  =  3 g^{2+1/a} e^{A/a} F_g / F_1 + 20 F_g \log(F_g) /a.
\end{equation}
For all $H', H, G\in \mcal{H}(a)$ and matrices $r' \in K(H')$ and $r \in K(H)$ verifying
\begin{enumerate}[label = (\roman*)]
\item $ spt(H') = spt(H),$
\item  $\vnorm{ r' - r } \leq 1,$
\item $\tau(H,r,G) \leq A,$
\item $\tau(H',r',G)  < \infty,$
\end{enumerate}
one then has 
\begin{equation} \label{eqholdertau}
\abs{\tau(H',r',G) - \tau(H,r,G)} \leq \eta \left[ \vnorm{ r' - r }^{1/2} + \vnorm{ H'- H } \right].
\end{equation}
\end{proposition}
\begin{proof}
Consider $H',H,G, r',$ and $r$ verifying (i)-(iv) above.  Apply \eqref{holdermut} with H\"{o}lder coefficient $\alpha=1/2$ to get 
\begin{equation}\label{deltamut}
\abs{mut(r', H') - mut(r, H) } \leq c_0  \left[  \vnorm{ r' - r }^{1/2} + \vnorm{ H'- H } \right ] 
\end{equation}
where $c_0 = \frac{ 20}{a} F_g \log(F_g).$  From \eqref{lipPsi}, we get 
\begin{equation} \label{lipPsi2}
\vnorm{ \Psi(H',r') - \Psi(H,r) } \leq \frac{3 g F_g}{F_1} \left( \vnorm{ r' - r } + \vnorm{ H' - H } \right). 
\end{equation}
By definition of $\tau(H,r,G)$ and item (iii) above, one has 
\begin{equation} \label{KL+2}
KL( G, \Psi(H,r) ) \leq \tau(H,r,G) \leq A. 
\end{equation} 
Since $\tau(H', r',G)$ and $\tau(H, r, G)$ are finite, both $spt(\Psi(H,r))$ and $spt(\Psi(H',r'))$ contain $spt(G)$. Therefore, we can apply \eqref{KL+2} and \eqref{lipKLJI} to get 
$$
\abs{KL( G, \Psi(H',r')) - KL(G, \Psi(H,r))}  \leq  c_1  \vnorm{ \Psi(H',r') - \Psi(H, r) }
$$
with $c_1 = g^{1+1/a} e^{A/a}.$  In view of \eqref{lipPsi2}, this yields
$$
\abs{ KL( G, \Psi(H',r') ) - KL( G, \Psi(H,r) ) }  \leq  c_2  \left[  \vnorm{ r' - r } + \vnorm{ H'-H } \right]
$$
with $c_2 = 3 c_1 g F_g / F_1 $.  Combining this last result with \eqref{deltamut} and setting $\eta = c_0 + c_2$ yields
$$
\abs{\tau(H',r',G) - \tau(H,r,G) }\leq \eta \left[\vnorm{ r' - r }^{1/2} + \vnorm{ H'-H } \right].
$$
Noting that $\eta$ here is precisely \eqref{eta} concludes the proof.
\end{proof}
We can now prove uniform continuity results for the one-step cost function $C(H,G)$, which will be required to prove Theorem \ref{transLD}.
\begin{theorem} \label{holderCost}
Fix $0< a < 1, d >0,$ and the parameter set $\mcal{P}$. Consider any histograms $H, G, H',G' \in \mcal{H}(a)$ verifying 
\begin{enumerate}[label = (\roman*)]
\item $spt(G') = spt(G),$ 
\item $spt(H) = spt(H'),$
\item $\vnorm{ H' - H } \leq a/F_g,$
\item $C(H,G) \leq d.$
\end{enumerate}
There exists a constant $c = c(d, a, \mcal{P})$ such that 
\begin{equation} \label{continuityCHG}
\abs{ C(H',G') - C(H,G) } \leq c  \left(\vnorm{ H' - H }^{1/2} + \vnorm{ G' - G }\right). 
\end{equation}
\end{theorem}
\begin{proof}
Consider histograms $G, G', H\in \mcal{H}(a)$ verifying $(i)$ and $(iv)$ above.
One can then select $r \in K(H)$ with $\tau(H,r,G) \leq D+1$ so that $I = \Psi(H,r)$ must verify $spt(G') = spt(G) \subset spt(I)$ and $ KL(G, I) \leq D+1$. Apply \eqref{holderKLG} to get the Lipschitz bound 
\begin{equation} \label{KLGG'}
\abs{KL(G',I) - KL(G,I)} \leq c_1  \vnorm{ G' - G }  
\end{equation} 
with $c_1 = \frac{g}{a} (3 + \log g  + D ).$ This implies, by definition of $\tau(H,r,G)$,
$$
\abs{\tau(H,r,G') - \tau(H,r,G)} = \abs{KL(G', I) - KL(G, I)} \leq c_1 \vnorm{ G' - G }. 
$$
Since $ C(H,G') \leq \tau(H,r,G'),$
\begin{equation} \label{Ctau}
C(H,G') \leq \tau(H,r,G) + c_1 \vnorm{ G' - G } \leq d_1
\end{equation} 
with $d_1 = d + 1 + c_1$. In \eqref{Ctau}, take the infimum of the middle term over $r \in K(H)$ to obtain
\begin{equation} \label{CCGG'}
C(H,G') \leq C(H,G) + c_1  \vnorm{ G' - G }. 
\end{equation} 
Define $c_2= \frac{g}{a} (3 + \log g + d_1)$. Since $ C(H,G') \leq D_1$, the generic result \eqref{CCGG'} can now also be rewritten by switching the roles of $G$ and $G'$ provided one also replaces $d$ with $d_1$, and $c_1$ with $c_2$.  This yields $C(H,G) \leq C(H,G') + c_2 \vnorm{ G' - G }$. Therefore, the bound in \eqref{CCGG'}yields
\begin{equation} \label{lipCHGG'}
\abs{ C(H,G) - C(H,G') } \leq c_2  \vnorm{ G' - G }.
\end{equation}

Consider now any $H', H, G \in \mcal{H}(a)$ verifying $(ii)$-$(iv)$ above. Define a linear mapping $\mcal{L}$ for all $\rho \in K(H)$ denoted $\hat{\rho} := \mcal{L}(\rho)$ given by
\begin{equation} \label{biject}
\hat{\rho}_{j,k}  = 
\begin{cases}
\rho_{j,k} H(j) / H'(j), & j \in spt(H) \text{ and any } k, \\
 0, &  \text{otherwise.}
\end{cases} 
\end{equation} 
By definition of $K(H')$ and $K(H)$, one readily verifies that $\hat{\rho} \in K(H)$. Thus, $\mathcal{L}$ maps $K(H')$ into $K(H)$ bijectively with an inverse mapping $ \hat{\rho} \to \rho$ defined by exchanging $\rho$ and $\hat{\rho}$ as well as $H'$ and $H$ in \eqref{biject}.  Definition \eqref{biject} forces $spt(\rho) = spt(\hat{\rho})$. Since $spt(H')= spt(H)$, then \eqref{sptpsiHh} yields $spt(\Psi(H',\rho))= spt(\Psi(H,\hat{\rho}))$. Hence, either these two supports contain $spt(G)$ or neither of them do. This implies the equivalence
\begin{equation} \label{bothfinite}
\tau(H',\rho,G)< \infty \iff \tau(H,\hat{\rho},G) < \infty.
\end{equation} 
For $\rho \in K(H')$, the bound $\vnorm{ \rho } \leq F_g$ holds due to \eqref{Kj}.  Thus, we get $\abs{ 1- H(j)/H'(j) } \leq \vnorm{ H' - H } / a $ for $j \in spt(H') = spt(H)$. Hence, for all $\rho \in K(H'),$
\begin{equation} \label{rohatro}
\vnorm{ \rho - \hat{\rho} }  \leq \vnorm{ \rho}\, \vnorm{ H' - H } / a  \leq  \frac{F_g}{a}  \vnorm{ H' - H } \leq 1.
\end{equation} 
Define $U =\vset{ r \in K(H) \,|\, \tau(H,r,G) \leq D}$.  This implies $C(H,G) = \min_{r \in U} \; \tau(H,r,G) $ since $C(H,G) \leq D$. For $\hat{\rho} \in U$, the cost $\tau(H, \hat{\rho},G)$ is finite so that $\tau(H',\rho,G)$ is also finite due to \eqref{bothfinite}. Apply Proposition \ref{holdertau} to get the constant $\eta= \eta(D,a,PAR)$ given by \eqref{eta} such that 
$$
\abs{ \tau(H',\rho,G) - \tau(H,\hat{\rho},G) }  \leq  \eta \left(\vnorm{\rho - \hat{\rho}}^{1/2} + \vnorm{ H'-H }\right).
$$ 
For $\hat{\rho} \in U$, the bound in \eqref{rohatro} yields
\begin{equation} \label{tauHh}
\abs{\tau(H',\rho,G) - \tau(H,\hat{\rho},G) }  \leq  \eta_0 \vnorm{ H'- H }^{1/2}
\end{equation} 
with $\eta_0 = 2 (F_g / a)^{1/2} \eta $.
Since $\mathcal{L}$ is a bijection and $U\subset K(H)$, any $r \in U$ is of the form $r = \hat{s}$ for some $s = s(r) \in K(H')$. Then, for any $r \in U$, applying \eqref{tauHh} to $\rho = s$ and $r = \hat{s}$ implies 
$$
C(H',G) \leq \tau(H', s ,G)  \leq  \tau(H, r ,G) + \eta_0  \vnorm{ H'- H }^{1/2}.
$$
Take the infimum of the right-hand side over all $r \in U$ to get 
\begin{equation}\label{ChCH}
C(H',G) \leq  C(H,G) + \eta_0 \vnorm{ H'- H }^{1/2} \leq D + \eta_0.
\end{equation}
Define $V = \vset{\rho \in K(H) \, | \, \tau(H',\rho,G) \leq D + \eta_0}$. This implies $C(H',G) = \min_{\rho \in V} \; \tau(H',\rho,G)$.  For $\rho \in V$, the cost $\tau(H', \rho,G)$ is finite so that $\tau(H,\hat{\rho},G)$ is also finite. Apply again Proposition \ref{holdertau} with $A= D + \eta_0$ to get $\eta_1 = \eta(d +\eta_0, a, \mcal{P})$ such that, with \eqref{rohatro}, 
\begin{equation} \label{tauHh2}
\abs{\tau(H',\rho,G) - \tau(H,\hat{\rho},G)}  \leq  \eta_1  \left(\vnorm{ \rho - \hat{\rho} }^{1/2} + \vnorm{ H' -H }\right)  \leq  \eta_2 \vnorm{ H' - H }^{1/2}
\end{equation} 
with $\eta_2 = 2 (F_g / a)^{1/2} \eta_1 $.  This implies, for all $\rho \in V,$
$$
\tau(H',\rho,G) \geq \tau(H,\hat{\rho},G) - \eta_2 \vnorm{ H' - H }^{1/2} \geq C(H,G) - \eta_2 \vnorm{ H' - H }^{1/2}.
$$
Take the infimum of the left-hand side over $\rho \in V$ to get 
$
C(H',G) \geq C(H,G) - \eta_2 \vnorm{ H' - H }^{1/2}.
$
Combine this with \eqref{ChCH}, and set $\gamma = \max\vset{\eta_2, \eta_0}$ to obtain
\begin{equation} \label{holderCHhG}
\abs{C(H',G) - C(H,G)} \leq \gamma \vnorm{ H' - H }^{1/2}.
\end{equation}
Finally, consider histograms $H,G,H',G' \in \mcal{H}(a)$ verifying $(i)$-$(iv)$ in the statement of the theorem. Then, as shown previously, 
\begin{equation} \label{holderH'HG}
\abs{C(H',G) - C(H,G)} \leq \gamma \vnorm{ H' - H }^{1/2}.
\end{equation}
This forces $C(H',G) \leq d + \gamma.$ Apply then \eqref{lipCHGG'} to obtain $c_3 = c_2(d + \gamma, a, \mcal{P})$ such that 
$$
\abs{C(H',G) - C(H',G')} \leq c_3  \vnorm{ G' - G }.
$$
Combining this with \eqref{holderH'HG} yields \eqref{continuityCHG} with a constant $c = c_3 + \gamma.$
\end{proof}
A final ingredient needed to obtain a large deviations result for daily transitions is the following lemma.
\begin{lemma} \label{integral}
Fix the parameters $\mcal{P}$ and $0< a < 1$. Then, $a$ and  $\mcal{P}$ determine constants $c_0$ and $N_0$ with the following properties. Let $(H,r,G)$ be any $N$-rational triple with $H,G \in \mcal{H}(a)$ and $r \in K_N(H)$. Given $H_n= H $ and $R_n/N = r$, the population histogram $J$ after mutations is a deterministic function $J = J(H, r)$ defined by \eqref{WWn}.  For $N > N_0$ and $r \in K_N(H),$ the histograms $I = \Psi(H,r)$ and $J = J(H,r)$ satisfy $spt(J)= spt(I)$ with
\begin{equation} \label{compareJI}
\vnorm{J - I} \leq  c_0 / N. 
\end{equation}
Let $E_N = E_N(H,G) = \vset{ r \in K_N(H) \; | \; KL( G, J )< \infty \text{ and }  KL(G, I)< \infty}$.  For $N > N_0$ and $r \notin E_N$, one has $KL(G, J) = KL(G, I) = \infty$. For $N > N_0$ and $r \in E_N$, one has
\begin{equation} \label{**} 
\abs{KL( G, J ) - KL( G, I)}  \leq  c_0  e^{\kappa / a} \vnorm{ J - I } 
\end{equation} 
with $\kappa = \kappa(H,r,G) = \min \vset{KL(G,J) , KL(G,I)}$.  For $r \in E_N$, define $T(H,r,G) = mut(r, H) + KL( G, J(H,r) ).$  Then, for $ N> N_0$, the Markov transition kernel $\mathcal{Q}(H,G)$ verifies 
\begin{equation} \label{trans4}
\mathcal{Q}(H, G)  =  \sum_{r \in  E_N} \exp \left(-N[T(H,r,G) + o(N)]\right) 
\end{equation}
with $\abs{ o(N) } \leq c_0  \log N / N.$
\end{lemma}
\begin{proof}
Consider any $N$-rational triple$(H,r,G)$ as above. Given $H_n= H $ and $R_n/N = r $, the histogram $J = J(H, r)$ after mutations is given by \eqref{WWn}. Apply \eqref{JJn} and the discussion in Section \ref{JnMut} to $I = \Psi(H,r)$ and $J = J(H,r)$ to obtain \eqref{compareJI} for $N > N_1 = 20 + 2 F_g / a$ and $c = 13 F_g / F_1$.

For $N>N_1$, either $ spt(G) \subset spt(I) = spt(J) $ so that $KL(G,I)$ and $KL(G,J)$ are both finite, or $spt(G)$ is not included in $spt(I) = spt(J) $ so that $KL(G,I) = KL(G,J) =  \infty$. Define $\kappa$ and $E_N$ as in the statement of the lemma. Set $c_1 = g^{1+1/a}$ and $c_2= c_1 \, c$. Combine \eqref{lipKLJI} with \eqref{compareJI} to obtain
$
\abs{KL( G, J ) - KL( G, I)}  \leq  c_1  e^{\kappa / a} \vnorm{ J - I } \leq c_2/ N 
$
for $N > N_1$ and $r \in E_N.$  This proves \eqref{**}.
Finally, the transition kernel $\mathcal{Q}$ verifies by construction
\begin{equation} \label{trans0}
\mathcal{Q}(H, G)  =  \sum_{r \in K_N(H)} P( R_n/N = r \; | \; H_n = H )\,  P( H_{n+1} = G \; | \; H_n = H, \, R_n/N = r ).
\end{equation}
where the set $K_N(H)$ of $N$-rational matrices $r \in K(H)$ satisfies $card(K_N(H)) \leq [ (N+1) F_g ]^{g^2}$ since $\vnorm{ r } \leq F_g.$  From \eqref{LDselect1} we get, for $N > N_2 = N_1 + 2 / a + \log F_g$, 
\begin{equation} \label{trans2}
P(H_{n+1} = G \; | \; H_n = H, \, R_n/N = r) = \exp\left( - N  [ KL(G, J(H, r) )+ o_1(N) ]\right)
\end{equation}
with $\abs{o_1(N)} \leq (5 g+1) \log N / N$.  In \eqref{trans2}, the right-hand side is $0$ unless $r \in E_N$. Therefore, the sum in \eqref{trans0} can be restricted to $r \in E_N$. Then, \eqref{LDmut1} yields, for $N > N_3 = N_2 + [ 4 + 2 \log(g) ] / a ,$
\begin{equation} \label{trans3}
P( R_n / N = r \; | \;  H_n = H ) = \exp\left(- N [ mut(r, H) + o_2(N)]\right) 
\end{equation} 
with $\abs{o_2(N)} \leq 4 g^2  \log N / N$.  For $r \in E_N$, define $T(H,r,G) = mut(r, H) + KL( G, J(H,r) )$. For $N > N_3$, substitute \eqref{trans2} and \eqref{trans3} into \eqref{trans0} to prove \eqref{trans4} with $o(N) = o_1(N) + o_2(N).$
\end{proof}

Finally, we are ready to state the important large deviations result for the one-step transition kernel controlling daily transitions.  The following theorem will be key later when we lift these results to the large deviations analysis of random histogram trajectories (see Section \ref{LDPathTheory}).

\begin{theorem} \label{transLD}
Fix any $0< a < 1 < d$ and the parameters $\mcal{P}$. One-step large deviations for the Markov chain $H_n$ are controlled as follows by two constants $c=c(d, a, \mcal{P})$ and $N_0=N_0(d, a, \mcal{P})$.  Consider any $N$-rational histograms $H, G \in \mcal{H}(a)$ with transition cost $C(H,G) \leq d$. Then, the transition kernel $\mathcal{Q}(H,G)$ has a uniform large deviations approximation, valid for all $N >N_0$ and $H,G\in \mcal{H}(a)$ as above, 
\begin{equation} \label{onestepLD}
\frac{1}{N} \log \mathcal{Q}(H, G) = -  C(H,G) + o(N)
\end{equation}
with  $\abs{ o(N) } \leq c / \sqrt{N}.$
\end{theorem}
\begin{proof}
Consider any $H,G\in \mcal{H}(a)$ with $C(H,G) < d$. All $c_i$ and $N_i$ defined below will be constants depending only on $d,a,$ and $\mcal{P}$. For $r \in K_N(H)$, let $J = J(H,r)$ be the population histogram after mutations given $H_n= H$ and $R_n/N = r$. Let $I = \Psi(H,r)$. Denote $\kappa(H,r,G) = \min\vset{KL( G, J) , KL( G, I)}.$ Let
$ T(H,r,G) = mut(H,r) + KL(G, J)$ and $\tau(H,r,G) = mut(H,r) + KL(G, I)$. Let $E_N$ be as in Lemma \ref{integral}. Then, Lemma \ref{integral} provides constants $c_0$ and $N_0$ such that for $N >N_0$,
\begin{align} 
spt(J) = spt(I) \;& \; \; \text{and} \; \; \vnorm{  J - I } \leq  c_0/ N \; \; \text{for} \; r \in K_N(H), \label{recall1}\\
\abs{KL( G, J ) - KL( G, I)} \;& \leq \; c_0 e^{\kappa / a} \vnorm{ J - I } \; \; \text{for} \; r \in E_N. \label{recall2}
\end{align} 
Furthermore, 
$$
\mathcal{Q}(H, G)  =  \sum_{r  \in  E_N}  \exp (- N [T(H,r,G) + o(N) ]),
$$
with $\abs{o(N)} \leq c_0  \log N / N$ for $N > N_0$.  Let $U_N = \exp(- N [ T(H,r,G) + o(N) ])$. Partition $E_N$ into two subsets $\mathcal{M}$ and $\mathcal{R}$ defined by 
\begin{equation} \label{MM} 
\mathcal{M} = \vset{ r \in E_N \; | \; KL( G, J) > 4 d} \qquad \text{and} \qquad  \mathcal{R} = \vset{ r \in E_N \; | \; KL( G, J) \leq 4 d}
\end{equation} 
so that $\mathcal{Q}(H, G)  =  S(\mathcal{R})+ S(\mathcal{M})$ with $S(\mathcal{R})  =  \sum_{r  \in  \mathcal{R} }  U_N$ and $S(\mathcal{M})  =  \sum_{r  \in  \mathcal{M} } U_N.$  Set $N > N_1 = N_0 + (c_0/ 2 D)^2$ to write, for $r \in \mathcal{M}$ and $N> N_1$, 
$$
U_N \geq KL(G, J) - c_0 \log N /N \geq 4 d - c_0 \log N / N \geq 2d.
$$ 
Since $card(\mathcal{M}) \leq (N+1)^{g^2}$, this yields 
\begin{equation} \label{S0bound}
S(\mathcal{M}) \leq (N+1)^{g^2} \; e^{- 2 N d } \leq e^{- N d }
\end{equation} 
for $N > N_2 = N_1 + (2 g^2/ d)^2.$ For $r \in \mathcal{R}$, one has $\kappa < KL( G, J ) \leq 4d$ , and \eqref{recall2} yields for $N > N_2$ 
$$
\abs{KL( G, J ) - KL( G, I)}  \leq  e^{\kappa / a}(c_0^2 / N) \leq c_2 / N \
$$
with $c_2 = c_0^2  e^{4d/a}.$ This implies $\abs{T(H, r, G) - \tau(H, r, G)} \leq  c_2 / N$ for $r \in \mathcal{R}$ and $N > N_2$ . This forces $U_N = \tau(H,r,G) + o_1(N)$ with $\abs{o_1(N) } \leq \abs{o(N) }+ c_2 / N \leq c_3  \log(N) / N$ with $c_3 = c_0 + c_2$.  By definition of $S(\mcal{R}),$ this implies for $N > N_2$ that
\begin{equation}\label{S3}
S(\mathcal{R}) = \sum_{r  \in  \mcal{R} }  \exp (- N  [ \tau(H,r,G) + o_1(N) ] ).
\end{equation}
By definition, one has $\tau(H,r,G) \geq C(H,G)$ so that \eqref{S3} gives
\begin{equation} \label{SS+}
S(\mathcal{R})  \leq  (N+1)^{g^2} \exp (- N \,C(H, G) + c_3 \log(N)) \leq \exp (- N \,C(H, G) + c_4 \log N )
\end{equation} 
for $N > N_2$ and some constant $c_4 > c_3.$  Equation \eqref{S0bound} gives $ S(\mathcal{M}) \leq e^{- N d} \leq e^{- N C(H,G)} $, since $C(H,G) \leq D$. Combine this with \eqref{SS+} to get
$$
\mathcal{Q}(H, G) = S(\mathcal{M}) + S(\mathcal{R}) \leq e^{- N \,C(H,G)} ( 1 + e^{c_4 \log{N}} ) \leq \exp (  - N \,C(H,G) + c_5 \log{N})
$$
for $N > N_2$ and $c_5 = 3 c_4.$  This yields the large deviations upper bound 
\begin{equation} \label{upLD}
\frac{1}{N} \log \mathcal{Q}(H, G) \leq - C(H,G) + c_5 \log N /N. 
\end{equation}
By definition of $C(H,G)$, there exists a matrix $w= w(H,G,N) \in K(H)$ such that 
\begin{equation}\label{wHGN}
C(H,G) \leq \tau(H,w,G) \leq C(H,G) + 1/ \sqrt{N}.
\end{equation}
Lemma \ref{density} shows that for $N > N_3 =N_2 + \frac{ g^2}{a F_1}$, there is an $N$-rational matrix $s \in K_N(H)$ such that $spt(s) = spt(w)$ and $\vnorm{ s - w } \leq g/N.$   Equation \eqref{wHGN} implies $\tau(H,w,G) \leq d + 1$ so that $\tau(H,s,G)$ must be finite since $ spt(s) = spt(w)$. A fortiori $KL(G, \Psi(H,s))$ is finite, implying $spt(\Psi(H,s)) \subset spt(G)$. But since $s \in K_N(H)$, equation \eqref{recall1} forces $spt(J(H,s)) = spt(\Psi(H,s))$, so that $KL(G, J(H,s))$ must also be finite.  Therefore, $s \in E_N$.  Apply \eqref{eqholdertau} to get $c_6 = \eta(d+1, a, \mcal{P})$ such that
\begin{equation}\label{tau.sw}
\abs{\tau(H,s,G) - \tau(H,w,G)} \leq c_6 \vnorm{ s-w }^{1/2} \leq c_6 \sqrt{g} / \sqrt{N}.
\end{equation}
Set $c_7 = 1 + c_6 \sqrt{g}$ and $N_4 = N_3 + c_7^2$. Using \eqref{wHGN}, we get for $N>N_4$
\begin{equation} \label{tauHsG}
\tau(H,s,G) \leq \tau(H,w,G) + c_6 \sqrt{g} ) / \sqrt{N} \leq C(H,G) + c_7 / \sqrt{N} \leq d+1.
\end{equation}
Then, $\kappa(H, s, G) \leq KL(G,\Psi(H,s)) \leq \tau(H,s,G)) \leq d + 1$.  For $N> N_4$, since $s \in E_N$, apply \eqref{recall2} to the triple $(H,s,G)$ to get 
\begin{equation} \label{***} 
\abs{KL( G, J(H,s) ) - KL( G, \Psi(H,s) ) }  \leq  c_8 /N   
\end{equation}
with $c_8 = c_0^2 e^{(d+1)/a}.$ Since $d \geq 1$, this implies
$$
KL( G, J(H,s) ) \leq KL( G, \Psi(H,s) ) + c_8 /N \leq d + 1 + c_8 /N < d + 2 < 4 d 
$$ 
for $N > N_5 = N_4 + 1/c_8.$ Hence for $N > N_5$, the matrix $s \in E_N$ must belong to $\mathcal{R}$, and \eqref{S3} provides the lower bound 
\begin{equation}\label{S4}
S(\mathcal{R}) = \sum_{r \in \mathcal{R} } \exp (  - N  [ \tau(H,r,G) + o_1(N) ]  ) \geq \exp ( - N  [ \tau(H,s,G) + o_1(N) ] ). 
\end{equation}
Due to the bound on $o_1(N)$ and \eqref{tauHsG}, 
$$
\frac{1}{N} \log S(\mathcal{R}) \geq - \tau(H,s,G) +o_1(N) \geq - C(H,G) - c_7 / \sqrt{N} - c_3 \log N /N 
$$
for $N > N_5.$  Finally, since $\mathcal{Q}(H, G) \geq S(\mathcal{R})$, setting $c_9 = c_3 + c_7$ yields
$$
\frac{1}{N} \log \mathcal{Q}(H, G) \geq - C(H,G) - c_9 / \sqrt{N}.
$$
In view of \eqref{upLD}, we now obtain for
\begin{equation} \label{trans}
- C(H,G) - c_9 / \sqrt{N} \leq  \frac{1}{N} \log \mathcal{Q}(H, G)  \leq - C(H,G) + c_5 / \sqrt{N}
\end{equation}
for $N > N_5,$ which concludes the proof.
\end{proof}



\subsection{Computation of the one step transition cost}
Since the main goal of this paper is to present an application of large deviations theory to computing most likely evolutionary paths linking an initial histogram to a desired target histogram, obtaining an expression for the one-step cost function will be important in Section \ref{mostlikely}.  By definition of the one-step cost given by \eqref{costHG}, one needs to minimize the convex function $\tau(H,r,G)$ over all $g \times g$ matrices $r \in K(H).$  Recall that matrices in the convex set $K(H)$ satisfy
\begin{align} 
r_{j,k} & \geq 0 \qquad \text{for all $(j,k)$,}\label{cons0} \\
\sum_k  r_{j,k} & <  Fj H(j) \qquad \text{for all $j$,}  \label{cons1} \\
r_{j,k} & =  0 \qquad \text{when $Q_{j,k} = 0$.}\label{cons2}
\end{align}
The interior $K(H)^\circ$ of $K(H)$ is the set of all $r$ verifying \eqref{cons1},  \eqref{cons2}, and  $r_{j,k}> 0$ whenever $Q_{j,k} >0$.  The following theorem gives the explicit computation of this one-step cost function.
\begin{theorem}\label{thCHG}
Fix the parameters $\mcal{P}$ and any $0 < a < 1$. Let $\Gamma(a)$ be the set of interior histograms $J$ with $b(J) > a$. 
There is a constant $c = c(a,\mcal{P}) > 0$ such that for all $H, G \in \Gamma(a)$ and $0 \leq m < c$, the transition cost $C(H,G)$ is a finite $C^{\infty}$ function of $(m, H, G)$. Moreover, $C(H,G)$ has an explicit first-order expansion in $m$ given by
\begin{equation}\label{formulaCHG}
C(H, G) =  KL(G,\Phi) + m \sum_{j,k}  F_j H(j) Q_{j,k} [1 - U_k/U_j ] + O(m^2)
\end{equation}
where  
\begin{align*} 
KL(G,\Phi) = \sum_j  G(j) \log(G(j) / \Phi(j)) > 0, \\
U_j =  \exp \left( \frac{G(j)}{F_j H(j)}\right), \qquad \Phi_j = F_j H(j) / \langle F,H\rangle.
\end{align*}
\end{theorem}
\begin{proof}
Recall that $\tau(H,r,G) = mut(H,r) + KL(G, \Psi(H,r)$ is given by, for $r \in K(H)$,
\begin{align} 
\Psi_i = \Psi_i(H,r) & =  \frac{1}{\langle F, H\rangle} \left( F_i H(i) - \sum_k  r_{i,k} + \sum_k r_{k,i} \right),  \label{eqpsi}\\
KL(G, \Psi(H,r)) & = \sum_i  G(i) \log ( G(i) / \Psi_i ),  \label{eqKL}\\
mut(H,r) & = \sum_{(j,k) \in spt(Q)} \left[ r_{j,k} \log r_{j,k} - r_{j,k} - rjk \log (m f_{j,k})\right].  \label{eqmut}
\end{align}
To minimize $\tau(H,r,G)$, we only need to consider $r \in K(H)$ with $\tau(H,r,G)$ finite, which holds iff $\Psi_i > 0$ for all $i$. We seek to minimize $\tau(H,r,G)$ over $r =r(m)\in K(H)^\circ$. Such an $r$ must verify 
\begin{equation} \label{sys}
\frac{\partial}{\partial r_{j,k}} \tau(H,r,G) = 0, \qquad Q_{j,k} >0.
\end{equation}
For $(j,k) \in spt(Q),$ basic derivations and algebra reduce \eqref{sys} to the system 
\begin{equation}\label{Dtau}
0 =  \log\frac{r_{j,k}}{m f_{j,k}} - \sum_i \; \frac{G(i)}{\Psi_i}\frac{\partial \Psi_i}{\partial r_{j,k}}.
\end{equation}
The coefficients $p(i,j,k) = \frac{\partial \Psi_i}{\partial r_{j,k}}$ are constants given by
\begin{equation}\label{DPsi}
\langle F,H \rangle  p(i,j,k) = 
\begin{cases} 
-1, &  i =j  \neq k,\\
1, &   j \neq k = i, \\
0 &  \text{otherwise.}
\end{cases}
\end{equation}
Let $q = card(spt(Q)).$ Define the vector $ x \in (R^+)^q $ by $x_{j,k} = r_{j,k}/m f_{j,k}$ for $Q(j,k) >0$. Substitute \eqref{DPsi} into \eqref{Dtau} to get  
\begin{equation}\label{xsolv}
\log x_{j,k} = - G(j)/A_j + G(k)/A_k,
\end{equation}
where the $A_i =A_i(m,x)$ are given by 
\begin{equation*}
A_i(m,x)  = F_i H(i) + \sum_k  m \left[ -  f_{i,k}\, x_{i,k} +  f_{k,i}\,x_{k,i} \right]. 
\end{equation*}
To enforce  \eqref{cons1}, the vector $x$  must verify, for all $1 \leq i \leq g,$ the set of strict linear constraints 
\begin{equation}\label{consx}
F_i H(i) - m \sum_k f_{i,k} x_{i,k} > 0.
\end{equation}
Call $\mathcal{O}$ the open set of all $x \in (R^+)^q$ verifying \eqref{consx}. Rewrite \eqref{xsolv} as an implicit equation $W(m,x) = 0$ for $x \in \mathcal{O}$ where the function $W(m,x)$ is given by
$$
W_{j,k}(m,x) = x_{j,k} - \exp\left( - G(j)/A_j + G(k)/A_k \right) 
$$
for all $(j,k) \in spt(Q), x\in \mcal{O}, $ and $m \geq 0.$ Then $W(m,x)$ is of class $C^{\infty}$ in $(m,x)$, as well as in $H,G \in \mathcal{H}^o$. \\
For $m=0$ and $(j,k) \in spt(Q)$, the system $W(m,x) = 0$ has a unique solution $\wt{x} \in \mathcal{H}^\circ$ given by
\begin{equation} \label{xsolution}
\wt{x}_{j,k} = \exp\left(- G(j)/[ F_j H(j) ] + G(k)/ [ F_k H(k) ]\right).
\end{equation}
For all $(s,t) \in spt(Q)$ and all $i$, the derivatives $a(i,s,t) = \frac{\partial A_i}{\partial x_{s,t}}(m,x)$ verify $a(i,s,t) =m b(i,s,t)$ where the only nonzero terms of $b(i,s,t)$ are $b(i,i,t) = -f_{i,t}$ and  $b(i,s,i) = f_{s,i}$.  The Jacobian $JAC(m,x)$ of $W(m,x)$ has coefficients $\frac{\partial W_{j,k}}{\partial x_{s,t}}$ given for all $(j,k)$ and $(s,t)$ in $spt(Q)$ by 
$$
\frac{\partial W_{j,k}}{\partial x_{s,t}} = Id_{(j,k) \, ; \, (s,t)} - m Y_{j,k} \exp\left(- G(j)/A_j + G(k)/A_k\right) 
$$
where $Y_{j,k} = b(j,s,t) G(j) / A_j^2 - b(k,s,t) G(k) /A_k^2,$ and $Id$ is the $q \times q$ identity matrix.  At the point $(0,\wt{x}),$ the Jacobian  $JAC(0,\wt{x})$ is hence equal to Id and thus invertible. The classical implicit function theorem then applies to $W(m,x) = 0$ and provides $c = c(a, \mcal{P}) > 0$ such that the equation $W(m,x) = 0$ has a unique solution $x(m) \in \mathcal{O}$ for $m \leq c$  and $H,G \in \Gamma(a)$. The same theorem implies that $x(m) = x(m,H,G)$ is of class  $C^{\infty}$ in $(m,H,G) $.
Define $r(m)$ by
\begin{equation*}
r_{j,k}(m) =
\begin{cases}
m x_{j,k}(m),&  Q(j,k) > 0,\\
0,&  \text{otherwise.}
\end{cases}
\end{equation*}
Then, $r(m)$ is a solution of \eqref{sys} and inherits from $x(m)$ the $C^{\infty}$ smoothness in $(m,H,G)$. Moreover, $r(m)$ verifies the constraints \eqref{cons0}--\eqref{cons2} since $x(m) \in \mathcal{O}$. The positivity of all coordinates of $x(m)$ implies $r_{j,k}(m) >0$ for $Q(j,k) >0$. 
Hence, $r(m) \in K(H)^\circ$ for $m< c$ and solves \eqref{sys}. The strict convexity of $\tau(H,r,G)$ on the open convex set $K(H)^\circ$ forces $r(m)$ to be the unique minimizer of $\tau(H,r,G)$ for $r \in K(H)$ so that $C(H,G) = \tau(H,r(m),G)$.  The function $(H,r,G) \to \tau(H,r,G)$ is $C^{\infty}$ for $H,G \in \mathcal{H}^\circ$ and $r$ in $K(H)^o$. Hence the function $(m, H,G) \to C(H,G)$ is also $C^{\infty}$ for $m < c$ and  $(H,G) \in \Gamma(a)$.
To now obtain the first-order expansion, for all genotypes $j$ and $k,$ define 
\begin{equation} \label{UjEjk}
U_j  =  \exp \left( G(j)/ [ F_j H(j) ]\right) \qquad \text{and} \qquad  E_{j,k} =  U_k / U_j
\end{equation}
so that $\wt{x}_{j,k} = f_{j,k} E_{j,k} = F_j H(j) Q_{j,k} U_k/U_j$. To simplify notation, write $u \simeq v$ whenever $u(m) = v(m) + O(m^2)$ as $m \to 0$.  As $m \to 0,$ the differentiability of $x(m)$ gives   
\begin{equation}\label{rsolv}
r_{j,k}(m) =  m x_{j,k}(m) \simeq m \wt{x}_{j,k} =  m  f_{j,k} E_{j,k}.
\end{equation}
Inserting \eqref{rsolv} into \eqref{eqmut} yields $mut(H,r(m)) \simeq  m \mu$ with 
$
\mu = \sum_{j,k}  f_{j,k} \left(1 - E_{j,k} + E_{j,k}\log E_{j,k} \right). 
$
Since $\log(E_{j,k}) = - G(j)/ [ F_jH(j) ] + G(k)/ [ F_k H(k)]$, we have 
\begin{equation} \label{eqmu}
\mu = \sum_{j,k} f_{j,k}( 1 - E_{j,k}) + \sum_{j,k}  Q_{j,k} E_{j,k}\left(- G(j)  + G(k) \frac{F_j H(j)}{F_k H(k)}\right). 
\end{equation}
Substitute \eqref{rsolv} in \eqref{eqpsi} to get
\begin{align*}
\Psi_j &=  \Psi_j(H,r(m)) \simeq  (F_j H(j) - m t_j) / \langle F,H\rangle; \\
t_j &= \sum_k  F_j H(j) Q_{j,k} E_{j,k} -  \sum_k  F_k H(k) Q_{k,j}E_{k,j}.
\end{align*}
This yields $\log(\Psi_j) \simeq  \log\left(F_j H(j) /\langle F,H\rangle\right) -  m \frac{t_j}{F_j H(j)}.$  Recall that $\Phi_j(H) = F_j H(j)/\langle F,H\rangle$ is the population histogram at the end of the daily deterministic growth starting with histogram $H$.  Substitute the expansion of $\log(\Psi_j)$ in \eqref{KL} to obtain
\begin{align*}
&KL(G, \Psi(H, r(m)) \simeq  \kappa + m \eta;\\
\kappa = KL(G, \Phi) =  \sum_j  G(j) &\log\left(\frac{G(j)}{\Phi_j}\right)  > 0; \qquad  \qquad \eta  = \sum_j  \frac{ G(j) t_j}{F_j H(j)}.
\end{align*} 
This implies $\eta  = \sum_{j,k}  \left[G(j) Q_{j,k} E_{j,k} - \frac{G(j)}{F_jH(j)} F_k H(k) Q_{k,j}\right] E_{k,j}$.  Exchange $j$ and $k$ in the second term  of the previous sum to get 
\begin{equation}
\label{eta1}
\eta = \sum_{j,k} Q_{j,k} E_{j,k}\left[G(j) - G(k) \frac{F_j H(j)}{F_k H(k)}\right].
\end{equation}
Hence, $C(H,G) = \tau(H,r(m),G) = \kappa + m (\eta +\mu) + O(m^2)$ is the first-order expansion of $C(H,G)$.  Combine \eqref{eqmu} and \eqref{eta1} to obtain
$
\eta +\mu = \sum_{j,k} F_j H(j) Q_{j,k} \left( 1 - U_k/U_j \right),
$
which concludes the proof of \eqref{formulaCHG} and the theorem.
\end{proof}

In a future paper, we will outline how these approximations of $C(H,G)$ as $m \to 0$ extend to histograms $H$ and $G$ that are allowed to have some coordinates equal to $0.$ The proofs and computations of such generic cost approximations are more  complicated than the case when $H,G \in \mathcal{H}^\circ.$ In addition, the differentiability of $C(H,G)$ is far more restricted when $H$ or $G$ lie on the boundary of $\mathcal{H}^\circ$. 

We have now completed the large deviations analysis of daily transitions, including the explicit calculation of the rate function given in the above theorem for the transition kernel.  With this analysis completed, we can now progress by lifting these results to path space in order to complete a large deviations analysis of random histogram trajectories, which will form the foundation for our application given in Section \ref{mostlikely}.






\section{Large Deviations for Evolutionary Trajectories}
\label{LDPathTheory}
Fix any \emph{time horizon} $T \geq 2$. For each population size $N$, the stochastic genetic evolution of the population is described here by the random histogram path $ \mbold{H}= [ H_1\, \ldots  H_T ]$ where $H_n$ is the Markov chain studied above.  Call $\Omega_T \equiv \mathcal{H}^T $ the \emph{path space} of all histogram trajectories $\mbold{H} = [ H_1\, \ldots H_T ]$ with all $H_n \in \mathcal{H}.$  We set the following definitions and notation for $\mbold{H}, \mbold{H}' \in \Omega_T:$ 
\begin{itemize}
\item The \emph{essential minimum} $b(\mbold{H})$ is defined by $b(\mbold{H}) = \min_{n =1 \ldots T}  b(H_n).$ 
\item The \emph{distance} between $\mbold{H}$ and $\mbold{H}'$ is given by $\vnorm{ \mbold{H} - \mbold{H}' }= \max_{n =1 \ldots T}  \vnorm{ H_n - H'_n }$.
\item For any $a > 0,$ let $\Omega_T(a) \subset \Omega_T$ be the compact set of all $\mbold{H} \in \Omega_T$ such that $b(\mbold{H}) \geq a,$
\item A trajectory $\mbold{H} = [H_1 \, H_2 \ldots H_T]$ will be called \emph{N-rational} if $H_n \in \mcal{H}_N$ for all $1\leq n\leq T$. 
\item For each trajectory $\mbold{H} \in \Omega_T$, define the \emph{N-rational ball} $B_N(\mbold{H})$ as the set of all $N$-rational trajectories $\mbold{H}'$ such that $\vnorm{ \mbold{H}'- \mbold{H}} < \frac{2}{3 N}$.  
\end{itemize}
We note that $B_N(\mbold{H})$ is a finite subset of the ball $\Omega_T$ of equal radius with cardinality bounded by $ 2^{g T}$.  

With these notations and definitions in mind, our goal in this section is to establish a large deviations framework for the path space $\Omega_T$.  The analysis and large deviations framework for daily transitions presented in Section \ref{LDcycles} will play a key role.  In fact, many of the results for the one-step cost function for daily transitions will lift to an analagous rate function for $\Omega_T.$ 

\subsection{Large deviations for a single trajectory}
We begin with the following definition that relates the one-step cost to the cost of a random path.
\begin{mydef}
For any $\mbold{H} \in \Omega_T$, we define the \emph{large deviations rate function} $\lambda : \Omega_T \to [ 0, \infty ]$ by
\begin{equation} \label{lambda}
\lambda(\mbold{H}) = \sum_{n=1}^{T-1} C(H_n, H_{n+1}).
\end{equation}
We then define the \emph{large deviations set functional} $\Lambda(F) \in [0, \infty]$ for any $F \subset \Omega_T$ by 
\begin{equation} \label{Lambda}
\Lambda(F) = \inf_{ \mbold{H} \in F} \; \lambda(\mbold{H}).
\end{equation}
\end{mydef}
Analogous to the H\"{o}lder property established for the one-step cost function in Theorem \ref{holderCost}, we will need a similar property for the large deviations rate function.  The following theorem essentially lifts Theorem \ref{holderCost} to the function $\lambda.$
\begin{theorem} \label{holderlambda}
Fix the parameter set $\mcal{P}$ along with positive constants $d$ and $a$. Consider any path $\mbold{H} \in \Omega_T $ such that $ b(\mbold{H}) \geq a$ and $\lambda(\mbold{H}) \leq d$. Let $\mbold{H}' \in \Omega_T $ be any path that satisfies the following for all integers $1 \leq n \leq T:$
\begin{equation} \label{lambda1}
b(\mbold{H}') \geq a, \qquad \vnorm{\mbold{H}' - \mbold{H}} \leq a/F_g, \qquad  spt(H'_n) = spt(H_n). 
\end{equation}
There is a constant $c = c(d, a, \mcal{P})$ such that for all $T ,\mbold{H},$ and $\mbold{H}'$ as above, one has
\begin{equation} \label{lambda2}
\abs{ \lambda(\mbold{H}') - \lambda(\mbold{H})} \leq c T\vnorm{ \mbold{H}' - \mbold{H}}^{1/2}.
\end{equation}
\end{theorem}
\begin{proof}
By Theorem \ref{holderCost}, there is a constant $c_1(d,a,\mcal{P})$ such that 
$$
\abs{ C(H'_n, H'_{n+1}) - C(H_n, H_{n+1}) } \leq c_1 [\vnorm{ H'_n - H_n }^{1/2} + \vnorm{ H'_{n+1} - H_{n+1}}]
$$
for each $ n=1, \ldots, T-1 .$ This implies $\abs{ C(H'_n, H'_{n+1}) - C(H_n, H_{n+1}) } \leq 2 c_1 \vnorm{ \mbold{H}' - \mbold{H}}^{1/2} $. By definition of $\lambda$, we then have
$$
\abs{\lambda(\mbold{H}') - \lambda(\mbold{H})} \leq \sum_{n=1}^{T-1}\, \abs{ C(H'_n, H'_{n+1}) - C(H_n, H_{n+1}) } \leq cT \vnorm{\mbold{H}' - \mbold{H}}^{1/2},
$$
which proves \eqref{lambda2} with $c = 2 c_1$. 
\end{proof}
Now, with the regularity of $\lambda$ in mind, we can now justify calling $\lambda$ the rate function with the next theorem.  This large deviations result for single trajectories will be important in extending to the large deviations result in Theorem \ref{thsetLD} for sets of paths.

\begin{theorem} \label{thtrajLD}
For any path length $T \geq 2$, denote $\mbold{H} = [H_1 \, \ldots H_T]$ as the random trajectory of population histograms. Fix the parameters $\mcal{P}$ and any positive constants $d >0$ and $a >0$. Then, the triple $(d,a,\mcal{P})$ determine positive constants $c$ and $N_0$ such that the following holds. For any N-rational path $\mbold{h}= [h_1 \, h_2\, \dots \,h_T] \in \Omega_T$ such that $\lambda(\mbold{h}) \leq d$ and $b(\mbold{h}) \geq a$, one has  
\begin{equation} \label{trajsingle}
\frac{1}{N} \log P( \mbold{H} = \mbold{h} \; | \; H_1) = - \lambda(\mbold{h}) + o(N)
\end{equation} 
for all $N > N_0$ with $\abs{ o(N) } \leq cT / \sqrt{N}.$
\end{theorem}
\begin{proof}
For all $n \leq T-1$, we must have $C(h_n, h_{n+1})\leq d$ since $\lambda(\mbold{h}) \leq d$. Hence, Theorem \ref{transLD} provides constants $N_0$ and $c$, determined by $(d, a, \mcal{P})$, such that the Markov transition kernel $\mathcal{Q}$ verifies 
$
\frac{1}{N} \log \mathcal{Q}(h_n , h_{n+1}) = - C(h_n , h_{n+1}) + \eta_n(N) 
$ 
with $\abs{ \eta_n(N) } \leq c / \sqrt{N}$ for all $\mbold{h}$ as above, $n \leq T-1$, and $N > N_0.$  The Markov property yields 
$
\log P( \mbold{H} = \mbold{h} \; | \; H_1) = \sum_{n =1}^{T-1} \log \mathcal{Q}(h_n , h_{n+1}). 
$
Since $\lambda(\mbold{h}) = \sum_{n =1}^{T-1} C(h_n , h_{n+1}),$ we obtain
\begin{equation} \label{trajh}
\frac{1}{N} \log P( \mbold{H} = h \; | \; H_1) = - \lambda(\mbold{h}) + o(N) 
\end{equation} 
with $\abs{o(N)} \leq \sum_{n =1}^{T-1} \abs{ \eta_n(N) } \leq cT / \sqrt{N}$ for $N > N_1$.
\end{proof}
This large deviations result for a single random trajectory now sets the stage for a large deviations result for sets of trajectories.

\subsection{Large deviations for sets of trajectories}
 We will show that for large $A$, the probability of observing random paths $\mbold{H}$ of population histograms such that $\lambda(\mbold{H}) > A$ is bounded above by $e^{-c A N}$ for some constant $c>0$.  This will naturally lead to the main result of this section given by Theorem \ref{thsetLD}.  We will first need a definition of particular open neighborhoods similar to the Definition \ref{balls} for histograms. 
\begin{definition} \label{VNBNpaths}
For any $\Gamma \subset \Omega_T$, define the open neighborhood $V_N(\Gamma)$ as the union of all balls $V_N(\mbold{H})$ with radius $\frac{2}{3 N}$ and arbitrary center $\mbold{H} \in \Gamma$.  Denote $B_N(\Gamma)$ as the (finite) set of $N$-rational paths in $V_N(\Gamma)$. Define also $b(\Gamma) = \inf_{\mbold{H} \in \Gamma} b(\mbold{H}).$
\end{definition}

\begin{theorem} \label{thlargecost}
Fix $0 < a < 1$, the parameters $\mcal{P}$, and any initial histogram $H_1 = H$ with $b(H) \geq a$ for the random path $\mbold{H}.$ There is a constant $c = c(a,\mcal{P}) >0$ such that for all $A > c T$ 
\begin{equation} \label{largecost}
P( \lambda(\mbold{H}) > A \; \; \text{and} \;\; b(\mbold{H}) > a \; | \; H_1 = H) \leq e^{- \frac{a}{2 T} N A}
\end{equation}
for all $N > c A/T.$
\end{theorem}
\begin{proof}
Assume the setting and notation of Lemma \ref{integral}.  Consider arbitrary $N$-rational $H,G \in \mcal{H}(a)$ and $r \in K_N(H)$.  Lemma \ref{integral} provides positive $c_0(a,\mcal{P})$ and $N_0(a,\mcal{P})$ such that for $N > N_0,$
\begin{align} 
&spt(J) = spt(I) \quad \text{with} \quad \vnorm{ J - I }  \leq  c_0 / N,    \label{JI1}\\
&\mathcal{Q}(H, G)  =  \sum_{r  \in  E_N}  \exp (- N  [ T(H,r,G) + o(N) ] )   \label{QHG1}
\end{align} 
where $\abs{o(N)} \leq c_0  \log N  / N$ and $T(H,r,G) = mut(H,r) + KL(G,J)$. 
 
Suppose that $C(H,G) > A$ so that $\tau(H,r,G) \geq C(H,G) > A $ for $r \in K(H)$.  From \eqref{mutbound}, we get $c= c(a,\mcal{P})$ such that $mut(H,r) \leq c$ whenever $b(H) \geq a$ and $r \in K(H)$.  We then have $KL(G,I) = \tau(H,r,G) - mut(H,r) > A - c$.  This yields $\beta(G,I) = \max_{k \in spt(G)}  {1/I(k)} \geq KL(G,I) > A - c$ by \eqref{betabound}.  Therefore, for some $j \in spt(G),$ one has $I(j) < 1 / (A-c)$ so that $J(j) < 1 / (A -c) + c_0 / N$ by \eqref{JI1}. This yields $1/J(j) > \frac{3}{4} (A - c)$ provided $N >  N_0 + 8 c_0 A $.
A fortiori, we get $\beta(G, J) = \max_{k \in spt(G)}  {1/J(k)} > \frac{3}{4} (A - c)$.  For $N > N_0 + 8 c_0 A $ and $r \in E_N$,
\begin{equation} \label{L(A)}
T(H,r,G) \geq KL(G,J) > b(G) \beta(G,J) - \log{g} > \frac{3 a}{4} (A - c) - \log(g) = \frac{3 a A}{4} - c_1 
\end{equation} 
by \eqref{betabound} with  $c_1 = \log(g) + 3 a c /4$. Since $card(E_N) \leq [ (N+1) F_g ]^{g^2}$, combining \eqref{L(A)} and \eqref{QHG1} yield 
\begin{equation} \label{Q1}
\mathcal{Q}(H, G)  \leq [ (N+1) F_g ]^{g^2} \exp \left( - N \left[ \frac{3 a A}{4} - c_1 \right] + c_0  \log N\right) \leq \exp \left( N \left[- \frac{3 a A}{4} + c_2 \right]\right)   
\end{equation} 
for $N > N_0 + 8 c_0 A$ with $ c_2  =  c_1 + c_0 +g^2 (\log F_g +1)$.

Let $\Gamma(A) = \vset{ \mbold{H}\in \Omega_T \; | \; \lambda(\mbold{H}) > A;\, b(\mbold{H})> a}$. Let $\Gamma_N(A)$ be the set of $N$-rational $\mbold{H} \in \Gamma(A)$.  Then, $P( \mbold{H} \in \Gamma_N(A) \; | \; H_1 =H) = P( \mbold{H} \in \Gamma(A) \; | \; H_1=H)$ since all random paths $\mbold{H}$ are $N$- rational. For $\mbold{H} \in \Gamma_N(A)$, the relation $\lambda(\mbold{H}) > A$ provides at least one time step $\nu = \nu(\mbold{H})  \leq T$ such that $C(H_\nu, H_{\nu+1}) > A/T$. Apply \eqref{Q1} to $(H_\nu, H_{\nu+1})$ to get
$$
P(\mbold{H} = \mbold{h} \; | \; H_1= H) \leq P(H_{\nu+1} = h_{\nu+1} \; | \; H_\nu = h_\nu) \leq \exp \left( N \left[- \frac{3 a}{4 T} A + c_2\right]\right) 
$$
for $N > N_0 + 8 c_0 A / T$ and $\mbold{H} \in \Gamma_N(A).$  Since $card(\Gamma_N(A)) \leq (N+1)^{T g}$, this yields
\begin{align*}
P(\mbold{H} \in \Gamma(A) \; | \; H_1=H) &= P( \mbold{H} \in \Gamma_N(A)\; | \; H_1=H) \\
&\leq (N+1)^{T g} \exp\left(N\left[- \frac{3 a}{4 T} A + c_2\right]\right) \\
&\leq \exp\left( N\left[- \frac{3 a}{4 T} A + c_3\right]\right)
\end{align*}
for $N > N_0 + 8 c_0 A / T$ with $c_3 = c_2 + T g$. Impose now $A/T  > 4 c_3 $ to get 
$$ 
P( \mbold{H} \in \Gamma(A) \; | \; H_1=H)  \leq \exp\left(- \frac{a}{2T}  N A \right) 
$$
for $N > N_0 + 8 c_0 A/T.$ Set $c = 4 c_3 + N_0 /8 c_0 +16 c_0$, so that the simpler constraints $A/T > c$ and $N > c A/T$ force $A/T > c_3$ and $N > N_0 + 8 c_0 A/T$. This completes the proof.
\end{proof}
\begin{theorem} \label{thsetLD}
Let $\Lambda$ be the large deviations set functional of our Markov chain of population histograms defined by \eqref{lambda} and \eqref{Lambda}. Fix the parameters $\mcal{P}$, the path length $T,$ and $0<a < 1$.  Let $\Gamma \subset \Omega_T$ be the set of histograms trajectories with $b(\Gamma) \geq a$ and finite $\Lambda(\Gamma) = L > 0$ where $b(\Gamma)$ is given in Definition \ref{VNBNpaths}.  Let $V_N(\Gamma)$ be as in Definition \ref{VNBNpaths}. Denote $\mbold{H}$ as the generic random histogram trajectory starting at some fixed $H_1$. Assume that all paths in $\Gamma$ also start at $H_1$.  We then have that $(T,a, \mcal{P},L)$ determine positive constants $c$ and $N_0$ such that for all $N > N_0$,the uniform large deviations result 
\begin{equation} \label{setLD}
\frac{1}{N} \log P( \mbold{H} \in V_N(\Gamma) \; | \;  H_1 ) = - \Lambda(\Gamma) + o(N)
\end{equation}
holds with $\abs{o(N)} \leq c / \sqrt{N}.$ This yields the asymptotic large deviations limit 
$$
\lim_{N \to \infty}  \frac{1}{N} \log P( \mbold{H} \in V_N(\Gamma) \; | \; H_1 ) = - \Lambda(\Gamma).
$$
\end{theorem}
\begin{proof}
Let $\Gamma$ be any subset of $\Omega_T$ with $b(\Gamma) \geq a$ and  finite $\Lambda(\Gamma) = L > 0$. For any $A>0$, let $E(A)$ be the set of all $N$-rational paths $\mbold{h} \in \Omega_T$ such that $\lambda(\mbold{h}) > A $ and $b(\mbold{h})> a/2$. Theorem \ref{thlargecost} provides $c = c(a,\mcal{P})$ such that  
$
P( \mbold{H} \in E(A) \; | \; H_1) \leq e^{- \frac{a}{2 T} N A } 
$
for all $A > c T$ and $N > c A/T.$ Set $d = T (c + 4 L / a) + L$. Set $A=d$, which forces $A >cT$ and $\frac{a}{2 T} A > 2 L$ so that 
\begin{equation} \label{PED}
P( \omega \in E(d) \; | \; H_1) \leq e^{- 2 N L } 
\end{equation} 
for all $N > N_1 = c d/T.$ Let $W_N(\Gamma) \subset B_N(\Gamma) \subset V_N(\Gamma)$ be the set of $N$-rational paths $\mbold{h} \in V_N(\Gamma)$ such that $\lambda(\mbold{h}) \leq d$.  Then, $B_N(\Gamma) \subset W_N(\Gamma) \cup E(d)$ so that \eqref{PED} yields
\begin{equation} \label{set0}
P( \mbold{H} \in B_N(\Gamma) \; | \; H_1) \leq  P( \mbold{H} \in W_N(\Gamma)  \; | \; H_1) +  e^{-2 L N}
\end{equation} 
for $N > N_1.$ Due to Theorem \ref{holderlambda}, there is a constant $c_1 > 0$ determined by $(d,a,\mcal{P})$, and thus by $(T, L, a, \mcal{P}),$ such that for any paths $\mbold{h}, \mbold{h}' \in \Omega_T,$ the inequality
\begin{equation} \label{lambda11}
\abs{\lambda(\mbold{h}') - \lambda(\mbold{h})} \leq c_1 \vnorm{\mbold{h}' - \mbold{h}}^{1/2}
\end{equation}
must hold whenever the following holds for all integers $1 \leq n \leq T:$
\begin{align} 
&\min\vset{\lambda(\mbold{h}), \lambda(\mbold{h}')} \leq d;    \label{lambda1a} \\
b(\mbold{h}') \geq a&; \qquad \vnorm{\mbold{h}' - \mbold{h}} \leq a/F_g; \qquad spt(h'_n) = spt(h_n) \label{lambda2a}.
\end{align}
Since $b(\Gamma) \geq a$, Lemma \ref{basicBN} implies $b(V_N(\Gamma)) \geq a/2$ so that  all $\mbold{h} \in W_N(\Gamma)$ verify $\lambda(\mbold{h}) \leq d$ and $b(\mbold{h}) > a/2$. Therefore, Theorem \ref{thtrajLD} provides $c_2$ and $N_2 > N_1$ determined by $(d, a, \mcal{P}),$ and thus by $(T, L, a, \mcal{P})$, such that
\begin{equation} \label{set1}
\frac{1}{N} \log P( \mbold{H} = \mbold{h} \; | \; H_1) = - \lambda(\mbold{h}) + o_1(N)
\end{equation}
for all $N > N_2$ and $\mbold{h} \in W_N(\Gamma)$ with $\abs{o_1(N)} \leq T c_2 / \sqrt{N}.$ For each $\mbold{h} \in W_N(\Gamma) \subset V_N(\Gamma),$ one can select a path $\hat{\mbold{h}} \in \Gamma$ such that $\vnorm{\hat{\mbold{h}}- \mbold{h}} \leq 1/N$. By Lemma \ref{basicBN}, the paths $\mbold{h}$ and $\mbold{h}' =\hat{\mbold{h}}$ verify both \eqref{lambda1a} and \eqref{lambda2a} so that \eqref{lambda11} applies and yields 
$
\abs{ \lambda(\hat{\mbold{h}}) - \lambda(\mbold{h})} \leq c_1 \vnorm{ \hat{\mbold{h}} - \mbold{h}}\leq c_1 / \sqrt{N}
$ 
provided $N > F_g / a.$  This implies $ - \lambda(\mbold{h}) +  o_1(N) = - \lambda(\hat{\mbold{h}}) + o_2(N)$ with $\abs{ o_2(N)} \leq (c_1 + T c_2) / \sqrt{N}$ for $N >N_3 =N_2 + F_g /a$.  For $N>N_3$ and $\mbold{h} \in W_N(\Gamma)$, this yields
$$
P( \mbold{H} = \mbold{h} \; | \; H_1 ) = \exp (- N \lambda(\hat{\mbold{h}}) + N o_2(N)) \leq \exp (- N L + c_3 \sqrt{N})
$$
with $c_3= c_1 + T c_2$ due to \eqref{set1} and $\lambda(\hat{\mbold{h}}) \geq \Lambda(\Gamma) = L.$ For $N > N_3$, sum over $\mbold{h} \in W_N(\Gamma)$ to get 
$$
P( \mbold{H} \in W_N(\Gamma) \; | \; H_1 ) \leq \exp [- N L + c_3 \sqrt{N} + gT  \log(N+1)] \leq \exp (- N L + (c_3 + gT) \sqrt{N})
$$
since $card(W_N(\Gamma)) \leq (N+1)^{gT}.$  Due to \eqref{set0}, this gives 
\begin{align}\label{PBNup}
P( \mbold{H} \in B_N(\Gamma) \; | \; H_1 ) &\leq \exp\left[- N L + (c_3 + T g) \sqrt{N}\right] + e^{-2 N L} \nonumber \\
&\leq 2 \exp\left[- N L + (c_3 + gT) \sqrt{N}\right]
\end{align}
for $N > N_3.$ One has $P( \mbold{H} \in V_N(\Gamma) \; | \; H_1 ) = P( \mbold{H} \in B_N(\Gamma) \; | \; H_1 )$ since the random paths $\mbold{H}$ are always $N$-rational. Hence for $N > N_3$, \eqref{PBNup} yields the upper bound,
\begin{equation} \label{logPVNup}
\frac{1}{N} \log P( \mbold{H} \in V_N(\Gamma) \; | \; H_1 ) \leq - L + c_4 / \sqrt{N}
\end{equation}
with  $c_4 =c_3 + gT +1.$ Set $c5 = \left[\frac{1+c_1}{T (c + 4 L/a)}\right]^2$and $N_4 = N_3 + c_5$. This yields $L +(1 + c_1) / \sqrt{N} \leq d = L +T (c + 4 L / a)$ for $N > N_4.$  For each $N,$ select $\mbold{h}^* = \mbold{h}^*(N) \in \Gamma$ such that $L \leq \lambda(\mbold{h}^*) \leq L + 1 / \sqrt{N}$. For $N > N_4$, one has $\lambda(\mbold{h}^*) \leq d $, and \eqref{lambda11} applies to the pair of trajectories $\mbold{h}^*$  and $\mbold{h}'$ for any $N$-rational $\mbold{h}' \in V_N(\mbold{h}^*)$ to give
$
\abs{\lambda(\mbold{h}') - \lambda(\mbold{h}^*) }  \leq c_1 / \sqrt{N}.
$
Hence, $\lambda(\mbold{h}') \leq L + (1+c_1) / \sqrt{N} \leq d$ which shows that $\mbold{h}' \in W_N(\Gamma)$.  Apply \eqref{set1} to $\mbold{h}' \in W_N(\Gamma)$ to obtain
\begin{equation}\label{lowPh'}
\frac{1}{N} \log P( \mbold{H} = \mbold{h}' \; | \; H_1 ) = - \lambda(\mbold{h}') + o_2(N) \geq - L - (1 + c_1 +T c_2) / \sqrt{N}
\end{equation}
for $N> N_4.$ Since $\mbold{h}' \in W_N(\Gamma) \subset V_N(\Gamma),$ one has 
$$
P( \mbold{H} \in V_N(\Gamma) \; | \; H_1 ) \geq P( \mbold{H} \in W_N(\Gamma) \; | \; H_1 ) \geq P( \mbold{H} = \mbold{h}' \; | \; H_1 ). 
$$
Combining this with \eqref{lowPh'} yields the lower bound
\begin{equation} \label{logPVNlow}
\frac{1}{N} \log P( \mbold{H} \in V_N(\Gamma) \; | \; H_1 ) \geq - L - T (1 + c_1 + T c_2) / \sqrt{N}
\end{equation}
for $N>N_4.$ Combining \eqref{logPVNup} and \eqref{logPVNlow} concludes the proof.
\end{proof}
Theorem \ref{thsetLD} naturally gives a notion of the likelihood of observing the random histogram trajectory $\mbold{H}$ in an generic set of paths $\Gamma \in \Omega_T$ quantified by its set functional value $\Lambda(\Gamma).$  Smaller values of $\Lambda(\Gamma)$ would correspond to a higher likelihood of observance. Consequently, Theorem \ref{thsetLD} is vital in applying large deviations to specific paths where we fix the beginning and ending histogram in order to obtain information on the likelihood of an initial population evolving over time to a particular state of interest.  In addition, the actual calculation of $\Lambda(\Gamma)$ involves a minimization of the rate function $\lambda(\mbold{h})$ over paths $\mbold{h}\in \Gamma.$  Thus, if $\mbold{h}^* \in \Gamma$ is the optimal path, we would expect the occurrence of the event $(\mbold{H} \in \Gamma)$ to be roughly quantified by $\lambda(\mbold{h}^*)$ and for the random histogram trajectory $\mbold{H}$ to follow closely to the minimizing path $\mbold{h}^*,$ giving a notion of the most likely way in which the system achieves this event.  We explicitly quantify this line of thinking in Section \ref{mostlikely}.  

For now, if $\Gamma$ happens to contain a path that has zero cost, then the minimization is trivial, and we would find the occurrence of the event $(\mbold{H} \in \Gamma)$ to be a highly likely event.  This motivates us to quantify these zero-cost paths.

\subsection{Mean evolution and zero-cost trajectories}\label{zerotraj}
Introduce the histogram valued  function $H \to \zeta(H)$ defined for all $H \in \mathcal{H}$ and genotypes $j$ by 
\begin{equation} \label{zeta}
\zeta_j(H)  = \frac{1}{\langle F, H \rangle} \left( F_j H(j) - m \sum_k  Q_{j,k} F_j H(j) + m \sum_k  Q_{k,j} F_k H(k) \right).
\end{equation}
Once an initial histogram is set, the histogram above recursively categorizes all zero-cost trajectories, which is detailed in the theorem below.
\begin{theorem} \label{thzerotraj}
Fix a path length $T$. A histogram path $\mbold{h}= [ h_1\, \dots, h_T] \in \Omega_T$  satisfies $\lambda(\mbold{h}) = 0$ if and only if  $h_{n+1} = \zeta(h_n)$ for $1 \leq n \leq T-1$, where the function  $\zeta$ is defined by \eqref{zeta}. Consequently, a zero cost path $\mbold{h}$ is uniquely determined by its starting point $h_1$.  Fix $h_1 \in \mathcal{H},$ which uniquely determines the zero cost path $\mbold{h}$ starting at $h_1$. Let $V_N(\mbold{h})$ be the ball of center $\mbold{h}$ and radius $\frac{2}{3 N}$  in $\Omega_T$. Then, the initial point $h_1$ and the parameters $\mcal{P}$ determine constants $c$ and $N_0$ such that the random path $\mbold{H} \in \Omega_T$ verifies
\begin{equation}  \label{V(h*)}
1 \geq P( \mbold{H} \in V_N(\mbold{h})  \; | \; H_1 = h_1 ) \geq 1 - e^{- c  \sqrt{N}} 
\end{equation}
for all $N > N_0.$
\end{theorem}
\begin{proof}
Given any two histograms $H,G\in \mcal{H}$ with finite cost $C(H,G)$,  Lemma \ref{finitecost} proves the existence of a matrix $r \in \overline{K(H)}$ such that $C(H,G) = \tau(H, r, G)$. Hence, $C(H,G) = 0$ if and only if $\tau(H,r,G) = 0,$ which is equivalent to $mut(r,H) = KL(G, \Psi(H,r) ) = 0$.  From \eqref{zeromutcramer}, one has $mut(r,H) = 0$ iff 
\begin{equation} \label{r.optimal}
r_{j,k} = m Q_{j,k} F_j H(j) = M_{j,k} F_j H(j) 
\end{equation}
for all $j,k.$ By definition of the Kullback-Leibler divergence, one has  $KL(G, \Psi(H,r) ) = 0$ if and only if $\Psi(H,r) = G$. Combine this relation with \eqref{psi} and \eqref{r.optimal} to conclude that $C(H,G) = 0$ if and only if one has
$$
G(j) = \frac{1}{\langle F, H \rangle}\left( F_j H(j) - \sum_k  M_{j,k} F_j H(j) + \sum_k  M_{k,j} F_k H(k) \right).
$$
for all $j.$ Hence $C(H,G) = 0$ if and only if $G= \zeta(H)$. Now, for any path $\mbold{h} \in \Omega_T$, the relation $\lambda(\mbold{h}) = 0$ holds if and only if $C(h_n, h_{n+1} ) = 0$ for all $1 \leq n \leq T-1$, which is equivalent to $h_{n+1}= \zeta(h_n)$ for all $1 \leq n \leq T-1.$  Finally, the bound in \eqref{V(h*)} is an immediate consequence of \eqref{setLD} applied to the set $F = \vset{\mbold{h} \;|\; \lambda(\mbold{h})= 0}.$
\end{proof}
We will show that as $n \to \infty,$ any infinite zero-cost path $\mbold{h} = [h_n]_{n=1}^{\infty}$ achieves near fixation of some explicitly determined genotype.  Before doing so, we make a quick definition.
\begin{definition}\label{reachgen}
For any non-empty set $S$ of genotypes, define the set $R(S)$ of \emph{genotypes reachable from S} as the set of all genotypes such that there is a genotype sequence $k_1, \ldots, k_T$ of arbitrary length $T$ such that $k_1 \in S, k_T = j,$ and all $M_{k_t, k_{t+1}} > 0.$
\end{definition}
\begin{theorem} \label{hinfty} 
Let $\mbold{h}$ be any zero-cost histogram path of infinite length starting at $h_1 \in \mathcal{H}$. Denote $S_n =spt(h_n)$. Then, the set $S_n$ increases with $n$ and thus stabilizes to $S_n= S_{ter}$ for all $n > n_0$ for some finite $n_0$. Moreover, $S_{ter} = R(S_1)$ where $R(S_1)$ is the set of all genotypes reachable from $S_1$ as in Definition \ref{reachgen}. Let $s \in S_{ter}$ be the fittest genotype within $S_{ter}$. Then, the initial point $h_1$ and the parameters $\mcal{P}$ determine $m_0 >0$ such that for all mutation rates $m \leq m_0$,    
\begin{equation}\label{fixation}
\lim_{n \to \infty}  h_n = h_{ter} 
\end{equation}
where $h_{ter}$ is the unique solution of $ h_{ter} = \zeta(h_{ter})$, has support $spt(h_{ter}) = S_{ter},$ and is a $C^{\infty}$ function of $m$ with first-order expansion
\begin{align} \label{fixedpoint}
h_{ter}(i) &\simeq m  Q_{s, i} / (F_s - F_i); \qquad  i \in (S_{ter} - s), \\
h_{ter}(s) &\simeq 1- m\sum_{ i \in (S_{ter} - s)} Q_{s, i} / (F_s - F_i).
\end{align}
\end{theorem}
\begin{proof}
Define  $c > 1$ and $0 < \eta <1 $ by $c = 1 + \max_j \left\{ \sum_k Q_{j,k} \right\}$ and $1- \eta = \max_{1\leq k\leq g-1}  F_k/ F_{k+1}.$   Fix temporarily a constant $0 < A <1/2$ to be selected later on. Impose $m c \leq A $ so that all $m_j = m  \sum_k Q_{k,j} \leq A.$  For $H \in \mathcal{H},$ set $U_j(H) = F_j H(j) (1 - m_j) + m \sum_k Q_{k,j}F_k H(k)$. Setting $G(j)= \zeta_j(H) = U_j(H) / <F, H>,$ we have that $G(j) = 0$ if and only if $H(j) = 0$ and $H(k) = 0$ for all $k$ such that $Q_{k,j}>0.$  This implies $spt(H) \subset spt(\zeta(H)).$ Consequently, since $h_{n+1} = \zeta(h_n)$, the set $S_n =spt(h_n)$ must increase with $n$. Therefore, there exists a fixed set $S$ of genotypes and an $n_0$ such that $S_n = S$ for $n \geq n_0$. An easy recurrence based on the conditions above for $G(j) = \zeta_j(H) = 0$ shows that $S$ is the set of all $j$ reachable by some finite sequence $k_1, k_2, \ldots, k_T = j$ with $k_1 \in spt(h_1)$ and all $M_{k_t, k_{t+1}} >0$.
 
Now, call a genotype $q$ \emph{dominant} in  $H$ if $H(q) = \max_k \; H(k) $. Select any dominant $q= q(H)$ in $H$. Then, we have $\langle F, H \rangle  \leq  g F_g H(q);$ $1/g \leq H(q);$ and
\begin{equation}\label{boundLj}
(1 - A) F_j H(j)  \leq  U_j (H)  \leq  F_j H(j) + A \langle F, H\rangle  \leq  F_j H(j)+ g A F_g H(q).
\end{equation}
for all $j$.  Since $G(j) / G(q) = U_j(H) / U_q(H)$, this implies
\begin{equation}\label{GjGi}
\left(\frac{(1 - A) F_j} {( F_q + g A F_g )}\right) \left(\frac{H(j)}{H(q)}\right)  \leq  \frac{G(j)}{G(q)} \leq \frac{F_j H(j) + g A F_g H(q)}{F_q H(q) (1 - A)}.
\end{equation}
For  $j < q$, one has $H(j) \leq H(q)$ and $F_j \leq (1 - \eta) F_q$ by definition of $c$ and $\eta.$ Then, \eqref{GjGi} yields 
\begin{equation}\label{upGG}
\frac{G(j)}{G(q)} \leq  \left(\frac{1-\eta}{1-A}\right)\frac{H(j)}{H(q)} + \frac{g A  F_g}{(1-A) F_q} \leq  \frac{1-\eta + A F_g/F_1}{1-A}.  
\end{equation} 
Impose $A < A_1 =\frac{\eta}{1+F_g/F_1}$ to get  
$
G(j)/G(q) \leq 1-\eta/2 < 1 
$
for all $j < q.$ Hence, the dominant genotype $q(G)$ in $G=\zeta(H)$ verifies $q(G)  \geq q(H)$.  Consequently, the dominant genotype $q_n = q(h_n)$ of $h_n$  verifies $q_{n+1} \geq q_n$ for all $n$ so that there exists a finite $n_1 \geq n_0$ and a genotype $q^*$ such that $q_n = q^*$ for all $n > n_1$. From \eqref{boundLj}, we get $h_n(q_n) \geq 1/g$ and hence $q^* \in S$.  Let $s = \max(S)$ be the fittest genotype in $S$ so that $q^* \leq s$.  Assume there is a $j \in S$ verifying $ q^* < j \leq s$. We will proceed by contradiction to show that such a $j$ cannot exist. For $n > n_1$ one has $h_n(j)  \leq  h_n(q^*)$ and, by \eqref{GjGi},
\begin{equation} \label{recur}
\theta \left(\frac{h_n(j)}{h_n(q^*)}\right) \leq \frac{h_{n+1}(j)}{h_{n+1}(q^*)}  \leq  1
\end{equation}
with $\theta  =  (1 - A) F_j / ( F_{q^*} + g A F_g ) $.  Since $j > q^*,$ we get $F_j  \geq  F_{q^*}/ (1-\eta)$ by definition of $\eta$ so that
$$
\theta  \geq  \frac{1 - A}{(1-\eta) ( 1 + g A  F_g/F_{q^*})}  \geq  \frac{1 - A}{(1-\eta) ( 1 + A g F_g/F_1 )}.
$$
Impose $A < A_2 =\frac{\eta^2  F_1}{2 g F_g}$ to  force $\theta > 1 +\eta$. By recurrence, \eqref{recur} implies for $s \geq j > q*$ and $n' > n > n_1$ that $(1+\eta)^{n' - n}\frac{h_n(j)}{h_n(q^*)} \leq  1.$  Fix $n >n_1$. Since $j \in S$, one has $h_n(j) >0.$ Letting $n' \to \infty$ yields a contradiction. Hence there is no $j \in S$ with $ q^* < j \leq s$, so that $q^* = s$.

For $i \in S$ with $i < s$ and $n > n_1$, the terms $y_n(i) = h_{n}(i) / h_n(s)$ verify $y_{n+1}(i) \leq  \mu y_n(i) + \nu$  by \eqref{upGG} and $q_n= s$ with $\mu = 1- \eta / 2$ and $\nu = 2 g A  F_g / F_1 .$  Iterating this inequality gives $y_n(i) \; \leq \; \mu^{n-n_1} y_{n_1}(i) + \frac{\nu}{1-\mu}$.
Select $n_2 = n_2(A) > n_1$ to force $(1-\eta / 2)^{n-n_1} y_{n_1}(i) < A$ for all $n > n_2$ and $i \in (S - s).$ Then, $y_n(i)  \leq c_1 A$, with $c_1 = 1 + \frac{4 g F_g }{\eta  F_1}$.  This yields $h_n(i)  \leq  c_1 \; A$ for all $n > n_ 2(A)$ and $i \in (S - s)$.  Hence $h_n(s)  \geq  1 -  c_1 (g-1) A$ for $n > n_2(A)$, since $spt(h_n) = S$.  We now fix  $A = A_3 = \min\vset{ A_2, \frac{1}{2 g c_1} }$ and set $n_3 = n_2(A_3)$ to get the fixed lower bound $h_n(s) \geq 1/2$ for all $n> n_3.$  For $n > n_3$ let $z_n$ be the vector of all the $z_n(i) = y_n(i)$ with $i \in (S - s)$. Since $spt(h_n) = S$, the nonzero $U_i(h_n)$ \emph{only depend on} $z_n$ and can  be denoted $U_i(z_n)$.  We have then $z_{n+1} = f(z_n)$ with $f_i(z)  = U_i(z) / U_s(z)$.  Set $p = card(S) -1.$ The rational fraction $f(z)$ is well-defined on the set $\Gamma$ of  all $z \in[0, 1]^p$ such that $U_s(z) \geq 1/2$ for $0 \leq m \leq A_3/c$.  For the limit case $m=0$ of no mutations, the function $f(z)$ has the form $\hat{f}_i (z) = \frac{F_i}{F_s}z_i$ for all $i \in (S - s)$. The Jacobian matrix $D_z \hat{f}$ obviously verifies
\begin{equation} \label{Dzf}
\vnorm{ D_z \hat{f}} \leq \max_{i \in (S - s)} \frac{F_i}{F_s}  \leq 1 - \eta
\end{equation}
for all $z \in [0,1]^p.$  For $z \in \Gamma$,  the numerator and denominator of each rational fraction $f_i(z)$ are separately affine in $m \leq A_3 / c$ and $z$, with denominator bounded below by $1/2$ and uniformly bounded coefficients. Elementary algebraic computations then prove that as $m \to 0$, the Jacobian matrix $D_z f (z)$ tends to $D_z \hat{f}(z)$ \emph{uniformly} over all $z \in \Gamma$.  Due to \eqref{Dzf}, this provides a constant $0 < A_4 < A_3/c$ such that $\vnorm{ D_z f } \leq (1 - \eta/2 )$ for all $m < A_4$ and all $z \in \Gamma$. This yields $\vnorm{ z_{n+1} - z_n } \leq (1-\eta/2) \vnorm{ z_n - z_{n-1} }$ for $m < A_4$ and $n > 1+n_3$. This contraction property classically shows that $z^*= \lim_{n \to \infty} z_n$ exists with $z^* \in \Gamma.$ Furthermore, $z^*$ is the \emph{unique} solution of $z^* - f(z^*) = 0$. Note that $z^*$ is an implicit  function of $m$ for $0 \leq  m < A_4$, with $z^*(0) = 0$.  Due to \eqref{Dzf}, the Jacobian $Id - D_z f$ of $z -f(z)$ is invertible for $m =0$. Since $f(z)$ is a $C^{\infty}$ function of $(m, z) \in [0, A_4) \times\Gamma$, the implicit function theorem applies to $z - f(z) =0$ and proves the existence of a constant $0 < A_5 < A_4$ such that $z^*$ is a $C^{\infty}$ function of $m \in [0, A_5).$  Since $( 1- h_n(s) ) / h_n(s) = \sum_{i \in (S-s)} z_n(i)$ and $h_n(i) = h_n(s) z_n(i)$ for $i \in (S-s)$, we see that $H= \lim_{n \to \infty} h_n$ exists with $spt(H) \subset S$ and verifies $H(s) = 1 / [ 1 + \sum_{i \in (S-s)}  z^*(i)]$ and $H(i) = H(s) z^*(i).$  Hence, $H$ is a $C^{\infty}$ function of $m \in [0, A_5)$ and $spt(H) = S$. For $m=0$, the solution of $z - \hat{f}(z) = 0$ is clearly $z^* = 0$, and the associated histogram $H$ verifies $H(s) = 1$ and $spt(H) = \vset{s }.$  The first-order Taylor expansion of $H$ as a function of $m$ is of the form $H(i) \simeq \;  m v_i$ for $i \in (S - s)$ and $H(s) \simeq \; 1- u m$. Substitute this into $H = \zeta(H)$ to get  $v_i = Q_{s, i} / (F_s - F_i)$ and $u = \sum_i v_i$. This concludes the proof.
\end{proof}

With our generic large deviations framework rigorously presented, we will now focus on the main application of interest: what is the most likely evolutionary trajectory linking a fixed initial histogram to a fixed terminal histogram?  In addition, what is the most likely time it takes for such a transition to occur?  









\section{Most Likely Evolution from Initial to Terminal Histograms} \label{mostlikely}
\input{sec5_LDApp.tex}

\section{Geodesic Computation by Reverse Shooting}
\input{sec_shooting.tex}


\section{Numerical Simulations}
\label{s:sims}

To illustrate our approach, 
we perform several sets of numerical simulations to verify conjectures presented in the previous section.
We consider various values of parameters which correspond to different evolutionary behavior of biological systems with $g=3$ genotypes. 
In particular, we consider a prototype example with three genotypes and the following parameter set $\mcal{P}$
\begin{equation}
\label{numpar}
F = (200, 200^{1.08}, 200^{1.12}), \quad N = 10^6,  
\quad H = (0.99, 0.005, 0.005), \quad \epsilon= 10^{-4}, \quad \delta = 50/N,
\end{equation}
where $F$ stores the fitnesses for the three genotypes, $N$ is the total population size,  
$H$ is the initial starting histogram of genotypes, $\epsilon$ is the discretization parameter used to generate the set of penultimate points $PEN$, and
$\delta$ is the boundary parameter that terminates trajectories if a histogram coordinate falls below this threshold. We also verified simulations presented here using the discretization parameter $d=10^{-3},$
yielding similar quantitative results.
The mutation matrix is given by
\[
M = mQ \quad \text{with} \,\, m=10^{-6} \,\, \text{and} \,\,
Q = \left(
\begin{array}{ccc}
0 & 1/2 & 1/2 \\
0 & 0 & 1 \\
0 & 0 & 0
\end{array}
\right).
\]
Entries of the matrix $M_{jk}$ are the mutation rates from
genotype $j$ to genotype $k$ and thus 
we only consider non-reversible advantageous mutations $Gen 1 \to Gen 2$, $Gen 1 \to Gen 3$, and $Gen 2 \to Gen 3$.



Recall that calculations of most likely trajectories and corresponding
costs are valid in the space of interior histograms, which is defined by all histograms
$w$ such that $w(j) \ge \delta$ for all $j = 1, 2, 3$.



Since genotype 3 has the highest fitness, the mean (zero-cost) trajectory corresponds to the fast emergence and growth of this genotype. The mean trajectory is depicted in the upper left part of Figure \ref{fig1}.
Recall, that for any histogram $w$, the frequency of the third genotype is computed as $w(3) = 1 - w(2) - w(1)$,
where $w(1)$ and $w(2)$ are frequencies of the first and second genotypes, respectively. 
Note that the genotype 2 never
exceeds $0.1$ and thus the population is quickly dominated by genotype 3.
\begin{figure}[t]
\centerline{\includegraphics[scale = 0.5]{fig1a.pdf}
\includegraphics[scale = 0.5]{fig1b.pdf}}
\centerline{\includegraphics[scale = 0.5]{fig1c.pdf}
\includegraphics[scale = 0.5]{fig1d.pdf}}
\caption{Evolutionary trajectory for the population starting at the histogram $H=[0.99, 0.005, 0.005]$ and three different target
histograms.
Upper left subplot - mean evolutionary trajectory.
Upper right subplot - most likely trajectory for the target histogram $G_1=[0.35, 0.25, 0.4]$,
lower left subplot - most likely trajectory for the target histogram $G_2=[0.35, 0.35, 0.3]$,
lower right subplot - most likely trajectory for the target histogram $G_3=[0.35, 0.5, 0.15]$.
Note the different in the vertical scale (genotype 2) between the upper-left subplot (mean trajectory) and
other sub-plots.}
\label{fig1}
\end{figure}



\begin{figure}[t]
\centerline{\includegraphics[scale = 0.5]{fig2a.pdf}
\includegraphics[scale = 0.5]{fig2b.pdf}}
\caption{Left sub-plot - Cost of the optimal trajectory vs $w(2)$ in the target histogram $G$.
Right sub-plot - number of steps in the optimal trajectory vs $w(2)$ in the target histogram $G$ 
for target histograms in \eqref{gtar}.}
\label{fig1a}
\end{figure}



\textbf{Simulations with varying target histogram $G$.}
In the first set of simulations, we consider target histograms
\begin{equation}
\label{gtar}
G = (0.35, w(2), 0.65 - w(2)), \text{~~with~} w(2) = 0.2,\ldots,0.5.
\end{equation}
Since histograms with relatively large genotype 2 populations ($w(2) \ge 0.2$) are far from the 
mean trajectory, populations with such frequencies are rare events.
In particular, target histograms in \eqref{gtar} correspond to realistic rare events where a genotype with a lower fitness emerges to 
some significant values before the genotype with the highest fitness dominates the whole population. 
For instance, 
the target histogram $G_1 = [0.35, 0.25, 0.4]$ corresponds to the event when the frequency of the genotype 2
is approximately 3 times the normal size (compared to mean trajectory) before the whole population is dominated by genotype 3. 
The target histogram $G_2 = [0.35, 0.35, 0.3]$  corresponds to the event when the
frequencies of all three genotypes are approximately equal, and the target histogram $G_3=[0.35, 0.5, 0.15]$
corresponds to the population with genotype 2 cells dominating. Most likely trajectories for $G_1$, $G_2$, and $G_3$ are 
depicted in top-right, bottom-left, and bottom-right parts of Figure \ref{fig1}, respectively. 

We also computed other most likely trajectories for 
target histograms in \eqref{gtar} using the step $\Delta w(2) = 0.01$. 
Since the frequency of the genotype 1 is fixed in \eqref{gtar}, this set of simulations corresponds to a vertical line $w(1) = 0.35$. Optimal trajectories for these target histograms are
conceptually similar (with a different cost and number of steps) to the most likely trajectories for 
target histograms $G_1$, $G_2$, and $G_3$. All most likely trajectories for target histograms in \eqref{gtar} involve very small shifts
in the genetic composition of the population during early stages (i.e. $H \to w_1 \to w_2 \to \ldots$). However, these small shifts are rather costly, and early steps in the most likely trajectory comprise a very large percentage of the total cost for the most likely trajectory.
These results are presented in detail later in this section and Appendix \ref{app2}.

Figure \ref{fig1a} depicts the cost and the number of steps in most likely trajectories as a function of $w(2)$ in target histograms \eqref{gtar}. Intuitively, events should become more rare as the frequency of the genotype 2 increases (i.e., target histogram moving away from the mean trajectory). This is confirmed in our simulations. In particular, the left sub-plot of Figure \ref{fig1a} demonstrates that the cost function is a monotonically increasing function of the frequency $w(2)$ in the target histogram.
Recall that the probability of each evolutionary trajectory is roughly $P\{ {\mathbf{w}} \} \sim e^{-N\Lambda({\mathbf{w}})}$;
thus, trajectories with smaller cost $\Lambda({\mathbf{w}})$ are more likely.
In addition, the number of steps in the most likely trajectories increases with $w(2)$.

Total trajectory cost as a function of the penultimate point is depicted in Figure \ref{figcost}.
The trajectory cost looks like a single-well potential, and the cost grows rapidly outside of a small region around $y^*$.
Therefore, one can potentially implement a gradient-descent algorithm starting at the point $y^*,v $ which
would significantly accelerate computations compared to a straightforward grid search.
Of course, the single-well shape of the cost function cannot be proven in general and needs to be verified for biological 
models with variables mutation rates, deleterious mutations, and a larger number of genotypes.
If the single-well shape of the cost function does not hold in higher dimensions and/or other parameter values, it is 
possible to implement a probabilistic minimization algorithm similar to the simulated annealing.
\begin{figure}[t]
\centerline{\includegraphics[scale = 0.9]{cost.pdf}}
\caption{Cost of trajectories connecting the initial histogram $H = [0.99, 0.005, 0.005]$ and the final histogram 
$G=[0.35, 0.35, 0.3]$ as a function of the penultimate point. Red point in the Gen1 - Gen2 plane indicates the $y^* \approx[0.475, 0.310, 0.215]$ computed from conjecture \ref{conj:PenPt}.}
\label{figcost}
\end{figure}







\textbf{Illustrating Conjectures \ref{conj:PenPt}, \ref{conj:Discretization}, \ref{conj:Boundary}.}
Next, we illustrate the applicability of Conjectures \ref{conj:PenPt}, \ref{conj:Discretization}, \ref{conj:Boundary}
for a particular target histogram $G=[0.35, 0.35, 0.3]$. We would like to emphasize that numerical results which are presented here
for this particular histogram are generic, and a similar behavior has been confirmed for several other target histograms in \eqref{gtar}. 
For the target histogram $G=[0.35, 0.35, 0.3],$ the solution of the linear system 
\eqref{e:sys} in Conjecture \ref{conj:PenPt} is $y^* \approx [0.47430, 0.31043, 0.21527]$ while the
penultimate point yielding the most likely trajectory for this target histogram is $[0.4752, 0.30821, 0.21659].$
The difference between these two points is $O(10^{-3})$. This is comparable with the coarse discretization $\epsilon=10^{-3}$
which produced quantitatively similar results. For all target points we investigated, solution of the linear system $y^*$ in 
Conjecture \ref{conj:PenPt} was very close to the penultimate point on the most likely trajectory and serves 
as a natural choice for the center of the set of penultimate points.

Figure \ref{fig2} depicts upper 5\% and 10\% quantiles for the values of the function $1/h(y,G)$ discussed in Conjecture 
\ref{conj:Discretization}. Table \ref{tab0} presents values of the constant $c$ for these quantiles for several different target 
histograms. Figure \ref{fig2} demonstrates that the upper 5\% quantile of the function $1/h(y,G)$ is an
adequate choice for the set of penultimate points in the backward search algorithm. The upper 10\% quantile set
appear to be excessive as a choice of penultimate points in the search algorithm since it covers two evolutionary 
points on the  trajectory.  This quantile also covers a considerable region where the frequency of 
genotype 1 is smaller than 0.35. For the target point with $w(1) = 0.35,$ the most likely trajectory would never go through the
region with $w(1) < 0.35.$ Recovering this genotype back to $0.35$ would be a very costly step due to the
fact that this is the genotype with the lowest fitness. 

We observe that values of the constant $c$ in Conjecture \ref{conj:Discretization} for all target points in Table \ref{tab0} 
are extremely small due to the fact that $h(y^*,G) = O(m)$
and quickly increases away from $y^*$. However, Table \ref{tab0} demonstrates that the constant $c$ is of the same order for
different target points and thus roughly independent of the choice of the target histogram $G$. Table \ref{tab0} also
demonstrates that constant $c$ increases slightly when the target histogram moves closer to the boundary.
In general, choosing $c \approx 5\times 10^{-6}$ is an adequate choice to define the initial adaptive set of 
penultimate points in the backward search algorithm for all target histograms $G$.

In practice, implementation of Conjecture 
\ref{conj:Discretization} requires some preliminary exploration of the state space of interior histograms, such as computing
quantiles for a few choices of target histograms, computing the corresponding values of $c$, and plotting regions of
penultimate points for the backward search algorithm. However, such a preliminary exploration of the state space is rather fast 
computationally since it does not require computing any trajectories and thus can be quite beneficial computationally for 
evolutionary problems with a large number of genotypes (e.g., 8 genotypes). We will develop a practical approach using the idea of the adaptive region in subsequent papers for problems with more than 3 genotypes. 

Backward search for the most likely trajectory is rather fast computationally for the evolutionary problem with 3 genotypes. 
Therefore, we also explored the cost of multi-stage trajectories bouncing from the boundary as 
described in Section \ref{ss:multi}. Typically, there is a small number of incomplete trajectories with an 
(incomplete) cost which is smaller than the cost of the one-stage optimal trajectory. In particular, we found only 44 incomplete trajectories with an (incomplete) cost which is lower than the cost of one-stage most likely trajectory for the target histogram 
$G=[0.35, 0.35, 0.3]$. Out of these 44 trajectories, most of them reach the boundary at a point (histogram) with the frequency
of genotype 1 larger than the frequency of genotype 1 in the initial histogram $H$ in \eqref{numpar}.
It is very costly to connect a histogram with $w(1) > 0.99$ to the initial histogram $H$ since this corresponds to the event
(forward in time) when the frequency of genotype 1 increases for a few steps and then starts decreasing again. This is an extremely unlikely event since genotypes 2 and 3 have higher fitnesses, making them more likely to grow faster than genotype 1.

For several incomplete trajectories terminating at $w_{l}$ in reverse time with $w_{l}(1) < 0.99$, the cost of connecting the point $w_{l-1}$ near the 
boundary with the initial histogram $H$ is very costly.
The details of the most likely trajectory connecting $H=[0.99, 0.005, 0.005]$ and $G=[0.35, 0.35, 0.3]$
are presented in Table \ref{tab1}. In particular, the cost of the first step $H \to w_1$ is the largest in the whole trajectory.
This suggests that events in the initial stage (i.e., $H \to w_1 \to w_2 \to \ldots$) of any trajectory essentially determine the 
overall cost of this trajectory. 

Initial costs are highly sensitive to the choice of $w_1$, $w_2$, and so on. 
This leaves a relatively small set of multi-stage trajectories which can have an overall cost
which is lower than the optimal one-stage trajectory. In other words, the numerical investigation of multi-stage trajectories
become computationally cheap after a good optimal one-stage candidate has been found. This implies that investigating 
multi-stage trajectories is a secondary task while the primary problem is finding the optimal one-stage trajectory.
In our simulations, we determined that multi-stage trajectories have costs which are significantly higher than the cost of the one-stage
most likely trajectory. There is a strong indication that events which correspond to ``bouncing away from the boundary'' 
generate costs which are significantly higher than the cost of the initial stages of the one-stage most likely trajectory.

\begin{figure}[th]
\centerline{\includegraphics[scale = 0.5]{Quant5Plot.pdf}
\includegraphics[scale = 0.5]{Quant10Plot.pdf}}
\caption{Upper 5\% (left sub-plot) and 10\% quantiles for the values of $1/h(y,G)$ for the target histogram $G=[0.35, 0.35, 0.3]$.
Note that $y^*$ computed as the solution for the linear system \eqref{e:sys} nearly coincides with the penultimate point for the 
most likely trajectory. The corresponding values of $c$ in Conjecture \ref{conj:Discretization} are 
$c=7.2 \times 10^{-6}$ and $c=4.9 \times 10^{-6}$ for the 
5\% and  10\% upper quantile, respectively.}
\label{fig2}
\end{figure}

\begin{table}
\begin{center}
\begin{tabular}{ l | c | c }
\hline
 & Quant 5\%   &   Quant 10\% \\
\hline
 G=[0.35, 0.35, 0.3]   & $7.2 \times 10^{-6}$  & $4.9 \times 10^{-6}$ \\
 G=[0.35, 0.45, 0.2]   & $9.7 \times 10^{-6}$  & $6.8 \times 10^{-6}$ \\
 G=[0.2, 0.35, 0.45]   & $5.6 \times 10^{-6}$  & $3.7 \times 10^{-6}$ \\
 G=[0.2, 0.5, 0.3]      & $7.7 \times 10^{-6}$  & $5.1 \times 10^{-6}$ \\
 G=[0.2, 0.6, 0.2]	   & $10.5 \times 10^{-6}$ & $7.2 \times 10^{-6}$ \\
\hline
\end{tabular}
\end{center}
\caption{Values of $c$ for the upper 5\% and 10\% quantiles of $1/h(y,G)$ for various target points.}
\label{tab0}
\end{table}



\textbf{Simulations with varying initial histogram $H$.}
We consider the target histogram $G=[0.35, 0.35, 0.3]$ and initial histograms
\begin{equation}
\label{gtar1}
H = (w(1), 1 - w(1)/2, 1 - w(1)/2), \text{~~with~} w(1) = 0.99,\ldots,0.999.
\end{equation}
This corresponds to initial populations with a decreasing frequency of genotypes 2 and 3.
Examples of four most likely trajectories and zoom-ins on the initial part of each trajectory are presented in Appendix \ref{app1}. 
The large-scale aspects of these trajectories look similar. In fact, the difference in penultimate points for these trajectories is
$O(10^{-3})$.  Since all trajectories are completely determined by their penultimate points, the behaviors
of these most likely trajectories near the target histogram $G=[0.35, 0.35, 0.3]$ are very similar. 
However, there are subtle differences between these trajectories 
near the corresponding initial histogram for each trajectory
(Figure \ref{fig3a} in Appendix \ref{app2}).
We also depict trajectory cost and number of steps in each trajectory vs. $w(1)$ in Figure \ref{fig3b}. 
Numerical values of all points (histograms) in the most likely trajectories for two starting histograms with 
$w_0(1)=0.99$ and $w_0(1)=0.999$ together with the cost of each step are presented in Tables \ref{tab1} and \ref{tab2} 
in Appendix \ref{app2}.
\begin{figure}[th]
\centerline{\includegraphics[scale = 0.5]{fig4a.pdf}
\includegraphics[scale = 0.5]{fig4b.pdf}}
\caption{Left sub-plot - Cost of the optimal trajectory vs. $w(1)$ in the initial histogram $H$ in \eqref{gtar1}.
Right sub-plot - number of steps in the optimal trajectory vs. $w(1)$ in the initial histogram $H$ in \eqref{gtar1}.}
\label{fig3b}
\end{figure}




The most likely trajectory reaching $G=[0.35, 0.35, 0.3]$ starting with $H = [0.99, 0.005, 0.005]$ is considerably more costly
compared with the most likely trajectory starting with $H = [0.999, 0.0005,$ $0.0005]$. Therefore, initial populations with 
a higher frequency of the genotype with the lowest fitness are more likely to exhibit rare events. There are
also similarities in both trajectories. In particular, the cost of the few initial steps (e.g., $H \to w_1 \to w_2$) 
contributes heavily to the overall cost of the trajectory. The fraction of the total cost for each step in the two trajectories 
for different initial histograms is depicted in Figure \ref{fig5}.
In particular, the first step $H \to w_1$ corresponds to the
event when the frequency of genotype 2 increases much more than the frequency of genotype 3. 
Obviously, this is a rare event with a low probability. 
This trend continues in subsequent steps initially so that the frequency genotype 2 reaches large values compared to genotype 3.
Such sequences of rare events determine the most likely trajectory describing the global rare event of reaching
the target histogram $G$.
\begin{figure}[th]
\centerline{\includegraphics[scale = 0.5]{fig5a.pdf}
\includegraphics[scale = 0.5]{fig5b.pdf}}
\caption{Fraction of the total cost per step for optimal trajectories starting at the initial histogram 
$H_0$ and reaching the final histogram  $G=[0.35, 0.35, 0.3]$. 
Left sub-plot - trajectory starting at the initial histogram $H=[0.99, 0.005, 0.005]$ 
(total cost of this trajectory = 3.630589e-03).
Right sub-plot - trajectory starting at the initial histogram $H=[0.999, 0.0005, 0.0005]$
(total cost of this trajectory = 5.721059e-04). }
\label{fig5}
\end{figure}
The costs of each step decay exponentially during the initial stages of the trajectory.
We would like to point out that first step in both trajectories is approximately 47\% of the total cost, and the first two steps
constitute approximately 70\% of the total cost. This implies that in most likely trajectories, radical shifts in populations occur early 
and determine the subsequent evolution of the population.

\section{Conclusions}

In this paper we presented a novel approach for computing the most likely evolutionary trajectories for rare events
in large bacterial populations described by the ``locked-box'' type models. The approach developed here
relies on Large Deviations Theory. In particular, the two most important theoretical results are
Theorems \ref{thCHG} and \ref{MinPathInt}. 
An explicit formula for the cost function for the transition between two interior histograms is derived in Theorem \ref{thCHG}.
This explicit formula allows for a subsequent construction of the backward search algorithm in Theorem \ref{MinPathInt}.
These two theorems lay a proper theoretical foundation for a practical implementation of the search algorithm for the most likely
evolutionary trajectory connecting two interior histograms.

Numerical implementation of the search algorithm is quite subtle as described in Sections \ref{ss:stage1} and \ref{ss:multi}. In particular, 
boundary effects need to be potentially taken into account, multi-stage trajectories significantly increase computational time, along with other issues.
However, we demonstrated in our numerical examples that all these potential obstacles can be overcome by developing a ``smart''
search algorithm. Conjectures \ref{conj:PenPt} - \ref{conj:Boundary} significantly accelerate the search algorithm, and the 
Large Deviations approach developed here yields an efficient practical tool for investigating various evolutionary scenarios.

We used our approach to investigate a biological system with three genotypes. In particular, our extensive numerical investigation 
of this system demonstrated the validity of Conjectures \ref{conj:PenPt} - \ref{conj:Boundary}. Moreover, this example 
elucidates very interesting behavior of the biological system under investigation. For instance, most likely trajectories computed for this example demonstrate that few initial steps correspond to a very large percentage of the total cost. This implies that radical shifts in the 
frequency of genotypes occur early for rare events and determine the subsequent evolution of the population and
emergence of genotypes with a lower fitness.

The numerical approach developed here has been verified on a prototype example with non-reversible advantageous mutations.
We also investigated the effect of allowing reversible mutations. We obtained quantitatively similar results 
(not presented here for the brevity of presentation). In particular, most likely trajectories look quite similar to the examples with 
non-reversible mutations presented here. We will address the importance of 
reversible mutations in biological systems with more genotypes in subsequent publications.

The backward search algorithm can be significantly accelerated using Conjectures \ref{conj:PenPt} - \ref{conj:Boundary}.
With such accelerated methodology, it stands ready to be applied to more realistic biological systems with up to 8 genotypes.
We plan to address the following important questions in subsequent papers - 
(i) importance of variable mutation rates for different genotypes, 
(ii) deleterious mutations leading to the decrease in fitness, 
(iii) mutation matrices representing realistic biological scenarios. 
In particular, biological examples with 8 genotypes represent bacterial systems with 3 possible mutations. 
Our numerical approach stands ready to tackle such computationally intensive models.  
Quantitative investigation of most likely evolutionary trajectories in such realistic examples can potentially shed light on
mechanisms for rare events (e.g. emergence or fixation of non-dominant genotypes) in these systems.






\clearpage


\subsection{Interior Histograms} \label{interior}
An histogram $H \in \mathcal{H}$ will be called \emph{interior} if $H(j) > 0$ for all $1 \leq j \leq g$. We denote $\mathcal{H}^\circ$ as the set of interior histograms.  Note that for $H,G\in\mathcal{H}^\circ,$  the transition cost $C(H,G)$ is always finite. A path $\mbold{h} = [h_1 \ldots h_T ] $ will be called an \emph{interior path} if $h_n\in \mcal{H}^\circ$ are interior histograms for all $1 \leq n \leq T$.  To develop explicit computational schemes, the remaining of this paper will from now on, focus only on interior paths and histograms.
\subsection{Sets of Thin Tubes realizing rare events}
For $\eta>0$ and $E^* \subset  \Omega_T$, define the $\eta$-neigborhood $U^{\eta}(E^*)$ of $E^*$ as the union of all open balls of radius $\eta$ and center in $E ^*$. For small $\eta$, the set  $U_{\eta}(E^*)$ is a set of thin tubes of paths with ``axes'' $h \in E^*$.  Let $\mbold{H} \in \Omega_T $ be the random trajectory of population histograms.  For any closed set of interior paths $E \subset \Omega_T$, one can easily prove that $a = b(E) >0$. When  $\Lambda(E) > 0$, Theorem \ref{thsetLD} provides the fast-vanishing bounds
$$
P(\mbold{H} \in E) \leq P(\mbold{H} \in V_N(E)) \leq 2 e^{- N A}
$$ 
for any $0 <A < \Lambda(E)$ and $N > N_0(a, A, \mcal{P})$ so that $\{\mbold{H} \in E\}$ and  $\{\mbold{H} \in V_N(E)\}$ are \emph{rare events}.
\begin{theorem} \label{conditioning} Fix the path length T and an initial histogram H. Denote $P_H$ as the probability distribution of random histogram paths $\mbold{H} \in \Omega_T$ starting at H. Let $E \subset \Omega_T $ be any closed set of interior paths starting at H satisfying $0 < \Lambda(E) < \infty $. Let $E^*$ be the set of all paths $\mbold{h}$ minimizing the rate function $\lambda(\mbold{h})$ over all $\mbold{h} \in E$. Then $E^*$ is a closed subset of E. For any fixed $ \eta > 0$, the $\eta$-neighborhood $U= U^{\eta}(E^*)$ verifies 
$$
\lim_{N \to \infty}  P_H \left(\mbold{H} \in U \; | \; \mbold{H} \in V_N(E) \right)   =  1
$$
with exponential speed of convergence.
\end{theorem}
\begin{proof}
Set $a = b(E) > 0 $ and $L = \Lambda(E) > 0$. Let  $E(c) = \vset{ \mbold{h} \in E \; | \; \lambda(\mbold{h}) \leq c }$. Then the function $\lambda(\mbold{h})$ is continuous on $E(2 L)$ due to Theorem \ref{holderlambda} applied to interior paths. Hence, $E(2 L)$ is closed and must contain any path minimizing  $\lambda(\mbold{h})$ over $E(2 L).$ However, the two sets of minimizers of $\lambda$ over $E$ and over $E(2 L)$ are obviously identical. This proves $E^* \subset E(2 L) \subset E$.

Theorem \ref{holderlambda} applied to interior paths gives a constant $c = c (T,L,a,\mcal{P})$ such that for all $\mbold{h} \in E^*$ and $\mbold{h}' \in E$ with $\vnorm{ \mbold{h}' - \mbold{h} } < a/F_g$, one has $\abs{ \lambda(\mbold{h}') - \lambda(\mbold{h})} \leq c \vnorm{\mbold{h}'-\mbold{h}}^{1/2}.$ Fix any $0 < \eta < a/F_g$. Let $K \subset E$ be the open $\eta$-neighborhood  of $E^*$ within $E,$  and set $W= E - K.$ For each $\mbold{h}' \in K$, there is one $\mbold{h} \in E^*$ with $\vnorm{\mbold{h}' - \mbold{h}} < \eta.$ Therefore, Theorem \ref{holderlambda} implies $\lambda(\mbold{h}') < L + c \sqrt{\eta}.$  This forces $\lambda(\hat{\mbold{h}}) \geq L + c \sqrt{\eta}$ for all $\hat{\mbold{h}} \in W$ so that $\Lambda(W) \geq L + c \sqrt{\eta}$.  Apply Theorem \ref{thsetLD} to get $N_0$ and $c_0$ determined by $(T,L,a,PAR,\eta)$ such that for $N >N_0$, 
\begin{align*}
\frac{1}{N}  \log P_H (\mbold{H} \in V_N(W) ) & =  - L - c \sqrt{\eta}+o_1(N), \\
\frac{1}{N}  \log P_H (\mbold{H} \in V_N(E) )  & =  - L +o_2(N), 
\end{align*} 
with  $\abs{o_1(N)} < c_0 / \sqrt{N}$ and $\abs{o_2(N)} < c_0 / \sqrt{N}.$ These  results yield
\begin{align}
P_H (\mbold{H} \in V_N(W)) & \leq \exp(-N L - N c \sqrt{\eta} +c_0 \sqrt{N}),\\
P_H (\mbold{H} \in V_N(E)) & \geq \exp( -N L  - c_0 \sqrt{N}),
\end{align}
so that
$$
\frac{P_H (\mbold{H} \in V_N(W))}{P_H (\mbold{H} \in V_N(E))}  \leq \exp(- N c \sqrt{\eta} + 2 c_0 \sqrt{N}).
$$
Impose $N > N_0 + {16 c_0^2}{c^2 \eta}$ to force  $N (c/2) \sqrt{\eta} > 2 c_0 \sqrt{N}$ to give 
$$
\frac{P_H (\mbold{H} \in V_N(W))}{P_H (\mbold{H} \in V_N(E))} \leq \exp(- N (c/2) \sqrt{\eta}).
$$
Since $E = W \cup K,$ one has $V_N(E) \subset  V_N(W) \cup V_N(K).$ Therefore, 
$$
P_H (\mbold{H} \in V_N(K)) \geq P_H (\mbold{H} \in V_N(E)) - P_H (\mbold{H} \in V_N(W)).
$$
From $ P_H (\mbold{H} \in V_N(K) \; | \; \mbold{H} \in V_N(E)) = \frac{P_H (\mbold{H} \in V_N(K))}{P_H (\mbold{H} \in V_N(E))}$, we now get
$$
1 \geq P_H (\mbold{H} \in V_N(k) \; | \; \mbold{H} \in V_N(E) ) \geq 1 - \frac{P_H (\mbold{H} \in V_N(K))}{P_H (\mbold{H} \in V_N(E))} \geq 1 - \exp\left(- N (c/2) \sqrt{\eta}\right)
$$ 
for $N > N_1.$ For $ N > N_1 +1/\eta$, one has $V_N(K) \subset U^{2 \eta}(E^*) \cap  V_N(E).$ Therefore,
$$
1 \geq P_H (\mbold{H} \in U^{2 \eta}\; | \; \mbold{H} \in V_N(E) ) \geq P_H (\mbold{H} \in V_N(K) \; | \;  \mbold{H} \in V_N(E)) \geq 1 - \exp\left(- N (c/2) \sqrt{\eta}\right),
$$
which concludes the proof.
\end{proof}

So for large population size $N$, rare evolutionary events $\{\mbold{H} \in E\}$ with finite $\Lambda(E) > 0$ can only be realized by population evolutions following very thin tubes around the paths $\mbold{h}^*$ minimizing $\lambda(\mbold{h})$ over paths $\mbold{h} \in E$. Computing such paths requires efficient numerical strategies discussed in Section \ref{s:numtraj}.
\subsection{Most likely path connecting two histograms} \label{mostlikelypath}
\begin{definition} \label{reachable}
Let $Q$ be the transfer matrix quantifying mutation rates $M = mQ$. For any non-empty set $S$ of genotypes, recall the set $R(S)$ of genotypes reachable from $S$ defined in Definition \ref{reachgen}.  For any histograms  $H$ and $G$, we say that \emph{G is reachable from H} if $spt(G) \subset R(spt(H))$.
\end{definition}
Note that if there exists some  power $p$ such that all coefficients of $Q^p$ are positive,  one has $R(S) = \vset{ 1, \dots, g }$ for any non- empty set of genotypes $S$, and hence, any  $G$ is reachable from any $H$.  
Our next theorem answers an important question for bacterial genetic evolution:  how can one reconstitute the \emph{most likely} evolutionary path starting at a known initial histogram $H$ and reaching a known terminal histogram $G$ after an unknown number $T$ of daily cycles. 
\begin{theorem}\label{tube}
Fix any histograms $H$ and $G$ such that $G$ is reachable from $H$. Fix $0 < a <1$. Define the set of paths $E_T$ 
$$
E_T = \vset{\mbold{h} \in \Omega(T) \,|\, b(\mbold{h}) \geq a ;  \, h_1 = H  ; \,  h_T = G  ; \, h_n \neq G \text{ for } 1 \leq n < T}.
$$ 
Let $E(H,G) = \bigcup_t E_T$ be the set of all paths starting at $H$ and hitting $G$ at some finite time $T$.  Then $L(H,G) = \Lambda (E(H,G) ) = \inf_T \, \inf_{\mbold{h} \in E_T} \, \lambda(h)$ must be finite. If $H, G\in \mcal{H}^\circ$ and the set $E^* = \vset{\mbold{h}^* \in E(H,G) \,|\, \lambda(h^*) = L(H,G)}$ contains only interior paths of lengths inferior to some finite $T^*$, then any open neighborhood U of $E^*$ verifies
$$
\lim_{N \to \infty}  \frac{1}{N} \log P (\mbold{H} \in U\; | \; \mbold{H} \in E(H, G)) = 1
$$
with convergence at exponential speed.
\end{theorem}
The proof is similar to the proof of the previous theorem and will be omitted.





\section{Computation of cost-minimizing histograms trajectories}\label{s:numtraj}
\subsection{Geodesics in the space of histograms}
To identify the most likely thin tubes of paths linking histograms $H$ and $G$ in a given finite time $T$, thus outlining potential bacterial evolution scenarios from $H$ to $G$, one needs to compute discretized paths $\mbold{h}^* = [h_1^* \, h_2^* \dots h_T^*]$ minimizing the large deviations cost $\lambda(\mbold{h})$ over all $\mbold{h} \in \Omega_T$ such that $h_1 = H$ and $h_T = G$. We call any such $\mbold{h}^*$ a \emph{geodesic} from $H$ to $G$ provided $\lambda(\mbold{h}^*)$ is finite.  When $h_n \in \mcal{H}^\circ$ for all $1 \leq n \leq T$, we call $\mbold{h}^*$ an \emph{interior geodesic}.  Computing geodesics presents numerical and mathematical challenges. We now develop an efficient theoretical approach to iteratively generate geodesics.
\subsection{Explicit computation of geodesics}
\begin{theorem}\label{MinPathInt}
Let $\mbold{h}^*$ be any interior geodesic in $\Omega_T$ with $T \geq 2$. Denote $a = b(\mbold{h}^*)$. There is a constant $m_0 = m_0(a,\mcal{P})$ such that for $m<m_0$ and any $1 \leq n \leq T-2$, the histogram $x:= h^*_n$ is fully determined by $y := h^*_{n+1}$ and $z := h^*_{n+2},$ which is given by $x = \chi(m,y,z)$ where $\chi$ is a $C^{\infty}$ function of (m,y,z) for $m<m_0$ and $y,z \in \mathcal{H}^\circ.$  Hence, for $m < m_0$, the geodesic $\mbold{h}^*$ is determined by its last two points $h^*_T$ and $h^*_{T-1}$ thanks to the reverse recurrence relation
\begin{equation}\label{reverse}
h^*_n = \chi(m, h^*_{n+1}, h^*_{n+2}), \qquad  1\leq n \leq T-2.
\end{equation}
Denote $x_s \simeq \wh{x}_s + m \wh{x}_s w_s$ for $s \in \vset{1, \dots, g}$ with remainder of order $m^2$ the first-order Taylor expansion of $x = \chi(m,y,z)$ in $m$ for $m < m_0$. The interior histogram $\wh{x}$ and the vector $w$ depend only on $y,z$ and are given below by the explicit formulas \eqref{X}, \eqref{ZeroOrder}, \eqref{w}, \eqref{alphas}, and  \eqref{betas}.
\end{theorem}
\begin{proof}
Any sub-segment $[h^*_n\, h^*_{n+1}\, h^*_{n+2}]$ of $\mbold{h}^*$ is also a geodesic from  $x =h^*_n $ to $z = h^*_{n+2}$. Hence, given $(x, z)$, the two-step cost function $u(x,y,z) = C(x,y) + C(y,z)$ is minimized in $y$ by $y = h^*_{n+1}$.  For any three histograms $x,y,z \in \mcal{H}^\circ,$ both $C(x,y)$ and $C(y,z)$ are finite and differentiable in $y$ by Theorem \ref{thCHG} and convex in $y$ by Lemma \ref{finitecost}. Hence, for fixed $x,z \in \mathcal{H}^\circ$, the function $f(y) = C(x,y) + C(y,z)$ is finite, convex, and differentiable for all $y$ in the open convex set $\mathcal{H}^\circ$.  If $y \in \mathcal{H}^\circ$ is a minimizer of $f(y)$ over all $\mathcal{H}^o$, any vanishingly small  modification $\Delta y$ of $y$ within $\mathcal{H}^\circ$ must verify $\sum_k \Delta y_k =0$ and $f(y+\Delta y) \geq f(y)$. Hence the gradient of $\left[f(y) + \mu  \sum_k  y_k \right]$ must be 0 for some Lagrange multiplier $\mu$. For each $s\in \vset{1, \dots, g},$ this yields the system
\begin{equation}\label{twosteps}
\frac{\partial}{\partial y_s} f(y) + \mu= 0.
\end{equation}
Denote $u(x,y,z) = f(y) = C(x,y) + C(y,z)$. For given  $y, z \in \mathcal{H}^\circ$, extend \eqref{twosteps} into the following system of $(1+g)$ equations to be solved for a histogram $x$ and a Lagrange multiplier $\mu$:
\begin{align}
\frac{\partial}{\partial y_s} u(x,y,z) + \mu &= 0, \qquad  1 \leq s \leq g,   \label{s1} \\
\sum_k  x_k  &=1.  \label{s2}
\end{align}
This provides $(g+1)$ equations for $(g+1)$ unknowns $(x,\mu)$. We now show that for $m=0,$ this system has a unique explicit solution $(\wh{x}, \nu)$ before applying the implicit function theorem. By Theorem \ref{thCHG}, for $m < m_0$ with $m_0 > 0$ small enough, the function $u(x,y,z)$ is $C^{\infty}$ in $(m, x,y,z)$, with  explicit first-order Taylor expansion in $m$ given by \eqref{formulaCHG}. Recall our earlier notations
\begin{align*}
E_{s,k}(x,y) &= \exp \left( -\frac{y_s}{F_s x_s} + \frac{y_k}{F_k x_k} \right),  \\
E_{s,k}(y,z) &= \exp \left( -\frac{z_s}{F_s y_s} + \frac{z_k}{F_k y_k} \right). 
\end{align*}
Taking derivatives in $y$ of the first-order expansions \eqref{formulaCHG} of $C(x,y)$ and $C(y,z)$ readily yield the following first-order Taylor expansions in $m$ with remainders of order $m^2 \leq 10^{-12}:$  
\begin{align} 
\frac{\partial}{\partial y_s} C(x,y) &\simeq A_s (x,y) + m \wh{A}_s(x,y), \label{DyCxy1}\\ 
A_s(x,y) &= 1 + \log\frac{y_s}{F_s} - \log x_s, \label{DyCxy2}\\
\wh{A}_s(x,y) &= \sum_k \left( Q_{s,k} E_{s,k}(x,y) - \frac{F_k x_k}{F_s x_s} Q_{k,s}E_{k,s} \right), \label{DyCxy3}\\
\frac{\partial}{\partial y_s} C(y,z) &\simeq D_s(y,z) + m \wh{D}_s(y,z),  \label{DyCyx1} \\
D_s(y,z)  &= \frac{F_s}{\langle F,y \rangle} - \frac{z_s}{y_s}, \label{DyCyx2} \\
\wh{D}_s(y,z) &= F_s \sum_k Q_{s,k} - \left( F_s + \frac{z_s}{y_s} \right) \sum_k  E_{s,k}(y,z) Q_{s,k} - \frac{z_s}{F_s y_s^2} \sum_k  F_k y_k Q_{k,s} E_{k,s}(y,z). \label{DyCyx3}
\end{align}
For  $m= 0$,  the system in \eqref{s1} becomes $A_s(x,y) +D_s(y,z) + \mu =0$ for each $s$, which 
yields
\begin{equation*}
1 + \log\frac{y_s}{F_s x_s} + \frac{F_s}{\langle F, y \rangle} - \frac{z_s}{y_s} + \mu = 0.
\end{equation*}
Hence, we have $x = e ^{1+\mu} X$ where the vector  $X$ is given by 
\begin{equation} \label{X}
X_s = \frac{y_s}{F_s}  \exp \left( \frac{F_s}{\langle F, y \rangle} - \frac{z_s}{y_s} \right) >0.
\end{equation}
The constraint given by \eqref{s2} gives $e^{1+\mu} = \frac{1}{\sum_t  X_t}$. Therefore, for $m=0$ and all $s,$ the unique solution $\wh{x} = x(0,y,z)$ and $\wh{\mu}= \mu(0,y,z)$ of the system given by \eqref{s1} and  \eqref{s2} is
\begin{equation}\label{ZeroOrder}
\wh{x}_s = \frac{X_s}{\sum_t  X_t}, \qquad \qquad  X_s = \frac{y_s}{F_s}  \exp \left( \frac{F_s}{\langle F, y \rangle} - \frac{z_s}{y_s} \right), \qquad \qquad 1+ \wh{\mu} = - \log\left(\sum_t X_t\right).
\end{equation}
Denote $U(s,t) = \frac{\partial}{\partial x_t}  \frac{\partial}{\partial y_s}  u(x,y,z) + \mu$. For $m = 0,$ one has 
$$
U(s,t) = \frac{\partial}{\partial x_t}  \frac{\partial}{\partial y_s} [ A_s(x,y) +D_s(y,z) ] =  - 1_{\vset{s = t}}  \frac{1}{\wh{x}_s}.
$$ 
The $g \times g$ matrix $U$ is thus diagonal with non-zero diagonal terms, making it invertible. Therefore, the implicit function theorem applies to the system given by \eqref{s1} and \eqref{s2}. Hence, for some fixed $m_1 >0$, there is a unique solution $(x(m, y,z), \mu(m,y,z))$ to the system \eqref{s1}--\eqref{s2}, and the functions  $(x,\mu)$ are $C^{\infty}$ in $(m, y,z)$.

Let $x \simeq \wh{x} + m V$ and $\mu \simeq \wh{\mu} + m \mu_1$ be the first-order Taylor expansions of $(x, \mu)$ in $m$.  Denote $V_s = \wh{x}_s w_s$ so that $x_s \simeq \wh{x}_s (1 + m w_s)$. The constraint \eqref{s2} then implies $\langle \wh{x}, w \rangle = 0$.  The first-order expansion of \eqref{s1} becomes
\begin{equation}\label{s3}
A_s(\wh{x} + mV , y) + D_s(y,z) + \wh{\mu} + m (\wh{A}_s(\wh{x} , y) +  \wh{D}_s(y,z) + \mu_1) \simeq 0.
\end{equation}
Since $A_s(\wh{x} + m V , y) \simeq A_s(\wh{x},y) - m w_s$, the zero-order term in \eqref{s3} vanishes due to the values of $\wh{x}$ and $\wh{\mu}$. The first-order term must vanish as well, which gives for all $s$
$$
w_s = \wh{A}_s(\wh{x} , y) +  \wh{D}_s(y,z) + \mu_1.
$$
Since $ \langle \wh{x}, w \rangle = 0$, this yields 
$$
\mu_1 = - \sum_t \; \wh{x}_t [ \wh{A}_t(\wh{x} , y) +  \wh{D}_t(y,z)].
$$
We then have
$$
w_s = \wh{A}_s(\wh{x} , y) +  \wh{D}_s(y,z) - \sum_t   \wh{x}_t [\wh{A}_t(\wh{x} , y) +  \wh{D}_t(y,z) ]. 
$$
Define vectors $\alpha$ and $\beta$  by $\alpha_s = \wh{A}_s(\wh{x},y)$ and $\beta_s = \wh{D}_s(y,z)$ so that
\begin{equation} \label{w}
w = \alpha + \beta - \langle \wh{x} , \alpha +\beta \rangle.
\end{equation}
The expressions for $\wh{A}_s, \wh{D}_s,$ and $E_{s,k}$ given above yield directly 
\begin{equation} \label{alphas}
\alpha_s =  \sum_k  \left( Q_{s,k} e_{s,k} - \frac{F_k X_k}{F_s X_s} Q_{k,s} e_{k,s} \right)
\end{equation}
with $ e_{s,k}=1/e_{k,s} = \exp \left[ -y_s/(F_s \wh{x}_s) + y_k/(F_k \wh{x}_k)\right]$ and
\begin{equation} \label{betas}
\beta_s = F_s \sum_k Q_{s,k} - \left( F_s + \frac{z_s}{y_s} \right) \sum_k  f_{s,k} Q_{s,k} - \frac{z_s}{F_s y_s^2} \sum_k  F_k y_k Q_{k,s} f_{k,s}
\end{equation}
with $f_{s,k}=1/f_{k,s} = \exp \left[-z_s/(F_s y_s) + z_k/(F_k y_k)\right]$.
Note that $\wh{x}, \alpha, \beta,$ and $w$ depend only on $y$ and $z$. The preceding formulas provides the explicit first-order expansion $x_s = \wh{x}_s + m \wh{x}_s w_s$, concluding the proof.
\end{proof}
\section{Most Likely Trajectories for the Target histogram $G=[0.35, 0.35, 0.3]$ and varying initial histogram $H$}
\label{app1}

\begin{figure}[th]
\centerline{\includegraphics[scale = 0.5]{fig3a.pdf}
\includegraphics[scale = 0.5]{fig3b.pdf}}
\centerline{\includegraphics[scale = 0.5]{fig3c.pdf}
\includegraphics[scale = 0.5]{fig3d.pdf}}
\caption{Most likely evolutionary trajectories for the population with the target histogram $G=[0.35, 0.35, 0.3]$
and four different starting histograms.
Upper left subplot - starting histogram $H=[0.99, 0.005, 0.005]$.
Upper right subplot - starting histogram $H=[0.992, 0.004, 0.004]$.
Lower left subplot - starting histogram $H=[0.996, 0.002, 0.002]$.
Lower right subplot - starting histogram $H=[0.999, 0.0005, 0.0005]$.}
\label{fig3}
\end{figure}
\begin{figure}[th]
\centerline{\includegraphics[scale = 0.5]{fig31a.pdf}
\includegraphics[scale = 0.5]{fig31b.pdf}}
\centerline{\includegraphics[scale = 0.5]{fig31c.pdf}
\includegraphics[scale = 0.5]{fig31d.pdf}}
\caption{Zoom-in on the initial part of most likely evolutionary trajectories for the population 
with the target histogram $G=[0.35, 0.35, 0.3]$
and four different starting histograms.
Upper left subplot - starting histogram $H=[0.99, 0.005, 0.005]$.
Upper right subplot - starting histogram $H=[0.992, 0.004, 0.004]$.
Lower left subplot - starting histogram $H=[0.996, 0.002, 0.002]$.
Lower right subplot - starting histogram $H=[0.999, 0.0005, 0.0005]$.}
\label{fig3a}
\end{figure}





\section{Most Likely Trajectories for Initial Histograms $[0.99, 0.005, 0.005]$ and $[0.999, 0.0005, 0.0005]$}
\label{app2}

\begin{table}
\begin{center}
\begin{tabular}{ l | c | c | c | c  }
\hline
Point \# &  Gen 1 & Gen 2 & Gen 3 & Cost(${\rm Point}_{i}$, ${\rm Point}_{i+1}$) \\
\hline
1 & 9.900000e-01  & 5.000000e-03  & 5.000000e-03  & 1.686764e-03 \\
2 & 9.858695e-01  & 9.477103e-03  & 4.653446e-03  & 8.536793e-04 \\
3 & 9.780684e-01  & 1.647496e-02  & 5.456620e-03  & 4.539794e-04 \\
4 & 9.652908e-01  & 2.709791e-02  & 7.611302e-03  & 2.515257e-04 \\
5 & 9.452312e-01  & 4.289898e-02  & 1.186983e-02  & 1.448337e-04 \\
6 & 9.144852e-01  & 6.584684e-02  & 1.966793e-02  & 8.682526e-05 \\
7 & 8.685775e-01  & 9.805712e-02  & 3.336541e-02  & 5.454721e-05 \\
8 & 8.026050e-01  & 1.410112e-01  & 5.638382e-02  & 3.627721e-05 \\
9 & 7.130923e-01  & 1.940756e-01  & 9.283203e-02  & 2.581977e-05 \\
10 & 6.011465e-01  & 2.527446e-01  & 1.461089e-01 & 1.983269e-05 \\
11 & 4.752100e-01  & 3.082000e-01  & 2.165900e-01 & 1.650544e-05 \\
12 & 3.500000e-01  & 3.500000e-01  & 3.000000e-01 &    \\
\hline
\end{tabular}
\end{center}
\caption{Most likely trajectory connecting the starting histogram $H = [0.99, 0.005, 0.005]$ and the target histogram 
$G = [0.35, 0.35, 0.3]$. Length of the trajectory = 12, total cost = 3.630589e-03.}
\label{tab1}
\end{table}
\begin{table}
\begin{center}
\begin{tabular}{ l | c | c | c | c  }
\hline
Point \# &  Gen 1 & Gen 2 & Gen 3 & Cost(${\rm Point}_{i}$, ${\rm Point}_{i+1}$) \\
\hline
1 & 9.990000e-01  & 5.000000e-04  & 5.000000e-04  & 2.690597e-04 \\
2 & 9.987909e-01  & 8.722365e-04  & 3.368221e-04  & 1.396945e-04 \\
3 & 9.982657e-01  & 1.461233e-03  & 2.730911e-04  & 6.773913e-05 \\
4 & 9.973505e-01  & 2.365339e-03  & 2.841971e-04  & 3.462357e-05 \\
5 & 9.958865e-01  & 3.746982e-03  & 3.664701e-04  & 1.879674e-05 \\
6 & 9.935968e-01  & 5.852370e-03  & 5.508047e-04  & 1.090245e-05 \\
7 & 9.900367e-01  & 9.051777e-03  & 9.115416e-04  & 6.818925e-06 \\
8 & 9.845097e-01  & 1.389603e-02  & 1.594310e-03  & 4.648979e-06 \\
9 & 9.759384e-01  & 2.119084e-02  & 2.870799e-03  & 3.472284e-06 \\
10 & 9.626815e-01  & 3.208268e-02  & 5.235793e-03  & 2.821475e-06 \\
11 & 9.423048e-01  & 4.812795e-02  & 9.567258e-03  & 2.450309e-06 \\
12 & 9.113674e-01  & 7.126889e-02  & 1.736372e-02  & 2.224294e-06 \\
13 & 8.654191e-01  & 1.035532e-01  & 3.102764e-02  & 2.066413e-06 \\
14 & 7.996234e-01  & 1.463457e-01  & 5.403090e-02  & 1.930059e-06 \\
15 & 7.105730e-01  & 1.988617e-01  & 9.056533e-02  & 1.787137e-06 \\
16 & 5.993759e-01  & 2.564621e-01  & 1.441620e-01  & 1.624884e-06 \\
17 & 4.743500e-01  & 3.102900e-01  & 2.153600e-01  & 1.445064e-06 \\
18 & 3.500000e-01  & 3.500000e-01  & 3.000000e-01 & \\
\hline
\end{tabular}
\end{center}
\caption{Most likely trajectory connecting the starting histogram $H = [0.999, 0.0005, 0.0005]$ and the target histogram 
$G = [0.35, 0.35, 0.3]$. Length of the trajectory = 18, total cost = 5.721059e-04.}
\label{tab2}
\end{table}

\subsection{Main focus of our numerical examples}
In our numerical examples, the goal is to compute geodesics starting at a given interior histogram $H$ and ending at a target histogram $G^*$ where $G^*$ implements the fixation of some specific genotype $j^*$.  More practically, we seek geodesics starting at a given $H \in \mathcal{H}^\circ$ and targeting interior histograms $G$ such that $G(j^*) > thr$ for some threshold $thr$. In our examples we take $thr = 0.95,$ but pragmatic values of interest may involve lower thresholds $thr \simeq 0.7$. We consider only $j^* < g$ so that $j^*$ has smaller  growth factor $F_j$ than the fittest genotype $g$.  Then for every $T$,  the random event $\vset{ H_1 = H } \cap \vset{ H_T = G }$ is a rare event with exponentially-fast vanishing probability as $N \to \infty$. The random paths realizing this rare event must then be very close to a geodesic $\mbold{h}^*$ linking $H$ to $G$ in $U$ with extremely high probability. We now outline our  numerical search for such geodesics. 
\subsection{First-stage Geodesic Search}\label{ss:stage1} 

{\bf Starting Zone:} In $\mathcal{H}^0$, fix an initial $H$ and a target $G$. By Theorem \ref{thzerotraj}, the zero cost trajectory $\mbold{v} = [v_1\, v_2\, \ldots] $ starting at $v_1 =H$ is recursively generated by $v_{n+1} = \zeta(v_n)$ for $n \geq1$. The function  $\zeta$ is explicitly given by \eqref{zeta}, which shows that $v_n \in \mathcal{H}^0$ for all $n$.  Any subsegment $[v_1\, v_2 \ldots v_n]$ of $\mbold{v}$ has zero cost. Generally, $\mbold{v}$ has an infinite number of steps, but for  $m < m_0$, Theorem \ref{hinfty} shows that  $\lim_{n \to \infty } v_n $ exists. Hence, for any $\rho,$ denote the \emph{starting zone} $STZ(\rho) \subset \mathcal{H}^\circ $ the open set of all $w \in \mathcal{H}^0$ such that $\vnorm{w - v_n} < \rho$ for some $n$.
 
Fix any very small  $\eta > 0$. The continuity of $\lambda$ on compact sets of interior paths provides $\rho = \rho(\eta) >0$ such that for any $w \in STZ(\rho),$ there is an $v_n$ with one-step cost $C(w,v_n) < \eta$.  Fix $\rho = \rho(\eta)$ and the starting zone $STZ = STZ(\rho).$ Fix a finite $n_{ter}$ such that $\vnorm{ v_n - v_{n_{ter}}} < \rho/2 $ for all $n > n_{ter}$. For any $w \in STZ,$ there is then an $v_n$ with $ n \leq n_{ter}$ such that the path $\boldsymbol{\phi} = [v_1 \ldots v_n  \, w]$ connects $H$ to $w$ at nearly zero cost $\lambda(\boldsymbol{\phi}) < \eta$.\smallskip \\
{\bf Truncated Reverse Geodesics:} Fix a finite subnet $S \subset \mathcal{H}^0$ with small mesh size  $\epsilon > 0$ and cardinal of order $1/\epsilon^g$ such that the balls of radius $\epsilon$ and centers in $S$ cover all of $\mathcal{H}$. Potential penultimate points will be chosen from a subset $PEN \subset S$ to be specified further on.  To each $y \in S$, associate the reverse geodesic $\mbold{z} = \mbold{z}(y)$ with $z_1=G$, $z_2 = y$, and interior histograms $z_k$ iteratively defined for all $k \geq 0$ by the recursive equation \eqref{reviterate}.

Since accuracy bounds of the form $c \log(N)/N$  emerged with various constants $c$ in all basic \emph{uniform} large deviations inequalities proved at the beginning of this paper, we call \emph{near-boundary histograms} all interior histograms $x$ with essential minimum $b(x) < \log(N)/N$.  For fixed small mutation rate $m < m_0$, the recursive equation  \eqref{reviterate} does not necessarily remain valid if $b(z_k) $ tends to 0. Therefore, we stop computing $z_k(y)$ at the first integer $k = k_{ter}$ such that  $z_{k+1}$ is a near-boundary histogram. If no such finite $k$ exists, we set $k_{ter} =  \infty$.

For each path $\mbold{z}(y)$, compute a finite \emph{truncation time} $\kappa(y) \leq k_{ter}$ and a \emph{jump time} $\nu(y) \leq n_{ter}$ as follows.\smallskip\\
{\bf Case (i):} If there is a finite $k \leq  k_{ter}$ such that $z_k(y) \in STZ$, the smallest such $k$ will be the \emph{truncation time} $\kappa(y)$ . There is then a \emph{jump time} $\nu(y) \leq n_{ter}$ such that the one-step jump from the zero-cost  trajectory point $v_{\nu(y)}$ to $z_{\kappa(y)}(y))$ has cost at most $\eta$. \smallskip\\
{\bf Case (ii):} If there is no finite $k \leq k_{ter}$ with $z_k(y) \in STZ$, compute a finite $\kappa(y) \leq  k_{ter}$ and  $\nu(y) \leq n_{ter}$ by minimizing the cost of a one-step jump from the zero cost path $\mbold{v}$ to the reverse geodesic $\mbold{z}(y)$ so that
$$
(\nu(y), \kappa(y)) = \arg \min_{(n,k)} C(v_n, z_k(y)), 
$$ 
where the minimization is restricted to $n \leq n_{ter}$ and $k \leq k_{ter}$ with $k$ finite.  Truncate a portion of $\mbold{z}(y)$ by keeping only the $z_k = z_k(y)$ with $1 \leq k \leq \kappa(y)$ to define the \emph{truncated reverse geodesic} $RG(y)$ by
$$
RG(y) = [G\, y \, z_2 \ldots z_{\kappa(y)}]. 
$$

We say that $RG(y)$ is \emph{complete} in Case (i) and \emph{incomplete} in Case (ii).  By time reversion, each such truncated $RG(y)$ becomes a potential \emph{terminal geodesic} $TG(y)$ starting at $z_{\kappa(y)}(y)$ and ending at $G,$ which is defined for $1 \leq k \leq \kappa(y)$ by 
$$
TG_k(y) = z_{\kappa(y) +1 - k}(y).
$$
{\bf Broken Geodesics from $H$ to $G$:}  For each $y \in PEN$, define a \emph{broken geodesic} $BG(y)$ linking $H$ to $G$ in $(\nu(y) + \kappa(y))$ steps and having penultimate point $y$ by concatenating two geodesic segments as follows.
Let  $BG_n(y) = v_n$ for $1 \leq n \leq \nu(y)$ be the zero-cost initial geodesic segment of $BG(y)$.  Define the terminal geodesic segment of $BG(y)$ by shifting time in $TG(y)$ so that
$$
BG_{\nu(y) + k}(y) = z_{\kappa(y) +1 - k}(y)
$$
for $1 \leq k \leq \kappa(y).$  Compute then $\lambda_1 = \inf_{y \in PEN} \lambda(BG(y))$ with penultimate $y_1 \in PEN$ realizing this infimum. The value $\lambda_1$ is a first upper bound for the cost $\lambda(\mbold{h}^*) $ of a geodesic $\mbold{h}^*$  linking H to G, and $\mbold{h}^1 = BG(y_1)$ is a first-stage approximation of $\mbold{h}^*$ by a broken geodesic.
Define $PEN_2 \subset PEN $ as the set of all $y \in PEN$ such that $RG(y)$ is incomplete and $\lambda(TG(y)) > CAP_1$.  If $PEN_2$ is empty, we end the numerical search for the geodesic $\mbold{h}^*$ from $H$ to $G$, and we consider $\mbold{h}^1$ as the best approximation of $\mbold{h}^*$.  If $PEN_2$ is not empty, we launch a second stage in the search for $\mbold{h}^*$.
\subsection{\bf Multi-Stage Geodesic Search} \label{ss:multi}
Each histogram $G^{y} =z_{\kappa(y)}$ with $y \in PEN_2$ is now considered as a new target for geodesics starting at $H$. For each such $G^y$,  implement the first-stage algorithm by replacing $G$ by $G^{y}$. This yields a broken geodesic $\boldsymbol{\phi}^y$ linking $H$ to $G$ with penultimate $y \in PEN$. By concatenating $\boldsymbol{\phi}^y$ with $TG(y)$, we obtain a broken geodesic $\wh{BG}(y)$ with $3$ geodesic segments. Let $\lambda_2 = \min(\lambda(\wh{BG}y))$. If $\lambda_2 \geq \lambda_1,$ we end the geodesic search with unchanged output $\mbold{h}^2 = \mbold{h}^1$. If $\lambda_2 <  \lambda_1,$ the second stage outputs a broken geodesic $\mbold{h}^2$ with minimal cost $\lambda_2.$   Similarly to the end of the first stage, one can proceed to a third stage, and so on. However, in all our numerical experiments for $g = 3$ and $g = 4,$ the first stage was the only stage necessary to complete the geodesic search (see Conjectures \ref{conj:Discretization}--\ref{conj:Boundary}).




\subsection{\bf Improving Selection of Penultimate points} \label{ss:numerics}

In our multi-stage reverse shooting algorithm, the mesh size of the $\epsilon$-net $S$ and choice of a set $PEN \subset S$ of penultimate points are key factors in controlling computational complexity and accuracy of geodesic search. For $g = 3$ genotypes, we have launched intensive numerical explorations with $\epsilon = 0.002$ so that  $card(S) =125,751$ and $PEN = S - \partial S$ to exclude from $S$ all near-boundary histograms. Our numerical results for $g=3$  indicate that nearly optimal broken geodesics from $H$ to $G$ are typically generated at first stage.

However, this brute-force choice $PEN \simeq S$ becomes much heavier computationally for $g \geq 4$.  For instance, if $\epsilon = 0.005$ and $g= 4$, the $\epsilon$-net $S$  already has cardinal $1,373,701$. Clearly, letting $PEN \simeq S$ quickly becomes inefficient as $g$ increases and causes major memory issues.  

To improve computation times for our  multi-stage geodesic search, we have developed ``smarter'' algorithms for more efficient selection of the set $PEN.$ The case $g = 3$ where a brute-force choice of PEN is fully feasible provides a pragmatic context to numerically validate these algorithmic selection tools.  As we will see in the following conjectures and examples, the gradient of the one-step cost function $C(y,G)$ becomes a key factor in developing a reasonable set $PEN$.\smallskip \\
{\bf Cost Gradient:} Our numerical tests indicate that one-step costs $C(y,G)$ from many penultimate points $y \in \mathcal{H}^0$ are already too high to possibly be part of a broken geodesic linking $H$ to $G$. This suggests a  pruning indicator to  reduce the size of $PEN.$ Since for $m$ small enough, the gradient $\nabla_{y} C(y,G)$ can be explicitly approximated by \eqref{DyCyx1}--\eqref{DyCyx3}, we should discard from $PEN $ all points $y$ for which the norm  $h(y,G) = \vnorm{\nabla_{y} C(y,G)}$ is large and favor the inclusion of $y \in PEN$ for $h(y,G)$ small. This leads to a pragmatic conjecture below.\smallskip\\
{\bf Adaptive Mesh Size:} Multi-scale discretizations offer natural approaches to upgrade the efficiency of geodesic search. Our numerical exploration has indicated that an adaptive discretization with \emph{local  mesh size} based on the norm of cost gradients is an efficient selection tool to reduce $card(PEN).$ We formalize this as a conjecture.
\begin{conjecture}\label{conj:Discretization}
For any finite set $PEN$ of interior histograms and each $y \in PEN,$ define the \emph{Local Mesh Size} $LMS(y)$ as the minimum of $\vnorm{ x-y }$ over all $x$  in $PEN - \{ y \}$. Given a  target histogram $G \in \mathcal{H}^\circ,$ sparse but efficient finite sets $PEN$ of penultimate points should have local mesh size $LMS(y) \simeq c / h(y,G)$ where $h(y,G)$ is the norm of the cost gradient $\nabla_{y}C(y,G)$ and $c$ is a constant.
\end{conjecture}
\noindent{\bf Preliminary Test Paths:} A natural computational booster is to first generate a single  reverse geodesic $\mbold{z}(y_0)$ with $z_0(y_0) = G$ and  $z_1(y_0) = y_0.$ We then truncate it if and when it reaches the near boundary and link it to $H$ by a zero-cost path as outlined in Section \ref{ss:stage1}. This very fast computation provides a broken geodesic from $H$ to $G$ with cost $\lambda_0$, which we call a \emph{test path} below.  In the first-stage geodesic search, the iterative computation of any reverse geodesic $\mbold{z}(y)$ can then be stopped at step $k$ as soon as the cost of the terminal geodesic segment $[z_k(y) \ldots z_1(y)]$ is larger than $\lambda_0$. In any such case, the penultimate point $y$ can  be discarded from $PEN$. We now outline a good choice for $y_0,$ which can be found independent of the choice of $\epsilon$-net $S$ and $PEN$.


\begin{conjecture}\label{conj:PenPt}
Given a target $G \in \mathcal{H}^\circ$, any $y_0 \in \mathcal{H}^0$ which approximately minimizes the norm  $h(y,G)$ of the cost gradient provides the  penultimate point $y_0$ of an efficient preliminary test path. For small $m$, one can approximate $y_0$ by the solution $y^*$ of the following linear system,   
\begin{equation}\label{e:sys}
\begin{pmatrix}
\left( 1 - \frac{1}{G(1)} \right) F_1 & F_2 & F_3 &\dots & F_g \\
F_1 & \left( 1 - \frac{1}{G(2)} \right) F_2 & F_3 &\dots & F_g\\
\vdots& \vdots &\vdots & \vdots& \vdots\\
F_1& F_2& \dots & \left( 1 - \frac{1}{G(g-1)} \right) F_{g-1}& F_g\\
1& 1& 1& \dots &1 
\end{pmatrix}
\begin{pmatrix}
y(1)\\
y(2)\\
y(3)\\
\vdots\\
y(g-1)\\
y(g)
\end{pmatrix}
= 
\begin{pmatrix}
0\\
0\\
0\\
\vdots\\
0\\
1
\end{pmatrix}.
\end{equation}
Fix a constant $0 <c <1$. Define PEN(c) as the set of all $y \in \mathcal{H}^0$ such that $h(y,G) \leq h(y^*,G) / c$. For adequate choices of $c$, PEN(c) is an efficient  set of penultimate points.
\end{conjecture}
\noindent
{\bf Geodesics and boundary points:}
We have proved earlier that geodesics linking two interior histograms must be interior histograms paths. We have also numerically checked in many cases with  $N >>10^8$ , $m <<10^{-8}$,  and $b(G) \approx 100/N$  that the function $h(y,G)$ tends to be very large when $b(y) \to 0$. This leads to the following conjecture.
\begin{conjecture}\label{conj:Boundary}
For $N$ and $1/m$ large enough, nearly-minimal broken geodesics linking two interior histograms never reach the boundary of $\mathcal{H}$. For efficient selection of $PEN$, we expect the first stage of our geodesic search to generate nearly-minimizing geodesics linking two interior histograms.
\end{conjecture}







\subsection{Summary of Numerical Results}

We provide evidence for these conjectures in Section \ref{s:sims}. For Conjectures \ref{conj:Discretization} and \ref{conj:PenPt}, we have calculated numerous geodesics for various pairs of histograms $H$ and $G$ and analyzed these geodesics along with the norm $h(y,G)$ of the cost gradient and the one-step cost $C(y,G)$. We found that no nearly-optimal broken geodesic bounced off the boundary of $\mathcal{H}$ and that minimizing broken geodesics were generated from penultimate points very close to the minimizers of $h(y,G)$.  We also compared these selective generation of efficient sets $PEN$ to the brute force approach and shown  significant savings in computation time from the order of several minutes in the brute force approach to the order of a few seconds in the accelerated approach.

By partitioning $\mathcal{H}^\circ$ through quantiles of the numerical values taken by $1/h(y,G)$, we identified adequate ranges for the constants $c$ involved in Conjectures \ref{conj:Discretization} and \ref{conj:PenPt}. We used a $2.5\%$ quantile to generate nearly-optimal broken geodesics for various pairs $H$ and $G$. The range $10^{-5} \leq c \leq 10^{-2}$ was adequate when the approximate minimizer $y^*$ of $h(y,G)$ provided $h(y^*,G)$ is of order $m$.

To test Conjecture \ref{conj:Boundary} for $g = 3$ genotypes and three distinct sets of process parameters $\mcal{P}$,  we generated $1030$ nearly-optimal broken geodesics targeting  $1030$ distinct histograms. Our numerical computations showed that none of these $3090$  paths had to bounce off the boundary. The total CPU time on a current standard laptop was roughly $3$ hours, yielding an average CPU time of approximately $3.6$ sec. per cost minimizing path. To further solidify this conjecture, we numerically verified  that  the cost gradient norm $h(y,G)$ becomes very large when both target and penultimate points are near the boundary.

Towards the end of Section \ref{s:sims} and in Appendix \ref{app1} and \ref{app2}, we discuss further some attributes of the approximate minimizing paths for various target and initial histograms along with the associated trajectory costs $\lambda.$

