A multilevel adaptive reaction-splitting simulation method for stochastic reaction networks
Alvaro Moraes\thanks{Mathematical and Computer Sciences and Engineering Division
King Abdullah University of Science and Technology (KAUST)
Thuwal
Saudi Arabia ({\tt alvaro.moraesgutierrez@kaust.edu.sa
Stochastic modeling of reaction networks is a framework used to describe the time evolution of many natural and artificial systems, including, biochemical reactive systems at the molecular level, viral kinetics, the spread of epidemic diseases, and wireless communication networks, among many other examples. In this work, we present a novel multilevel Monte Carlo method for kinetic simulation of stochastic reaction networks that is specifically designed for systems in which the set of reaction channels can be adaptively partitioned into two subsets characterized by either ``high'' or ``low'' activity. Adaptive in this context means that the partition evolves in time according to the states visited by the stochastic paths of the system. To estimate expected values of observables of the system at a prescribed final time, our method bounds the global computational error to be below a prescribed tolerance, $TOL$, within a given confidence level. This is achieved with a computational complexity of order $\Ordo{TOL^{-2$, the same as with an exact method, but with a smaller constant. We also present a novel control variate technique based on the stochastic time change representation by Kurtz, which may dramatically reduce the variance of the coarsest level at a negligible computational cost. Our numerical examples show substantial gains with respect to the standard Stochastic Simulation Algorithm (SSA) by Gillespie and also our previous hybrid Chernoff tau-leap method.
main_mixed.tex
% flatex input: [main_mixed.tex]
%\documentclass{amsart}

\documentclass[final,leqno]{siamltex}
%\documentclass[final,leqno]{siamltex1213}

% Packages
\usepackage{amsmath,amssymb}
\usepackage{mathtools}
%\usepackage{epstopdf}
%\usepackage{subfigure}
%\usepackage{makecell}
%\usepackage{multicol}
%\usepackage{enumerate}
%\usepackage{float}
\usepackage{hyperref} % must be last package
\hypersetup{
 %draft = true, % turns off all (for printing)
 %bookmarks=true,
 colorlinks=true,
 linktocpage=true,
 bookmarksnumbered=true
}
\usepackage{algorithm,algorithmic}
%\usepackage[all]{xy}
\usepackage[pdftex]{graphicx,color}
%\usepackage{marginnote}

% Theorems
\newtheorem{thm}{Theorem}[section]
%\newtheorem{proof}[thm]{Proof}
%\newtheorem{lemma}[thm]{Lemma}
\newtheorem{rem}[thm]{Remark}
\newtheorem{cor}[thm]{Corollary}
\newtheorem{ex}[thm]{Example}
\newtheorem{assu}[thm]{Assumption}
\newtheorem{alg}[thm]{Algorithm}
\newtheorem{defn}[thm]{Definition}

% Commands
\newcommand{\ceil}[1]{\lceil {#1} \rceil}
\newcommand{\seqof}[3]{(#1)_{#2}^{#3}}

\newcommand{\ssa}{\text{SSA}}
\newcommand{\hyb}{\text{Hyb}}
\newcommand{\MNRM}{\text{MNRM}}
\newcommand{\TL}{\text{TL}}
\newcommand{\ML}{\text{ML}}

%\newcommand{\NSSA}{N_{\ssa}}
\newcommand{\NSSAKone}{N_{\ssa,K1}}
\newcommand{\NSSAKtwo}{N_{\ssa,K1'}}

%\newcommand{\NMNRMKone}{N_{\MNRM,K1}}
%\newcommand{\NMNRMKtwo}{N_{\MNRM,K2}}
%\newcommand{\NMNRMKoneC}{N_{\MNRM,K1}^{(c)}}
%\newcommand{\NMNRMKtwoC}{N_{\MNRM,K2}^{(c)}}

\newcommand{\NMNRMKone}{N_{K1}}
\newcommand{\NMNRMKtwo}{N_{K2}}
\newcommand{\NMNRMKoneC}{N_{K1}^{(c)}}
\newcommand{\NMNRMKtwoC}{N_{K2}^{(c)}}

\newcommand{\NMNRM}{N_{\MNRM}}
\newcommand{\NMNRMC}{N_{\MNRM}^{(c)}}
%\newcommand{\NMNRMP}{N_{{\MNRM}^*}}

\newcommand{\NSSAP}{N_{{\ssa}^*}}
\newcommand{\NTL}{N_{\TL}}
\newcommand{\NTLP}{N_{{\TL}^*}}
\newcommand{\NTLC}{N_{\TL}^{(c)}}
%\newcommand{\NTLF}{\dbar N_{\TL}}

\newcommand{\EE}{\mathcal{E}}

\newcommand{\WE}{\EE_I}
\newcommand{\WEC}[1]{\EE_{I,#1}}
\newcommand{\WEH}[1]{\hat \EE_{I,#1}}
\newcommand{\WSSA}{\text{Cost}_{\ssa}}
\newcommand{\WHYB}{\text{Cost}_{\hyb}}


\newcommand{\avg}[2]{\mathcal{A}\left(#1;#2\right)}
%\newcommand{\SD}[2]{\mathcal{S}(#1;#2)}
\newcommand{\svar}[2]{\mathcal{S}^2\left(#1;#2\right)}

\newcommand{\dbar}[1]{\Bar{\Bar{#1}}}

\newcommand{\ie}{\emph{i.e.}}
\newcommand{\eg}{\emph{e.g.}}
\newcommand{\cf}{\emph{cf.}}
\newcommand{\ud}{\mathrm{d}}
\newcommand{\prob}[1]{\mathrm{P}\left(#1\right)}
\newcommand{\expt}[1]{\mathrm{E}\left[#1\right]}
\newcommand{\expth}[1]{\hat{\mathrm{E}}\left[#1\right]}

\newcommand{\var}[1]{\mathrm{Var}\left[#1\right]}
\newcommand{\cov}[2]{\mathrm{Cov}\left[#1, #2\right]}
\newcommand{\norm}[1]{\left\|#1\right\|}
\newcommand{\abs}[1]{\left|#1\right|}
\newcommand{\rset}{\mathbb{R}}
\newcommand{\nset}{\mathbb{N}}
\newcommand{\zset}{\mathbb{Z}}
\newcommand{\JAC}{\mathbb{J}}
\newcommand{\one}{\mathbf{1}}
\newcommand{\poisson}[1]{\mathcal{P}\left(#1\right)}
\newcommand{\normal}[1]{\mathcal{N}\left(#1\right)}
\newcommand{\binomial}[1]{\mathcal{B}\left(#1\right)}
\newcommand{\Ordo}[1]{{\mathcal{O}}\left(#1\right)}
\newcommand{\OrdoP}[1]{{\mathcal{O}_P}\left(#1\right)}
\newcommand{\ordo}[1]{{o}\left(#1\right)}
\newcommand{\abstr}[1]{\begin{abstract}#1\end{abstract}}
\newcommand{\thx}[1]{\thanks{#1}}
\newcommand{\red}[1]{{\color{red}#1}}
\newcommand{\blue}[1]{{\color{blue}#1}}
\newcommand{\green}[1]{{\color{magenta}#1}}

\newcommand{\latt}{\mbox{$\zset_+^d$}}


%Alvaro's notations

\newcommand{\indicator}[1]{\mathbf{1}_{#1}} 

\newcommand{\PERIOD}{.}
\newcommand{\COMMA}{,}

% Removed this command since it does not work properly with spacing
%\newcommand{\ChB}{{Chernoff bound }}
\newcommand{\LP}{\left(}
\newcommand{\RP}{\right)}

\newcommand{\LPe}{\left.}
\newcommand{\RPe}{\right.}

%\newcommand{\SEP}{\biggm|}
\newcommand{\SEP}{\, \big| \,}
\newcommand{\BX}{\bar X(t)}
\newcommand{\BXi}{\bar X_i(t)}
\newcommand{\BXitau}{\bar X_i(t{+}\tau)}

\newcommand{\sj}{\sum_{j=1}^J}

\newcommand{\Qi}{Q_i(t,\tau_i)}
\newcommand{\lji}{a_j(\bar X(t))\tau_i}
\newcommand{\Yj}{Y_j\left( a_j(\bar X(t))\tau_i \right)}
\newcommand{\ldi}{\log(\delta_i)}
\newcommand{\taui}{\tau_i(s)}
\newcommand{\DEN}[1]{- a_0 \LP \bar X(t) \RP +  \sj a_j(\bar X(t))e^{-#1\, \nu_{ji}}}
\newcommand{\aj}{ a_j(\bar X(t))}

\newcommand{\anu}{\sj \aj \nu_{ji}}
%\newcommand{\anu2}{\sj \aj (\nu_{ji})^2}
\newcommand{\tsi}{\tilde{s}_i}
\newcommand{\hs}{\hat{s}}
%\newcommand{\tti1}{\tilde{t}_{i,1}} 
\newcommand{\az}{a_0(\bar X(t))}
\newcommand{\dnux}{\delta_i^{\,\nu_{ji}/{\BXi}}}

\newcommand{\bj}{b_{ji}(\bar X(t))}
\newcommand{\mE}[1]{\mathcal{E}_{#1}}
\newcommand{\mV}[1]{\mathcal{V}_{#1}}
\newcommand{\hmV}[1]{\hat{\mathcal{V}}_{#1}}

\newcommand{\orden}{\alpha}

%Pedro's notations
\newcommand{\X}[1]{\bar X(t_{#1})}
\newcommand{\ti}[1]{t_{#1}}
%\newcommand{\iff}{\Leftrightarrow}

% Document info
%\title{Mixed Chernoff Tau-leap}
\title{A multilevel adaptive reaction-splitting simulation method for stochastic reaction networks}

\author{Alvaro Moraes\thanks{Mathematical and Computer Sciences and Engineering Division,
 King Abdullah University of Science and Technology (KAUST),
 Thuwal, Saudi Arabia ({\tt alvaro.moraesgutierrez@kaust.edu.sa}).}
        \and Raul Tempone\thanks{Mathematical and Computer Sciences and Engineering Division,
 King Abdullah University of Science and Technology (KAUST),
 Thuwal, Saudi Arabia ({\tt raul.tempone@kaust.edu.sa}).}
\and Pedro Vilanova\thanks{Mathematical and Computer Sciences and Engineering Division,
 King Abdullah University of Science and Technology (KAUST),
 Thuwal, Saudi Arabia ({\tt pedro.guerra@kaust.edu.sa}).}}
 
%\date{} 

% 60J75, Jump processes
% 60H10, Stochastic ordinary differential equations 
% 65C30, Stochastic differential and integral equations
% 60H35, Computational methods for stochastic equations
% 65Y20, Complexity and performance of numerical algorithms
% 65L50, Mesh generation and refinement
% 92C40, Biochemistry, molecular biology

\numberwithin{equation}{section}
\numberwithin{figure}{section}
\numberwithin{table}{section}

\begin{document}
\maketitle
%\green{  
%\begin{enumerate}
%\item{Quick intro}
%\item{What is done here and what is novel}
%\item{Numerical examples to show how powerful are our techniques}
%\end{enumerate}}
\begin{abstract}
Stochastic modeling of reaction networks is a framework used to describe the time evolution of many natural and artificial systems, including, biochemical reactive systems at the molecular level, viral kinetics, the spread of epidemic diseases, and wireless communication networks, among many other examples.
In this work, we present a novel multilevel Monte Carlo method for kinetic simulation of stochastic reaction networks  that is specifically designed for systems in which the set of reaction channels can be adaptively partitioned into two subsets characterized by either ``high'' or ``low'' activity. 
Adaptive in this context means that the partition evolves in time according to the states visited by the stochastic paths of the system.
To estimate expected values of observables of the system at a prescribed final time, our method bounds the global computational error to be below a prescribed tolerance, $TOL$, within a given confidence level. This is achieved with a  computational complexity of order $\Ordo{TOL^{-2}}$, the same as with an exact method, but with a smaller constant. We also present a novel control variate technique based on the stochastic time change representation by Kurtz, which may dramatically reduce the variance of the coarsest level
at a negligible computational cost.
Our numerical examples show substantial gains with respect to the standard Stochastic Simulation Algorithm (SSA) by Gillespie and also our previous hybrid Chernoff tau-leap method. 
\end{abstract}


\begin{keywords}
Error estimates, error control, control variates, weak approximation, hybrid algorithms, multilevel Monte Carlo, Chernoff tau-leap, reaction splitting
\end{keywords}

\begin{AMS}
60J75, 60J27, 65G20, 92C40 
\end{AMS}

\setcounter{tocdepth}{1}
%\tableofcontents

\section{Introduction}
Stochastic reaction networks (SRN) are mathematical models that employ Markovian dynamics to describe the {time evolution of interacting particle systems} where one particle interact with the others through a finite set of reaction channels.
Typically, there is a finite number of interacting chemical species $(S_1,S_2,\ldots,S_d)$ and a stochastic process, $X$, such that its $i$-th coordinate is a non-negative integer number $X_i(t)$ that keeps track of the abundance of the $i$-th species at time $t$.  Therefore, the state space of the process $X$ is the lattice $\latt$.
%Among the {phenomena usually modeled by SRN  we mention: biochemical reactions at thermal equilibrium and constant volume, spread of epidemic diseases in homogeneously mixed populations, natural and artificial neural networks, pharmacokinetics, virus kinetics and dynamics of social networks. This work presents a novel multilevel Monte Carlo Chernoff tau-leap method  based on our previous work \cite{ourSL,ourML}.

Our main goal is to estimate the expected value $\expt{g(X(T))}$, where $X$ is a non-homogeneous Poisson process describing a SRN, and $g:\rset^d\to\rset$ is a given real observable of $X$ at a final time $T$. 
%The process $X$ is a continuous-time Markov chain {defined through a finite number of reaction channels and its corresponding propensity functions},  taking values in the lattice of non-negative integers, $\zset_+^d$.  
Pathwise realizations can be simulated
exactly using the {Stochastic Simulation Algorithm} (SSA),
introduced by Gillespie in \cite{gillespie_ssa} (also known as Kinetic Monte Carlo among physicists, see \cite{KMC} and references therein), or the {Modified Next Reaction Method} (MNRM) introduced by Anderson in \cite{anderson2007modified}, among other methods. Although these algorithms generate exact realizations of $X$, they may be  computationally expensive for systems that undergo high activity. For that reason, Gillespie proposed in
\cite{gillespie_tau_leap} the tau-leap method to approximate the SSA by evolving the process with fixed time
steps while freezing the propensity functions at the beginning of each time step. %In fact, the tau-leap method can be seen as a forward Euler method for a stochastic differential equation driven by Poisson random measures (see, \eg, \cite{li}).

A drawback of the tau-leap method is that the simulated paths may take negative values, which is a nonphysical consequence of the approximation and not a qualitative feature of the original process. For that reason, we proposed in \cite{ourSL,ourML} a Chernoff-based hybrid method that switches adaptively between the tau-leap and an exact method. This allows us to control the probability of reaching negative values while keeping the computational work substantially smaller than the work of an exact method.
%\blue{TODO: Advantages, sell hybrid.}
%Let us consider the sum of the propensity functions as a measure of the total activity of the process $X$ at a state $x$. 
The hybrid method developed in \cite{ourSL,ourML} can be successfully applied  to systems where the state space, $\zset_+^d$, can be decomposed into two regions   according to the activity of the system; where all the propensities are uniformly low or uniformly high, \ie, non-stiff systems.
To handle stiff systems, we first measure the total activity of the system at a certain state by the total sum of the propensity functions evaluated at this state. 
The activity of the system is low when all the propensities are uniformly low, but a high level of activity can be the result of a high activity level in one single channel.   
This observation suggests that to reduce computational costs, we should adaptively split the set of reaction channels into two subsets according to the individual high and low activity levels.  
It is natural to evolve the system in time by applying the tau-leap method to the high activity channels and an exact method to the low activity ones. 
This is the main idea we develop in this work.

Reaction-splitting methods for simulating stochastic reaction networks are treated for instance in \cite{Harris2006, Puchalka2004, Haseltine2002, Plyasunov2005}, but our work is, to the best of our knowledge, the first that  
i) achieves the computational complexity of an exact method like the SSA by using the multilevel Monte Carlo paradigm, 
ii) explicitly uses a  decomposition of the global error to provide all the simulation parameters needed to achieve our goal with minimal computational effort,
iii) effectively controls the global probability of reaching negative populations with the tau-leap method, and 
iv) needs only two user-defined parameters that are natural quantities - the maximum allowed relative global error or tolerance and the confidence level.  

In \cite{Harris2006}, the authors propose an adaptive reaction-splitting scheme that considers not only the exact and tau-leap methods but also the Langevin and Mean Field ones. 
Their main goal is to obtain fast hybrid simulated paths, and they do not try to control the global error. 
The efficiency of their method is measured a posteriori using smoothed frequency histograms that should be close to the exact ones according to the distance defined by Cao and Petzold in \cite{CaoPetzold2006}. %This measure of fit is not incorporated into the process. %(in the sense that running parameters are inputs of the algorithms)
In their work, the tau-leap step is chosen according to the ``leap condition'' (as in \cite{Cao2006}) but they do not perform a rigorous control of the global discretization error.
In order to avoid negative populations, the authors reverse population updates if any value is found to be negative after accounting for all the reactions. Then, the tau-lep step size is decremented and the path simulation is restarted. 
This approach introduces bias in the estimations, and even by controlling the small reactant populations, a tau-leap step always may lead to negative populations subsequently increasing its computational work. Our Chernoff -based bound is a fast and accurate procedure to obtain the correct tau-leap step size. 
Finally, the method in \cite{Harris2006} needs to define three parameters that quantify the speed of the reaction channels, which, in principle, are not trivial to determine for a given problem.

Puchalka and Kierzek's approach \cite{Puchalka2004} seems to be closest to our approach in spirit since they also explore the idea of adaptively splitting the set of reaction channels using the tau-leap method for the fast ones and an exact method for the slow ones. 
They seek to simulate fast approximate paths while maintaining qualitative features of the system. 
The quantitative features are checked a posteriori against an exact method.  
Regarding their tau-leap step size selection, Puchalka and Kierzek consider a user-defined maximal time step empirically chosen by numerical tests instead of controlling the discretization error.
Their classification rule is applied individually to each reaction channel.  It takes into account both the percentage of individual activity and the abundance of the species consumed. In a certain sense it can be seen as a way of controlling the probability of negative populations and an ad-hoc manner to split the reaction channels by optimizing the computational work.
%in this work, we want to approximate the expected value of an observable of the process at a final time $T$ within a prescribed global error (tolerance) up to a given confidence level. In \cite{Puchalka2004}, . 
%Another difference is that for determining the tau-leap step size, we use the aforementioned Chernoff bound that globally controls the probability of negative populations while the authors 
%Regarding the splitting strategy: in our case we split the reaction channels by optimizing the computational work of reaching the next point in the time grid, while 

In \cite{Haseltine2002} and \cite{Plyasunov2005}, the reaction-splitting issue is addressed but the partition method is not adaptive, \ie, fast and slow reaction channels are identified offline and are inputs of the algorithms.
We note that these works do not provide any measure or control of the resulting global error. %  in their algorithms.
 Furthermore, they do not control the probability of attaining negative populations.   
%In order to do that, we first develop a computational work based decision rule that chooses a suitable partition at each time step, in terms of computational work. Then, using our previous results, we develop an algorithm which generates mixed paths, that is, it applies the tau-leap method and the modified next reaction method at the same time.

In the remaining of this section, we introduce the mathematical model and the path simulation techniques used in this work.
In Section \ref{sec:mixedpaths}, we present an algorithm to generate mixed trajectories; that is, the algorithm generates a trajectory using an exact method for the low activity channels and the Chernoff tau-leap method for the high activity ones. Then, inspired by the ideas of Anderson and Higham \cite{Anderson2012}, we propose an algorithm for coupling two mixed Chernoff tau-leap paths. This algorithm uses four building blocks that result from the combination of the MNRM and the tau-leap methods.
In Section \ref{sec:mlmcerr}, we propose a mixed MLMC estimator. Next, we introduce a global error decomposition and show that the computational complexity of our method is of order $\Ordo{TOL^{-2}}$.
Finally, we show the automatic procedure that estimates our quantity of interest within a given prescribed relative tolerance, up to a given confidence level.
Next, in Section \ref{sec:cvar}, we present a novel control variate technique to reduce the variance of the quantity of interest at level 0. 
%and finally, we indicate a Bayesian strategy for estimating parameters involved in the computational work optimization problem. 
In Section \ref{sec:examples}, the numerical examples illustrate the advantages of the mixed MLMC method over the hybrid MLMC method presented in \cite{ourML} and to the SSA. Finally, Section \ref{sec:conclusions} presents our conclusions.

\subsection{A Class of Markovian Pure Jump Processes}
\label{sec:pjp}
In this section, we describe the class of Markovian pure jump processes, $X:[0,T]\times \Omega \to\zset_+^d$, frequently used for modeling stochastic biochemical reaction networks.  
 
Consider a biochemical system of $d$ species interacting
through $J$ different reaction channels. For the sake of brevity, we write $X(t,\omega) {\equiv} X(t)$.  Let $X_i(t)$ be the
number of particles of species $i$ in the system at time $t$. We 
study the evolution of the state vector, $X(t) = (X_1(t), \ldots, X_d(t)) \in
  \zset_+^d$,
modeled as a continuous-time Markov chain starting at $X(0) \in \zset_+^d$. Each reaction can be described by the
vector $\nu_j\in \zset^d$, such that, for a state vector $x\in
\zset_+^d$, a single firing of reaction $j$ leads to the change
$x \to x+ \nu_j$.
The probability that reaction $j$ will occur during the small interval
$(t,t{+}\ud t)$ is then assumed to be
\begin{equation}\label{eq:prob}
  \prob{X(t+\ud t)=x+\nu_j | X(t) = x} =  a_j(x) \ud t + \ordo{\ud t}
\end{equation}
for a given non-negative {polynomial propensity function}, $a_j:\rset^d \to
\rset$. We set $a_j(x){=}0$ for those $x$ such that $x{+}\nu_j\notin\latt$.
%A process, $X$, that satisfies \eqref{eq:prob}, is a continuous-time, discrete-space
%Markov chain that 
The process $X$ admits the following random time change representation by Kurtz \cite{kurtzmp}:
\begin{equation}
  \label{eq:exact_process}
  X(t) = X(0) + \sum_{j=1}^J  \nu_j Y_j \LP \int_0^t a_j(X(s))\, \ud s \RP  \COMMA
\end{equation}
where $Y_j:\rset_+ {\times} \Omega \to\zset_+$ are independent unit-rate Poisson processes.
Hence, $X$ is a non-homogeneous Poisson process.

%\begin{rem}
%In this setting, the solution of the following system of ordinary differential equations
%\begin{align}\label{eq:meanfield}
%\left\{ \!
%\begin{array}{lll}
%\dot x (t)  &=&  \nu a(x(t)), \,\, t \in \rset_+ \\
%x(0) &=& x_0 \in \rset_+ \PERIOD
%\end{array}
%\right.
%\end{align}
%is called the Mean Field of the process $X$, where $\nu$ is the matrix with columns $\nu_j$ and $a(x)$ is the column vector of propensities. In this work, we use the Mean Field in Section \ref{sec:cvar}.%path for scaling and preprocessing constants. %associated with the computational work of the SSA and Chernoff tau-leap steps.
%\end{rem}
\subsection{The Modified Next Reaction Method (MNRM)}\label{sec:MNR}
%\subsection{Description of the Modified Next Reaction Method}
The MNRM, introduced in \cite{anderson2007modified}, and based on the Next
Reaction Method  \cite{GibBruck}, is an exact simulation algorithm like Gillespie's SSA that explicitly uses  representation \eqref{eq:exact_process} for simulating exact paths and generates only one exponential random variable per iteration. The reaction times are modeled with firing times of Poisson processes, $Y_j$, with
internal times given by the integrated propensity functions.
%Such a representation allows to consider the system as a whole, avoiding the
%computation of the next reaction and its time. 
The randomness is now
separated from the state of the system and is encapsulated in the $Y_j$'s.
%For each reaction, $j$, the internal time %in \eqref{eq:internal_time} 
%is crucial for
%understanding this method. % ,
%is defined as $R_j(t) {=} \int_0^t a_j(X(s))ds$. There are $J{+}1$ time frames in the system, the absolute one, $t$, and one for each Poisson
%process, $Y_j$.
Computing the next reaction and its time is equivalent to computing how
much time passes before one of the Poisson process, $Y_j$, fires, and
which process fires at that particular time, by taking the minimum of such times. %TheNRM and MNRM make use of the internal times to reduce the number of simulated random variables by half.  
%In the following, we describe the MRNM and then we present its implementation in Algorithm \ref{alg:mnr}.

%Given $t$, we have the propensity $a_j{=}a_j(X(t))$ and the internal
%time $R_j{=}R_j(t)$. Now, let $\Delta R_j$ be the remaining
%time for the reaction, $j$, to fire, assuming that $a_j$ stays constant
%over the interval $[t,t{+}\Delta R_j)$. Then, $t{+}\Delta R_j$ is the time
%when the next reaction, $j$, occurs. The next internal time at which the
%reaction, $j$, fires is then given by $R_j {+} a_j \Delta R_j$. When
%simulating the next step, the first reaction that fires occurs after
%$\Delta {=} \min_j \Delta R_j$. We then update the state of the system according to
%that reaction, add $\Delta$ to the global time, $t$, and then update the internal times by
%adding $a_j \Delta$ to each $R_j$.  We are left to
%determine the value of $\Delta R_j$, \ie, the amount of time until the
%Poisson process, $Y_j$, fires, taking into account that $a_j$ remains
%constant until the first reaction occurs. Denote by $R_j$ the first
%firing time of $Y_j$ that is strictly larger than $R_j$, \ie,
%$P_j {:=} \min  \{s{>}R_j:Y_j(s){>}Y_j(R_j)\}$ and finally $\Delta R_j {=} \frac{1}{a_j}(P_j-R_j)$.
%where $P_j$ are exponentially distributed random variables that are used to update the internal times.

%\begin{algorithm}
%\caption{The Modified Next Reaction Method. Inputs: the initial state, $X(0)$, the next grid point, $T_0> t_0$, the propensity functions, $\seqof{a_j}{j=1}{J}$, the stoichiometric vectors, $\seqof{\nu_j}{j=1}{J}$.
%Outputs: the history of system states, $\seqof{X(t_k)}{k=0}{K}$. Here, we denote $S\equiv \seqof{S_j}{j=1}{J}$, $P\equiv \seqof{P_j}{j=1}{J}$, and $R\equiv \seqof{R_j}{j=1}{J}$.
%}
%\label{alg:mnr}
%\begin{algorithmic}[1]
%\STATE $k \leftarrow 0$, $t_k \leftarrow 0$, $X(t_k) \leftarrow X
%(0)$ and $R \leftarrow \bf{0}$
%\STATE Generate $J$ independent, uniform$(0,1)$ random numbers, $r_j$
%\STATE $P \leftarrow \seqof{\log(1/r_j)}{j=1}{J}$
%\WHILE {$t_k<T_0$}
%\STATE $S \leftarrow \seqof{a_j(X(t_k))}{j=1}{J}$
%\STATE $\seqof{\Delta R_j}{j=1}{J} \leftarrow \seqof{(P_j - R_j)/S_j}{j=1}{J}$
%\STATE ${\mu} \leftarrow \text{argmin}_j \{\Delta R_j\}$
%\STATE $\Delta \leftarrow \min_j\{ \Delta R_j \}$ 
%\STATE $t_{k+1} \leftarrow t_k+\Delta$
%\STATE $X(t_{k+1}) \leftarrow X(t_k) + \nu_\mu$
%\STATE $R \leftarrow R + S \Delta$
%\STATE $r \leftarrow$ uniform$(0,1)$
%\STATE $P_{\mu} \leftarrow P_{\mu} + \log(1/r)$
%\STATE $k \leftarrow k{+}1$
%\ENDWHILE
%\RETURN $\seqof{X(t_l)}{l=0}{k-1}$
%\end{algorithmic}
%\end{algorithm}

%Here $R_j=R_j(t)$ are internal times, and $\Delta R_j$ is the amount of time for the $j$-th reaction to fire, assuming $a_j$ is constant over $[t,t{+}\Delta R_j)$. The internal
%time at which reaction $j$ fires is $R_j+ a_j \Delta R_j$. 

%Among the advantages already mentioned, we can easily modify Algorithm 1 to generate paths in the cases where the rate functions depend on time and also when there are reactions delayed in time.
It is important to mention that the MNRM is used to simulate correlated exact/tau-leap paths as well as nested  tau-leap/tau-leap paths, as in \cite{ourML,Anderson2012}.
%In \cite{ourML}, this technique is used to develop a uniform step, unbiased multilevel Monte Carlo (MLMC) algorithm. 
In Section \ref{cmp}, we use this feature for coupling two mixed paths.


\subsection{The Tau-Leap Approximation}\label{sec:num_approx}

%\red{TODO: Here we can add the Kurtz representation.}
In this section, we define $\bar X$, the tau-leap approximation of the process, $X$, which follows from applying the forward-Euler approximation to the integral term in the random time change representation \eqref{eq:exact_process}.

The tau-leap method was proposed
in \cite{gillespie_tau_leap} to avoid the computational drawback of the exact methods, \ie,  when many reactions occur
during a short time interval. The tau-leap process, $\bar X$, starts from $X(0)$ at time $0$, and 
given that $\bar X(t) {=} \bar x$ and a
time step $\tau{>}0$, we have that
$\bar X$ at time $t{+}\tau$ is generated by
\begin{equation*}
  %\label{eq:tau_leap}
  \bar X(t+\tau) = \bar x + \sum_{j=1}^J 
  \nu_j \mathcal{P}_j\left( a_j(\bar x)\tau \right) \COMMA
\end{equation*}
where $\{\mathcal{P}_j(\lambda_j)\}_{j=1}^J$ are independent
Poisson distributed random variables with parameter $\lambda_j$,
used to model the number of times that the reaction $j$ fires
during the $(t,t{+}\tau)$ interval.  Again, this is nothing else than a 
forward-Euler discretization of the stochastic differential equation formulation of the pure jump process
\eqref{eq:exact_process}, realized by the Poisson random measure with
state-dependent intensity (see, \eg, \cite{li}).

In the limit, when $\tau$ tends to zero, the tau-leap method gives the same solution
as the exact methods \cite{li}. The total number of firings in each channel is 
a Poisson distributed stochastic variable depending only on the
initial population, $\bar X(t)$.  The error thus comes from the
variation of $a(X(s))$ for $s\in(t,t{+}\tau)$. 
%We observe that the computational work of a tau-leap step involves the generation of $J$ independent Poisson random variables. This is in contrast to the computational work of an exact step, which only involves the work of generating two uniform random variables, in the case of the SSA, and only one in the case of MNRM. 
%In \cite{kt}, the authors assume that there exists a vector, $w\in\rset_+^d$, such that $(w,\nu_j)\leq 0 $, for any reaction $\nu_j$.
%Therefore, every reaction, $\nu_j$, must have at least one negative component. This means that the species can be either transformed into other species or be consumed during the reaction. As a consequence, the space of states is contained in a simplex with vertices in the coordinate axis. This assumption excludes, for instance,  birth processes.
%In our numerical examples, we allow the set of possible states of the system to be infinite, but we explicitly avoid cases in which one or more species grows exponentially fast or blows up in the time interval $[0,T]$.

\subsection{The Hybrid Chernoff Tau-leap Method}
\label{sec:chernoff}
In \cite{ourSL}, we derived a Chernoff-type bound that allows us to guarantee that
the one-step exit probability in the tau-leap method is less than a
predefined quantity, $\delta{>}0$. 
The idea is to find the largest possible time step, 
$\tau$, such that, with high probability, in the next step, the approximate process, $\bar{X}$, will take a value in the lattice, $\latt$, of non-negative integers.
%, \ie,
%\begin{equation}
%  \label{eq:taucond}
%  \prob{ \BX + \sj \nu_{j} Y_j\left( a_j(\BX)\tau \right) \in \latt \Biggm| \BX}\geq 1- \delta,
%\end{equation}
%for some small $\delta>0$. Observe that this value of $\tau$ depends
%on $\BX$. %\Biggm|
%
This can be achieved by solving $d$ auxiliary
problems, one for each $x$-coordinate, $\bar X_i(t)$,
$i=1,2,\ldots,d$ as follows. Find the largest possible $\tau_i\geq 0$, such that
\begin{equation}
\label{eq:taui-deltai}
\prob{  \BXi + \sj \nu_{ji} \mathcal{P}_j\LP a_j\LP \bar X(t)\RP\tau_i\RP < 0 \Biggm| \BX}\leq \delta_i \COMMA
\end{equation}
where $\delta_i {=} \delta/d$, and $\nu_{ji}$ is the $i$-th coordinate of
the $j$-th reaction channel, $\nu_j$. Finally, we let $\tau{:=}\min\{\tau_i: i=1,2,\ldots,d\} \PERIOD$
%We refer to \cite{ourSL} for details on
%In order to find the largest time steps, $\tau_i$, let $\Qi {:=} \sj (-\nu_{ji})\mathcal{P}_j\LP a_j\LP \bar X(t)\RP\tau_i\RP$. Then, for all $s{>}0$, 
%using the Markov inequality, we obtain an upper bound for the probabilities we want to control:
%\begin{equation}
%\label{eq:Markov}
%\begin{aligned}
%  \prob{ \Qi>\BXi \SEP \BX} &= \prob{ \exp{\LP s \Qi \RP} >  \exp{\LP s \BXi \RP} \Bigm| \BX}  \\
%  &\leq \dfrac{\expt{\exp\LP s\Qi \RP}}{\exp\LP s \BXi \RP} \PERIOD
%\end{aligned}
%\end{equation}
%Observe that the independent Poisson random variables, $\Yj$, have
%moment-generating functions,
%\begin{equation*}
%M_j(s) = \exp\LP \lji (e^s-1)\RP \COMMA
%\end{equation*}
%and, therefore,
%\begin{equation}
%  \label{eq:MGF}
%  \begin{aligned}
%    \expt{\exp\LP s\Qi \RP}&= \prod_{j=1}^J M_j(-s \nu_{ji})\\
%    &= \exp\LP \tau_i \sj \aj (e^{-s \nu_{ji}}-1) \RP \PERIOD
%  \end{aligned}
%\end{equation}
%
%By combining \eqref{eq:Markov} and \eqref{eq:MGF}, 
%we have the Chernoff bound: % for the multi-reaction case,
%\begin{equation*}
%\prob{ \Qi{>}\BXi \SEP  \! \BX} \leq \inf_{s>0}  \exp \! \LP \! \! -s \BXi + \tau_i \sj \aj (e^{-s \nu_{ji}}{-}1) \! \RP \! .
%\label{eq:chboundmulti}
%\end{equation*}
%To avoid the computational problem of finding exactly the above infimum and to
%guarantee that
%\begin{equation*}
%\prob{ \Qi>\BXi \SEP \BX} \leq \delta_i \COMMA
%\end{equation*}
%we proceed as follows.
%First, according to \eqref{eq:taui-deltai} and \eqref{eq:chboundmulti},
%\begin{equation*}
%-s \BXi + \tau_i \sj \aj (e^{-s \nu_{ji}}-1) = \ldi \PERIOD
%\end{equation*}
%Using this fact, we can 
%Expressing $\tau_i$ as a function of $s$, we write
%\begin{equation*}
%\label{eq:taui}
%\taui = \dfrac{\ldi + s \BXi}{\DEN{s}}=:\dfrac{R_i(s)}{D_i(s)}\COMMA%{\underbrace{\DEN{s}}_{=:D_i(s)}},
%\end{equation*}
%where 
%\begin{equation*}
%  \az := \sj \aj \PERIOD
%\end{equation*}
%\subsubsection{Study of $\tau_i(s)$} \label{sec:taui}
%In this section, 
%We want to maximize $\tau_i$ while satisfying condition \eqref{eq:taui-deltai}. 
%Obviously, it is satisfied for $\tau_i=0^+$. By a continuity argument,  
%Let $\tau_i^*$ be this maximum.
%  such that every point of the interval $[0,\tau_i]$ satisfies \eqref{eq:taui-deltai}. %Note that $\tau_i^*$
%could be $+\infty$. 
%
%We discuss how depending on certain relations among the pairs
%$\{(\aj,\nu_{ji})\}_{j=1}^J$, we can conclude that $\tau_i^*$
%is either a real number or $+\infty$. First of all, if $\nu_{ji}>0$, $\forall
%j$, then $\tau_i^*$ must be $+\infty$, since no reaction is pointing
%to zero. From now on, we assume that, given the coordinate $i$, there
%is at least one reaction pointing to zero, \ie,
%\begin{equation}
%\exists j \mbox{ such that } \nu_{ji}<0 \label{eq:as1} \PERIOD
%\end{equation}
%
%
%The denominator of \eqref{eq:taui} is the function, 
%\begin{equation}
%\label{eq:DEN}
%D_i(s) := \DEN{s},
%\end{equation}
%which is convex since it is a positive linear combination of the convex functions $e^{-s \nu_{ji}}$ plus the constant term $-\az$. We also notice that $D_i(0)=0$ and $D_i(+\infty) = +\infty$ when \eqref{eq:as1} holds. 
%
%On the other hand, the numerator of \eqref{eq:taui}, 
%\begin{equation*}
%  R_i(s) := \ldi +
%  s \BXi,  
%\end{equation*}
%is a straight line crossing the vertical axis at $\ldi <0$, and we can
%assume that its slope, $\BXi$, is positive. Otherwise, the $\bar
%X(t)$ process is at the boundary of $\latt$, and therefore no reaction is
%pointing outside the lattice, $\latt$. We therefore set
%$\tau_i^* = +\infty$.
%
%Let us define $s_{0,i}$ as the root of the numerator $R_i(s)$, \ie,
%
%\begin{equation}
%\label{eq:root}
%s_{i} := -\ldi/\BXi \PERIOD
%\end{equation}
%By direct substitution of \eqref{eq:root} into \eqref{eq:DEN}, we
%obtain
%\begin{equation}
%\label{eq:DD}
%D_i(s_{i}) = -\az + \sj \aj\dnux\COMMA
%\end{equation}
%and
%\begin{equation}
%\label{eq:DDP}
%D'_i(s_{i}) =  -\sj \aj \nu_{ji} \dnux\PERIOD
%\end{equation}
%In order to determine whether $\tau_i^*<\infty$ or $\tau_i^*=\infty$,
%we have to analyse all possible cases regarding the pair
%$(R_i(s),D_i(s))$. 
%
%Indeed, note that 
%\begin{equation*}
%  D'_i(0) = - \sj \aj \nu_{ji},
%\end{equation*}
%and if $D'_i(0) \geq 0$, which could be interpreted as a drift
%pointing to the boundary, then $D_i(s)$ is monotonically
%increasing in $[0,+\infty)$. This situation is illustrated in Figure \ref{fig:RD}: in the left panel, we see the pair $(R_i(s),D_i(s))$;
%in the right panel, we see the quotient $\tau_i(s) = R_i(s)/D_i(s)$.  The
%function $\tau_i$ achieves its maximum, $\tau_i^*$, at a unique point,
%$\tsi$.
%\begin{figure}[h!]
%\begin{minipage}[b]{0.48\linewidth}
%
%%\centering
%\includegraphics[scale=0.6]{Figures/RD1.pdf}
%%\vspace{-12pt}
%%\label{fig:cherleaps}
%\end{minipage}
%%\hspace{0.02\linewidth}
%\begin{minipage}[b]{0.48\linewidth}
%%\centering
%\includegraphics[scale=0.6]{Figures/tauR.pdf}
%%\vspace{-15pt}
%\end{minipage}
%\caption{Left: Numerator $R_i(s)$ and denominator
%  $D_i(s)$. Right: Quotient $\tau_i(s) =
%  R_i(s)/D_i(s)$. Both plots are for the case $D'_i(0) \geq 0$.}
%\label{fig:RD}
%\end{figure}
%If $D_i'(0)<0$, which can be interpreted as a drift pointing to $+\infty$, the value of $\tau_i^*$
%depends on $\BXi$, \ie, on the size of the slope of $R_i(s)$. Observe
%that $D_i(s)$ is then negative in an interval $(0,d_i)$, with
%$D_i(d_i)=0$, and in general there is not a closed form for $d_i$.
%Also, since $D_i(s)$ and $R_i(s)$ may have opposite signs for some
%$s\leq \max(s_{i},d_i)$, this allows for artificially negative
%values of $\tau_i$, which should not be taken into account.
%
%The value of $\tau_i^*$ is finite or $+\infty$ according to the sign of  $D_i(s_{i})$. These three cases are shown in the left panel of the Figure \ref{fig:tauRR2}. 
%When $\BXi$ is large enough, \ie, when
%$D_i(s_{i})<0$, we can see in the right panel of Figure
%\ref{fig:tauRR2} that $\tau_i* = +\infty$. This is true because the
%limit of $\tau_i(s)$, as $s\rightarrow d_i^+$, is
%$+\infty$. Therefore, if $\BXi$ is far from the boundary and the drift
%is pointing to $+\infty$, we can take $\tau_i$ to be as large as we wish.
%-----
%\begin{figure}[h!]
%\begin{minipage}[b]{0.48\linewidth}
%
%%\centering
%\includegraphics[scale=0.6]{Figures/RD2.pdf}
%%\vspace{-12pt}
%%\label{fig:cherleaps}
%\end{minipage}
%%\hspace{0.02\linewidth}
%\begin{minipage}[b]{0.48\linewidth}
%%\centering
%\includegraphics[scale=0.6]{Figures/tauR1.pdf}
%%\vspace{-15pt}
%\end{minipage}
%\caption{{In this case $\sj \aj \nu_{ji} > 0$. Left: relative positions of $(R_i(s),D_i(s))$, depending on the sign of $D_i(s_{i})$.  Right: $\tau_i(s) = R_i(s)/D_i(s)$ in the case $D_i(s_{i})<0$.}}
%\label{fig:tauRR2}
%\end{figure}
%The two other cases are as follows: if $D_i(s_{i})>0$, it means that
%$\BXi$ is, in a certain sense, close to the boundary, and even if the
%drift is pointing to $+\infty$, there exists an upper bound for
%$\tau_i$. This is illustrated in the left part of Figure
%\ref{fig:tauR23}, where $\tau_i^*$ is the maximum to the right of
%$s_{i}$.  Finally, if $D_i(s_{i})=0$, then $\tau_i^*$ can be
%obtained as the limit of $\tau_i(s)$ as $s\rightarrow d_i^+$. By
%l'H\^{o}pital's rule, we have that $\tau_i^* = \BXi/D'_i(s_{i})$.
%----
%\begin{figure}[h!]
%\begin{minipage}[b]{0.48\linewidth}
%
%%\centering
%\includegraphics[scale=0.60]{Figures/tauR2.pdf}
%%\vspace{-12pt}
%%\label{fig:cherleaps}
%\end{minipage}
%%\hspace{0.02\linewidth}
%\begin{minipage}[b]{0.48\linewidth}
%%\centering
%\includegraphics[scale=0.65]{Figures/tauR3.pdf}
%%\vspace{-15pt}
%\end{minipage}
%\caption{{The other two cases for $\tau_i(s)$ when $\sj \aj \nu_{ji} > 0$. Left: $D_i(s_{i})>0$. Right: $D_i(s_{i}) = 0$.}}
%\label{fig:tauR23}
%\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% Summarize
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%We can summarize the previous discussion as follows: 
%We then have the following possibilities: If $\nu_{ji}\geq
%0$, for all $j$, then naturally $\tau_i^*{=}+\infty$; otherwise, we have the
%following three cases:
%\begin{enumerate}%[(I)]%[(a)]
%\label{summarizeChernoff}
%\item \label{eq:tsi} $D_i(s_{i}){>}0$.  In this case,
%  $\tau_i(s_{i}){=}0$ and $D_i(s)$ is positive and increasing as 
%  $\forall s \geq s_{i}$. Therefore, $\taui$ is equal to the ratio of two
%  positive increasing functions. The numerator, $R_i(s)$, is a linear
%  function and the denominator, $D_i(s)$, grows exponentially
%  fast. Then, there exist an upper bound, $\tau_i^*$, and a unique
%  number, $\tsi$, which satisfy $\tau_i(\tsi){=}\tau_i^*$. We developed an algorithm in  \cite{ourSL} for approximating $\tsi$, using the relation
%  $\tau'_i(\tsi){=}0$.
%\item If $D_i(s_{i}){<} 0$, then  $\tau_i^*{=}+\infty$.
%\item If $D_i(s_{i}){=}0$, then $\tau_i^*{=}\BXi/D'_i(s_{i})$.
%\end{enumerate}
%Here $s_{i} {:=}-\ldi/\BXi$.
%\subsection{Hybrid Path Generation}
%\label{sec:hybrid}
Using the exact pre-leap method we developed in \cite{ourSL,ourML} for single-level and multilevel hybrid schemes, allows us to switch adaptively between the tau-leap and an exact method.
%In this section, we briefly summarize our previous work on hybrid path generation, presented in \cite{ourSL} and \cite{ourML}. The main idea behind the hybrid algorithm is the following. 
%A path generated by an exact method (like SSA or MNRM) never exits the lattice, $\latt$, although the computational cost may be unaffordable due to many small inter-arrival times typically occurring when the process is ``far'' from the boundary. A tau-leap path, which may be cheaper than an exact one, could leave the lattice at any step. The probability of this event depends on the size of the next time step and the current state of the approximate process, $\BX$. This one-step exit probability could be large, especially when the approximate process is ``close'' to the boundary. We developed in \cite{ourSL} a Chernoff-type of bound to control the mentioned one-step exit probability. Even more, 
By construction, the probability that one hybrid path exits the lattice, $\latt$, can be estimated by
\begin{align*}%\label{eq:pac}
\prob{ A^c} \leq \expt{1-(1-\delta)^{N_{\TL}}} = \delta \expt{N_{\TL}} - \frac{\delta^2}{2}(\expt{N_{\TL}^2} - \expt{N_{\TL}}) + o(\delta^2),
\end{align*}
where  $\bar{\omega} \in A$ if and only if the whole hybrid path, $\seqof{\bar X(t_k,\bar \omega)}{k=0}{K(\bar{\omega})}$, belongs to the lattice, $\latt$, $\delta {>}0$ is the one-step exit probability bound, and $N_{\TL}(\bar \omega) {\equiv} N_{\TL}$ is the number of tau-leap steps in a hybrid path. Here, $A^c$ is the complement of the set $A$.

To simulate a hybrid path, given the current state of the approximate process, $\BX$, we adaptively determine whether to use an exact or the tau-leap method  for the next step. This decision is based on the relative computational cost of taking an exact step versus the cost of taking a Chernoff tau-leap step.
%We show the switching rule in Algorithm \ref{alg:sel}.
Instead, in the present work, at each time step, we adaptively determine which reactions are suitable for using the exact method and which reactions are suitable for the Chernoff tau-leap method.

%\subsection{Outline of this Work}
%In Section \ref{sec:mixedpaths}, we present an algorithm to generate mixed trajectories; that is, the algorithm generates a trajectory using an exact method for the low activity channels and the Chernoff tau-leap method for the high activity ones. Then, inspired by the ideas of Anderson and Higham \cite{Anderson2012}, we propose an algorithm for coupling two mixed Chernoff tau-leap paths. This algorithm uses four building blocks that result from the combination of the MNRM and the tau-leap methods.
%In Section \ref{sec:mlmcerr}, we propose a mixed MLMC estimator. Next, we introduce a global error decomposition and show that the computational complexity of our method is of order $\Ordo{TOL^{-2}}$.
%Finally, we show the automatic procedure that estimates our quantity of interest within a given prescribed relative tolerance, up to a given confidence level.
%Next, in Section \ref{sec:cvar}, we present a novel control variate technique to reduce the variance of the quantity of interest at level 0. 
%%and finally, we indicate a Bayesian strategy for estimating parameters involved in the computational work optimization problem. 
%In Section \ref{sec:examples}, the numerical examples illustrate the advantages of the mixed MLMC method over the hybrid MLMC method presented in \cite{ourML} and to the SSA. Last, Section \ref{sec:conclusions} presents our conclusions and suggestions for future work.

\section{Generating Mixed Paths}
\label{sec:mixedpaths}
In this section we explain how mixed paths are generated. First, we present the splitting heuristic; that is, we discuss how to partition the set of reaction channels at each decision time. Then, we present the one-step mixing rule, which is the main building block for constructing a mixed path. Finally, we show how to couple two mixed paths.

\subsection{The Splitting Heuristic}
\label{sec:splitting}
In this section, we explain how we partition the set of reaction channels, $\mathcal{R}{:=}\{1,...,J\}$, into $\mathcal{R}_{\TL}$ and $\mathcal{R}_{\MNRM}$.

Let $(t,x)$ be the current time and state of the approximate process, $\bar X$, and $H$  be the next decision (or synchronization) time (given by the Chernoff tau-leap step size $\tau_{Ch}=\tau_{Ch}(x,\delta)$ and the time mesh).
%Let $\mathcal{R}$ be the set of reaction channels, in this case $|\mathcal{R}|{=}J$. 
We want to split $\mathcal{R}$ into two subsets, $\mathcal{R}_{\MNRM}$ and $\mathcal{R}_{\TL}$, such that 
the expected computational work of reaching $H$, starting at $t$, is minimal for all possible splittings.

The idea goes as follows. First, we define a  linear order on $\mathcal{R}$, based on the basic principle that we want to use tau-leap for the $j$-th reaction if its activity is high. This linear order determines $J{+}1$ possible splittings, out of $2^J$. 
In order to measure the activity, it turns out that using only the propensity functions evaluated at $x$, that is, $a_j(x)$, is not enough. This is because the $j$-th reaction could affect components of $x$ with small values. If this is the case, this determines small Chernoff tau-leap step sizes. In order to avoid this scenario, we penalize the $j$-th reaction channel if it has a high exit probability. We approximate this exit probability using a Poisson distribution for each dimension of $x$. For example, let $\nu_{ji}$ be the $i$-th component of the $j$-th reaction channel. If $\nu_{ji}<0$, then the probability that a Poisson distributed random variable with rate $a_j(x)(H{-}t)$ is greater than $x_i/\nu_{ji}$ measures how likely species $x_i$ can become negative in the interval $H{-}t$, independently of reactions $j'{\in} \mathcal{R}$, $j{\neq} j'$.
Let $I_j{:=} \{i:\nu_{ji} <0\}$,
%$$\theta_{ji}{:=} \left\{ 
%\begin{array}{ll}
%\prob{x_i + \nu_{ji}\, \mathcal{P}(a_j(x)(H{-}t)) <0\SEP x} & \text{ for } \nu_{ji} <0, \\
%0 &\text{otherwise}
%\end{array}
%\right.\PERIOD
%$$
\begin{align}\label{eq:thetas}
\theta_{j}:= \left\{ 
\begin{array}{ll}
\prob{\mathcal{P}(a_j(x)(H{-}t)) > \min_{i\in I_j} \{-\frac{x_i}{\nu_{ji}}\}\SEP x} & \text{if } I_j \neq \emptyset\\
0 &\text{otherwise}
\end{array}
\right. .
\end{align}
Then, the penalty weight for $a_j(x)$ is $1{-}\theta_j$. We define $\tilde{a}_j(x){:=}(1{-}\theta_j)a_j(x)$. The linear order is then a permutation, $\sigma$, over $\mathcal{R}$ such that
$$\tilde{a}_{\sigma(j)}(x)>\tilde{a}_{\sigma(j+1)}(x), \,\,\,  j{=}1,...,J{-}1 \PERIOD$$
 
Second, we find among the $J{+}1$ partitions the one with optimal work. This  is the computational work incurred when performing one step of the algorithm using tau-leap for the reactions $\mathcal{R}_{TL}$ and the MNRM for the reactions $\mathcal{R}_\MNRM$. The work corresponding to $\mathcal{R}_{\TL}$ is 
\begin{align}\label{eq:worktl}
\text{Work}(\mathcal{R}_{\TL},x,t):=\frac{H{-}t}{\min\{\tau_{Ch},H{-}t\}} \left(C_{s} + \sum_{j\in \mathcal{R}_{\TL}}C_P(a_{j}(x) \tau_{Ch})\right) \COMMA
\end{align}
where $C_s$ is the work of computing the split (see Section \ref{sec:cs}), and $C_P(\lambda)$ is the work of a Poisson random variate with rate $\lambda$. 
The factor $\frac{H{-}t}{\min\{\tau_{Ch},H{-}t\}}$ takes into account the number of steps required to reach $H=H(t)$ from $t$.
For the Gamma simulation method developed by Ahrens and Dieter in \cite{poivariate}, which is the one used by MATLAB, $C_P$ is defined as
\begin{align*}
C_P(\lambda):=
\left\{ \!
\begin{array}{ll}
b_1{+}b_2 \ln \lambda & \text{ for } \lambda >15 \\
b_3{+}b_4 \lambda & \text{ for } \lambda \leq 15
\end{array}
\right. \PERIOD
\end{align*}
In practice, it is possible to estimate $b_i$, $i{=}1,2,3,4$ using Monte Carlo sampling and a least squares fit. 
%, as shown in Figure \ref{fig:ahrens}. 
For more details, we refer to \cite{ourSL}.
%\begin{figure}[h!]
%\centering
%\begin{minipage}{0.49\textwidth}
%\includegraphics[width=\textwidth]{Figures/ahrens.pdf}\end{minipage}
%\hfill
%\begin{minipage}{0.49\textwidth}
%\includegraphics[width=\textwidth]{Figures/ahrens_lin.pdf}
%\end{minipage}
%%\caption{Left: The computational work (runtime) model for generating a Poisson random variate, using the Gamma method by Ahrens and Dieter \cite{poivariate}. Right: Linear growth detail, for $\lambda \in[0,15]$.}
%\label{fig:ahrens}
%\end{figure}

Similarly, the work corresponding to $\mathcal{R}_\MNRM$ is
$$\text{Work}(\mathcal{R}_{\MNRM},x,t):=\frac{H{-}t}{\min\{\tau_{\MNRM},H{-}t\}} \, C_{\MNRM} \COMMA$$
where the constant $C_{\MNRM}$ is the work of an MNRM step and $\tau_{\MNRM}{=}\LP \sum_{j\in \mathcal{R}_\MNRM} a_j(x)\RP^{-1}$.
%To estimate the computational cost (cost or work from now on), we use a Poisson cost model, presented in \cite{ourSL}. The cost corresponding to the tau-leap reactions is $$(H{-}t)/\tau_{Ch} \sum_{j\in \mathcal{R}_{TL}} C_P(a_{j}(x) \tau_{Ch})$$
%where $C_P(\lambda)$ is the cost of a Poisson random variate with rate $\lambda$.
%$$(H{-}t)/\min\{\tau_{Ch},H{-}t\} \sum_{j\in \mathcal{R}_{TL}} ahrens(a_{j}(x) \tau_{Ch})$$
%A similar calculation is performed for the cost of the MNRM part.
\subsection{On the Work required to the Splitting Heuristic, $C_s{=}C_s(J)$}
\label{sec:cs}
The work required to perform the splitting includes the work required to determine $\text{Work}(\mathcal{R}_{\TL})$ and $\text{Work}(\mathcal{R}_{\MNRM})$, both defined in Section \ref{sec:splitting}. The linear order previously defined determines $J{+}1$ possible splittings, $\mathcal{S}_i$, $i{=}0,...,J$, as follows:
\begin{center}
\begin{tabular}{lc|c}
& $\mathcal{R}_{\TL}$ & $\mathcal{R}_\MNRM$ \\
\hline 
$\mathcal{S}_0$ & $\emptyset$ & $\mathcal{R}$ \\
$\mathcal{S}_1$ & $\{\sigma^{-1}(1)\}$ & $\{\sigma^{-1}(2),..,\sigma^{-1}(J)\}$ \\
$\mathcal{S}_2$ & $\{\sigma^{-1}(1),\sigma^{-1}(2)\}$ & $\{\sigma^{-1}(3),..,\sigma^{-1}(J)\}$ \\
$\vdots$ & & \\
$\mathcal{S}_{J}$ & $\mathcal{R}$ & $\emptyset$
\end{tabular} .
\end{center}
The cost of computing each of the $J{+}1$ splits is dominated by the cost of determining the Chernoff tau-leap step size, $\tau_{Ch}$ (see \eqref{eq:worktl}). As we observe in \cite{ourSL}, the work of computing a single $\tau_{Ch}$ is linear on $J$.
Then,  in order to avoid  $J^2$ complexity of the splitting rule, we implement a local search instead of computing $J$ $\tau_{Ch}$'s, to keep the complexity of $C_s$ linear on $J$.
The main idea is to keep track of the last split at each decision time, assuming that the propensities do not vary widely between.
If that is the case, we can just evaluate the previous split, $\mathcal{S}_\kappa$, and its neighbors, $\kappa{-}1$ and $\kappa{+}1$. Then, the cost of the splitting rule is on the order of three computations of a Chernoff step size. It turns out that this local search is very accurate for the examples we worked on. In order to avoid being trapped in local minima, a randomization rule may be applied.%, but we pospone that for future work.

\begin{rem}[Pareto Splitting rule]
\label{rem:pare}
Instead of computing a cost-based splitting at each decision time, the following rule can be applied:
$$\mathcal{R}_{\TL} \text{ is defined s.t. } \frac{\sum_{j\in \mathcal{R}_{\TL}} \tilde{a}_{\sigma(j)}}{\sum_{k=1}^J \tilde{a}_k} \geq \nu \COMMA$$
where $\nu$ is a problem-dependent threshold, which can be estimated using the cost-based splitting rule.
The idea is to use the tau-leap method for a $(100\times \nu) \%$  of the penalized activity (measured as before using the $\tilde{a}_j$'s), and an exact method for the other channels.

This rule is adaptive because it depends on the current state of the process, but it does not take into account the computational cost of the resulting partition of $\mathcal{R}$.
The advantage of this rule is that it is three times faster than the previous one. For the examples we worked on, the overall average gain in terms of computational work in a whole mixed path is about 45\% of the total work.
%\begin{itemize}
%\item $\text{Work}(\text{splitting rule}) \geq 3 C_{\TL}$, for some problems this can be too expensive.
%\item For example, if the optimal split does not change ``too much'' along a path, we can get $\text{Work}(\text{splitting rule}) \geq C_{\TL}$, using a fixed splitting rule based on $\tilde{a}_{\sigma(j)}$ only.
%\end{itemize}
%\textbf{Idea}: {\sl ``Use tau-leap for a $100\nu$ percent of the activity, and an exact method for  $100(1{-}\nu)$.''}

\end{rem}
\subsection{The one-step Mixing Rule}
In this section we present the main building block for simulating a mixed path. Let $x {=} \BX$ be the current state of the approximate process, $\bar X$. Therefore, the expected time step  of the MNRM is  given by $1/{a_0(x)}$. %Let $\tau_{Ch}{=}\tau_{Ch}(x,\delta)$  be the Chernoff tau-leap step.
To move one step forward using the MNRM, we should compute at least $a_0(x)$ and sample a uniform random variable.
On the other hand, to move one step forward using the mixed Chernoff tau-leap method, we need first to compute the split, then compute the tau-leap increments for the reactions in the tau-leap set, $\mathcal{R}_\TL$, and finally compute the MNRM steps for the reactions in the set $\mathcal{R}_\MNRM$, as discussed in Section \ref{sec:cs}.

To avoid the overhead caused by unnecessary computation of the split, we first estimate the computational work of moving forward from the current time, $t$, to the next grid point, $\tilde{T}$, by using the MNRM only. If this work is less than the work of computing the split, we take an exact step. 
%This motivates us to define  $K_1$ as the ratio between the work of computing the split and the work of perfoming and exact step using the MNRM.
%Otherwise, we compute the split and decide 
%whether to take an SSA step or a tau-leap one, according to the comparison between $\tau_{Ch}$ and  $K_2/{a_0(x)}$. Here $K_2{=}K_2(x,\delta)$ is defined as the work of taking a Chernoff tau-leap step given the current state of the process, divided by the work of taking an SSA step plus the work of computing $\tau_{Ch}$. 
%As we mentioned, associated with each type of step, there is computational work. In the first case, when $K_1/a_0(x)>T_0{-}t$, the work is $C_1$, and includes the computation of $1/a_0(x)$ and the generation of two uniform random variates. In the same way, when $K_1/a_0(x) >T_0{-}t$ and $K_2/a_0(x)> \tau_{Ch}$, the work is $C_2$, and  involves the work contained in $C_1$ and of computing $\tau_{Ch}(x,\delta)$, which is denoted by $C_3$.
%On the other hand, when a Chernoff tau-leap step is taken, we have not only the constant work, $C_3$, but also \emph{variable} work, which is the work of generating the Poisson random variates. %The work of a Chernoff tau-leap step is composed by the work of computing $\tau_{Ch}(x,\delta)$ which is constant, and one variable work, which is 
%The latter work is a function of the propensities of all the reaction channels, namely, $a(x)\tau_{Ch}(x,\delta)$. We model the computational work of generating \emph{one} Poisson random variate according to \cite{poivariate}, and this work is denoted by $C_P(\cdot)$. 
%%We now briefly describe $C_P(\cdot)$. 
%In the Gamma simulation method developed by Ahrens and Dieter in \cite{poivariate}, which is used by MATLAB, the work grows like $b_1{+}b_2 \ln \lambda$ where $\lambda >15$ is the rate of the Poisson random variable. For $\lambda \leq 15$, the growth is linear. In practice, it is possible to estimate $b_1$ and $b_2$ using a Monte Carlo method with a least squares fit, as shown in Figure \ref{fig:ahrens}. 
%
%Summarizing,
%$K_1:= \frac{C_3}{C_1}$, and $K_2(\BX,\delta) := \frac{C_3 + \sum_{j=1}^JC_P(a_j(\BX)\tau_{Ch}(\BX,\delta))}{C_1+C_3}.$
%
%Observe that $K_2(x,\delta) \rightarrow \frac{C_3+J b_1}{C_1+C_3}=: \tilde{C}>0$ as $\delta \rightarrow 0$.
\begin{algorithm}[h!]
\caption{The one-step mixing rule. Inputs: the current state  of the approximate process, $\BX$, the current time, $t$, the values of the propensity functions evaluated at $\BX$, $\seqof{a_j(\BX)}{j=1}{J}$, the one-step exit probability bound $\delta$, the next grid point, $\tilde{T}$, and the previous optimal split, $\kappa$. Outputs: the tau-leap set, $\mathcal{R}_\TL$, the exact set, $\mathcal{R}_\MNRM$, and the new optimal split $\kappa$. %Notes: 
%based on $\expt{\tau_{\ssa}(\BX)}=1/a_0(\BX)$ and $\tau_{Ch}(\BX,\delta)$, this algorithm adaptively selects which  method to use: MNRM or Chernoff tau-leap (TL).  %Here $K_1$ and $K_2$ are positive constants. 
%We denote by $\tau_{\MNRM}$  ($\tau_{Ch}$) the step size when the decision is to use the $\MNRM$  (tau-leap) method. 
%For example, select $f_{\ssa}$ if $C \tau_{\ssa}<\tau_{Ch}$, where $C$ is a fixed constant.
}
\label{alg:sel}
\begin{algorithmic}[1]
	\REQUIRE $a_0 \leftarrow \sum_{j=1}^J a_j > 0$
	\IF{$K_1/a_0 < \tilde{T}-t$} 
		\STATE Compute $\theta_{j}$, {\small $j{=}1,..,J$} (see \eqref{eq:thetas})
		\STATE $\tilde{a}_{\sigma(j)} \leftarrow \text{Sort} \{(1{-}\theta_{j})a_j\}$ descending, {\small$j{=}1,..,J$}
		\STATE $\mathcal{S}_i \leftarrow $ Compute the splits taking into account the previous optimal split $\kappa$
		\STATE $(\mathcal{R}_\TL, \mathcal{R}_\MNRM, \kappa) \leftarrow $ Take the minimum work split
		%\STATE $\tau_{Ch} \leftarrow $ Compute Chernoff step size of $\mathcal{R}_{\TL}$
		\RETURN $(\mathcal{R}_\TL, \mathcal{R}_\MNRM,\kappa)$
	\ELSE
		\RETURN $(\emptyset, \mathcal{R},\kappa)$
	\ENDIF
\end{algorithmic}
\end{algorithm}
In order to compare the mentioned computational costs, we define $K_1$ as the ratio between the cost of computing the split, $C_s$, and the cost of computing one step using the MNRM.
%For further details on the switching rule, we refer to \cite{ourSL}.
\begin{rem}[Comparison with the one-step hybrid rule]
In \cite{ourSL} we developed a hybrid method, which, at each decision point, determines which method, exact or tau-leap, is cheaper to apply to the whole set of reactions. That is, in the hybrid method, we have either $\mathcal{R}_\TL =\emptyset$ and $\mathcal{R}_\MNRM =\mathcal{R}$ or $\mathcal{R}_\TL =\mathcal{R}$ and $\mathcal{R}_\MNRM =\emptyset$. Then, the mixed method can be seen as a generalization of the hybrid one. The key difference is in the cost of the decision rule, which, as we saw in Section \ref{sec:cs}, in the mixed method is on the order of three times the computation of the Chernoff step size. This difference can be significant in some problems. A Pareto splitting rule may be able to recover the cost of the hybrid one-step decision rule.
\end{rem}

\subsection{The Mixed-Path Algorithm}
In this section, we present a novel algorithm (Algorithm \ref{alg:mixed}) that combines the approximate Chernoff tau-leap method and the exact MNRM to generate a whole hybrid path. This algorithm combines the advantages of an exact method (expensive but exact) and the tau-leap method (may be cheaper but has a discretization error and a positive probability of exiting the lattice). 
This algorithm automatically and adaptively partitions the reactions into two  subsets, $\mathcal{R}_\TL$ and $\mathcal{R}_\MNRM$, using a computational work criterion.
%On the other hand, a tau-leap path, which may be cheaper than an exact one, could leave the lattice at any step. 
%Regarding the exit error, it depends on the size of the next time step and the current state of the approximate process, $\BX$. This one-step exit probability could be large, but, %especially when the approximate process is ``close'' to the boundary. 
%as we mention in Section \ref{sec:chernoff}, we can control this one-step exit probability adaptively, by adjusting the tau-leap step size. 
Since a mixed path consists of a certain number of exact/approximate steps, it may also exit the lattice, except in those steps in which the tau-leap method is not applied; that is, when $\mathcal{R}_\TL$ is empty. 
The idea of this algorithm is to apply, at each decision point, the one-step mixing rule (Algorithm \ref{alg:sel}) to determine the sets $\mathcal{R}_\TL$ and $\mathcal{R}_\MNRM$, and then to apply the corresponding method.

\begin{algorithm}[h!]
\caption{The mixed-path algorithm. 
Inputs: the initial state, $X(0)$, the propensity functions, $\seqof{a_j}{j=1}{J}$, the stoichiometric vectors, $\nu{=}\seqof{\nu_j}{j=1}{J}$, the final time, $T$, and the one-step exit probability bound, $\delta$. 
Outputs: a sequence of states, $\seqof{\bar{X}(t_k)}{k=0}{K}$, % such that $t_K \leq T$. If $t_K<T$, then the path exited the $\latt$ lattice before the final time $T$. 
and the number of times, $\NTL$, that the tau-leap method was successfully applied
(\ie, $\bar{X}(t_k)\in \latt$, we applied the tau-leap method and we obtained an $\bar{X}(t_{k+1})\in \latt$). %, the number of SSA steps such that $K_1/a_0(\bar X(t))>t_k{-}t$ is true, $\NSSAKone$, and the number of MNRM steps such that $K_1/a_0(\bar X(t)) <t_k{-}t$, $\NSSAKtwo$. 
Notes: given the current state, $next_{\MNRM}$ computes the next state using the MNRM method. Here, $t_i$ denotes the current time at the $i$-th step, and $\tau_{Ch}(\mathcal{R}_\TL)$ is the Chernoff step size associated with $\mathcal{R}_\TL$.
}
%\caption{The hybrid tau-leap algorithm. Inputs: the initial state, $X(0)$, the propensity functions, $\seqof{a_j}{j=1}{J}$, the stoichiometric vectors, $\nu{=}\seqof{\nu_j}{j=1}{J}$, and %a time mesh,  $\seqof{t_k}{k=0}{K}$, such that $t_K$ is equal to $T$, 
%the final time, $T$. Outputs: a sequence of states, $\seqof{\bar{X}(t_k)}{k=0}{K}$, %the last system state, $\X{K}$ (and its respective time  $t_K$), 
%such that $\X{K} \in \latt$  and $t_K \leq T$. It also returns the number of times, $N_{TL}$, the tau-leap method was successfully applied (that is, when $\X{K} \in \latt$), the number of SSA steps such that $K_1/a_0(\bar X(t))>t_k{-}t$ is true, $\NSSAKone$, and the number of SSA steps such that $K_1/a_0(\bar X(t)) >t_k{-}t$ is false and $K_2(\bar X(t))/a_0(\bar X(t))> \tau_{Ch}$ is true, $\NSSAKtwo$ (see Algorithm \ref{alg:sel}). Notes: given the current state, $next_{\ssa}$ computes the next state using the SSA method. Here, $t_i$ denotes the current time at the $i$-th step. 
%}
\label{alg:mixed}
\begin{algorithmic}[1]
\STATE $i \leftarrow 0$, $t_i \leftarrow t_0, \bar{X}(t_i) \leftarrow X(0), \bar Z \leftarrow X(0)$
\STATE $\mathcal{S}_j \leftarrow$ Compute splits, {\small $j{=}0,...,J$} 
\STATE $\kappa \leftarrow \arg \min_j \text{Work}(\mathcal{S}_j)$
\WHILE {$t_i<T$}
\STATE $\tilde T \leftarrow $ next grid point greater than $t_i$
%\STATE $\theta_{j} \leftarrow \max_{i, \nu_{ij}{<}0}\{\text{Klar}(a_j(x)\tau_{SSA},\lceil -x_i/\nu_{ij}\rceil)\} $,  $j{=}1,..,J$
%	\STATE $\tilde{a}_{\rho(j)} \leftarrow \text{Sort} \{(1{-}\theta_{j})a_j(x)\}$ desc., s.t. $\tilde{a}_{\rho(j)}{\geq} \tilde{a}_{\rho(j)+1}$,  $j{=}1,..,J$
\STATE $(\mathcal{R}_{\TL},\mathcal{R}_{\MNRM},\kappa) \leftarrow $ Algorithm \ref{alg:sel} with $(\bar Z,t_i,\seqof{a_j(\bar Z)}{j=1}{J},\delta, \tilde T,\kappa)$
%\STATE $\tilde{T} \leftarrow $ next grid point
%\STATE $\tau_{Ch} \leftarrow \min \{\text{compute tau Chernoff for }\mathcal{R}_{\TL},\tilde{T}{-}t\}$
\IF {$\mathcal{R}_{\TL}\neq \emptyset$}
%\STATE $\tau_{Ch} \leftarrow $ Compute Chernoff step size using $j {\in} \mathcal{R}_{\TL}$
\STATE $\Delta_{\TL} \leftarrow \mathcal{P}(a_j(\bar Z) \tau_{Ch}(\mathcal{R}_\TL)) \nu_{j}$,   for {\small$j{\in}\mathcal{R}_{\TL}$}
\STATE $H \leftarrow t_i+ \tau_{Ch}(\mathcal{R}_\TL)$
\ELSE
\STATE $H \leftarrow \min \{t_i{-}\log(r)/\sum_{j}a_j,T\}$,   $r{\sim} \text{Unif}(0,1)$
\ENDIF
\IF {$\mathcal{R}_{\MNRM}\neq \emptyset$}
\WHILE {$t_i<H$}
\STATE $(\bar{Z},t_i) \leftarrow next_{\MNRM}(\bar{Z},\mathcal{R}_e,t_i,H)$
%\STATE $t\leftarrow next_time(t,H)$
\ENDWHILE
\ENDIF
\STATE $\bar{Z} \leftarrow \bar{Z}+\Delta_{\TL}$
\IF {$\bar Z \in \latt$}
	\STATE $\NTL \leftarrow \NTL+1$
	\STATE $t_{i+1} \leftarrow H$
\ELSE
	\RETURN $(\seqof{\bar{X}(t_k)}{k=0}{i},\NTL)$
\ENDIF

\STATE $i \leftarrow i+1$
\STATE $\bar X(t_i) \leftarrow \bar Z$
\ENDWHILE
\RETURN $(\seqof{\bar{X}(t_k)}{k=0}{i},\NTL)$
\end{algorithmic}
\end{algorithm}

\subsection{Coupled Mixed Paths}\label{cmp}
In this section, we explain how to couple two mixed paths. This is essential for the multilevel estimator. The four algorithms that are the building blocks of the coupling algorithm were already presented in \cite{ourML}. The novelty here comes from the fact that the coupled mixed algorithm may have to run the four algorithms concurrently in the sense of the time of the process, $t$. In this section, we denote with a bar $\bar \cdot$ and a double bar $\dbar \cdot$ coarse and fine grid-related quantities.

We now briefly describe the mixed Chernoff coupling algorithm, \ie,  Algorithm \ref{alg:coupled}. Let $\bar X$ and $\bar{\bar X}$ be two mixed paths, corresponding to two nested time discretizations, called coarse and fine, respectively.
Assume that the current time is $t$, and we know the states, $\bar X(t)$ and $\bar{\bar X}(t)$, the next grid points at each level, $\bar t$, $\dbar t$, and the corresponding one-step exit probabilities, $\bar \delta$ and $\dbar \delta$.
Based on this knowledge, we have to determine the four sets $(\bar{\mathcal{R}}_{\TL},\bar{\mathcal{R}}_{\MNRM},\dbar{\mathcal{R}}_{\TL},\dbar{\mathcal{R}}_{\MNRM})$, that correspond to four algorithms, B1, B2, B3 and B4, that we use as building blocks. Table \ref{tab:algs} summarizes them.
\begin{table}[h]
\centering
\begin{tabular}{l|c|c}
  & $\bar{\mathcal{R}}_{\TL}$ & $\bar{\mathcal{R}}_{\MNRM}$ \\ \noalign{\smallskip} \hline\noalign{\smallskip} 
$\dbar{\mathcal{R}}_{\TL}$ &B1 & B2 \\
$\dbar{\mathcal{R}}_{\MNRM}$ &B3 & B4
\end{tabular}
%\label{tab:algs}
\medskip
\caption{Building blocks for simulating two coupled mixed Chernoff tau-leap paths.  
Algorithms B1 and B2 are presented as Algorithms 2 and 3 in \cite{Anderson2012}. Algorithms B3 and B4 can be directly obtained from Algorithm B2 (see \cite{ourML}).}\label{tab:algs}
\end{table}
In order to do that, the algorithm computes, independently, the sets $\mathcal{R}_\TL$ and $\mathcal{R}_\MNRM$ for each level, and the time until the next decision is taken, $H$, using Algorithm \ref{alg:timehor}. 
%At the beginning of each step of Algorithm \ref{alg:coupled} four set of reactions, which correspond to algorithms B1-B4, are determined. 
Next, it computes concurrently the increments due to each one of the sets (storing the results in $\Delta \bar X$ and $\Delta \dbar X$ for the coarse and fine grid, respectively). We note that the only case in which we use a Poisson random variates generator for the tau-leap method is in Algorithm B1 (Algorithm \ref{alg:bb1}). 
For Algorithms B2, B3 and B4, the Poisson random variables are simulated by adding independent exponential random variables with the same rate, $\lambda$, until exceeding a given time final time, $T$. The only difference in the latter blocks are the time points at which the propensities, $a_j$, are computed. For B2, the coarse propensities are frozen at time $t$, whereas for B3 the finer are frozen at $t$. In B4, the propensities are computed at each time step.
After arriving at time $H$, the four sets $(\bar{\mathcal{R}}_{\TL},\bar{\mathcal{R}}_{\MNRM},\dbar{\mathcal{R}}_{\TL},\dbar{\mathcal{R}}_{\MNRM})$ and the time until the next decision is taken, $H$, are determined again, and then all procedures are repeated until the simulation reaches the final time, $T$.

%We now explain algorithm B1 (TL-TL). Assume that tau-leap is chosen at the coarse and at the fine level, thus obtaining two time horizons, one for the coarse level, $\bar{H}$, and another for the fine level, $\dbar{H}$. In this case, the global time horizon will be $H{:=}\min\{\bar{H},\dbar{H}\}$. Since the chosen method in both grid levels is tau-leap, we need to freeze the propensities at the beginning of the corresponding intervals. In the coarse case, during the interval $[t,\bar{H})$ (the propensities  are equal to $a(\bar{X}(t)) {=:}\bar{a}$), and in the fine case during the interval $[t,\dbar{H})$ (the propensities  are equal to $a(\dbar{X}(t)) {=:}\dbar{a}$). 
%Suppose that $\bar{H}<\dbar{H}$. Then, we couple two Poisson random variables at time $t {=} \bar{H}$, using the idea described in Section \ref{sec:couplingpoiss}. 
%When time reaches $\bar{H}$, the decision of which method to use (and the corresponding step size) at the coarse level must be made again. Note that the propensities of the process at the fine grid will be kept frozen until $\dbar{H}$. The case when $\bar{H}>\dbar{H}$ is analogous to the one we described, but the decisions on the method and step size are made at the finer level, when time reaches $\dbar{H}$. It can also be possible that $\bar{H}=\dbar{H}$. In that case, the decision of which method to use (and the corresponding step size) must be made at the coarse and at the fine level.
%
%In the case of algorithm B2 (TL-MNRM), we assume that tau-leap is chosen at the coarse level, and MNRM at the fine level, obtaining two time horizons, one for the coarse level, $\bar{H}$, and another for the fine level, $\dbar{H}$. The only difference in how we determine the time horizons between algorithm B1 and B2 is that the time discretization grid points at the fine grid are not taken into account to determine $\dbar{H}$. Algorithm B2 is then applied until the simulation reaches $H{:=}\min\{\bar{H},\dbar{H}\}$. Suppose that $\dbar{H}<\bar{H}$. In this case, the process $\dbar{X}$ could take more than one step to reach $\dbar{H}$. At each step, the propensity functions $a(\dbar{X}(\cdot))$ are computed, but not the propensities for the coarse level, because in that case the tau-leap method is used. Note that the decision on which algorithm to use (B2 or another) is not made at those steps, but only when time reaches $\dbar{H}$. When time reaches $\dbar{H}$, the decision of which method to use (and the corresponding step size) at the fine level must be made again. In this case, the propensities at the coarse grid will be kept frozen until $\bar{H}$. The reasoning for the cases $\dbar{H}>\bar{H}$ and $\dbar{H}=\bar{H}$ are similar as before.
%
%The other two cases, that is, B3 and B4, are the same as B2. The only difference resides is when to update the propensity values, $\bar{a}$ and $\dbar{a}$. See Algorithm \ref{alg:coupled} for more details.
%It is clear in the preceding paragraphs, that the decision on which algorithm should be used for a certain time interval is made only at the horizon points.
%\begin{algorithm}[h!]
%\caption{Computes the Chernoff tau-leap step size. 
%}
%\label{alg:Ch}
%
%\end{algorithm}
\newcommand{\SWITCH}[1]{\STATE \textbf{switch} #1 \begin{ALC@g}}
\newcommand{\ENDSWITCH}{\end{ALC@g}\STATE \textbf{end switch} }
\newcommand{\CASE}[1]{\STATE \textbf{case} #1\textbf{:} \begin{ALC@g}}
\newcommand{\ENDCASE}{\end{ALC@g}}

\begin{algorithm}[h!]
\caption{\small Coupled mixed path. Inputs: the initial state, $X(0)$, the final time $T$, the propensity functions, $\seqof{a_j}{j=1}{J}$, the stoichiometric vectors, $\seqof{\nu_j}{j=1}{J}$, and  
two time meshes, one coarser $\seqof{t_i}{i=0}{N}$, such that $t_N {=} T$ and a finer one, 
$\seqof{s_j}{j=0}{N'}$, such that $s_0 {=} t_0$, $s_M {=} t_N$, and $\seqof{t_i}{i=0}{N} {\subset} \seqof{s_j}{j=0}{N'}$.
%Outputs: the last system state, in the coarse and in the fine time meshes, $\bar{X}$, $\dbar{X}$ (and its respective time $t$), such that $\bar{X}$, $\dbar{X} \in \latt$  and $t \leq T$. 
Outputs: a sequence of states evaluated at the coarse grid, $\seqof{\bar{X}(t_k)}{k=0}{K} \subset \latt$, %the last system state, $\X{K}$ (and its respective time  $t_K$), 
such that $t_K \leq T$, 
a sequence of states evaluated at the fine grid $\seqof{\dbar{X}(s_l)}{l=0}{K'} \subset \latt$,
such that $\bar{X}(t_K) \in \latt$ or $\dbar{X}(s_{K'}) \in \latt$.
%the last state in the coarse time mesh and the last state in the fine time mesh, $\bar{X}$, $\dbar{X}$ (and their respective times $\bar{t}$, and $\dbar{t}$). 
If $t_K<T$, both paths exit the $\latt$ lattice before the final time, $T$. 
It also returns the number of times the tau-leap method was successfully applied at the fine level and at the coarse level and the number of exact steps at the fine level and at the coarse level.
For the sake of simplicity, we omit sentences involving the recording of current state variables, counting of the number of steps, checking if the path jumps out of the lattice, the updating of the current split, $\kappa$, and the return sentence.}
\label{alg:coupled}
\begin{algorithmic}[1]
\begin{small}
\STATE $t \leftarrow t_0$; $\bar{X} \leftarrow X(0)$; $\dbar{X} \leftarrow X(0)$ %, out $\leftarrow$ \FALSE
\STATE $\bar{t} \leftarrow$ next grid point in $\seqof{t_i}{i=0}{N}$ larger than $t$ 
\STATE $(\bar{H},\bar{\mathcal{R}}_{\TL},\bar{\mathcal{R}}_{\MNRM},\bar{a})  \leftarrow$ Algorithm \ref{alg:timehor} with ($\bar{X}$,$t$,$\bar{t}$,$T$,$\bar{\delta}$)
\STATE $\dbar{t} \leftarrow$ next grid point in $\seqof{s_i}{i=0}{N}$ larger than $t$ 
\STATE $(\dbar{H},\dbar{\mathcal{R}}_{\TL},\dbar{\mathcal{R}}_{\MNRM},\dbar{a})  \leftarrow$ Algorithm \ref{alg:timehor} with ($\dbar{X}$,$t$,$\dbar{t}$,$T$,$\dbar{\delta}$)

	\WHILE{$t<T$}
		\STATE $H \leftarrow \min\{\bar H, \dbar H\}$
		\STATE $(B_1,B_2,B_3,B_4) \leftarrow$ split building blocks from $(\bar{\mathcal{R}}_{\TL},\bar{\mathcal{R}}_{\MNRM},\dbar{\mathcal{R}}_{\TL},\dbar{\mathcal{R}}_{\MNRM})$
		\STATE Algorithm \ref{alg:bb1} (compute state changes due to block $B_1$)
		\STATE Initialize internal clocks $R, P$ if needed (see \cite{ourSL,ourML})
		\STATE $\Delta \bar{X} \leftarrow 0$; $\Delta \dbar{X} \leftarrow 0$
		\FOR {$\mathcal{B} = B_2,B_3,B_4$}
		\STATE $t_r \leftarrow t$
		\STATE $\bar{X}_r \leftarrow \bar{X}$; $\dbar{X}_r \leftarrow \dbar{X}$
		\WHILE{$t_r<H$}
		\STATE update $P_{j \in \mathcal{B}}$
		\SWITCH {$\mathcal{B}$}
			\CASE {$B_2$}
			 	\STATE $\bar{d} \leftarrow \bar{a}_{j\in \mathcal{B}}$
				\STATE $\dbar{d} \leftarrow a_{j\in \mathcal{B}}(\dbar{X})$
				\STATE $\tau_r \leftarrow $ Compute the Chernoff tau-leap step size using $(\bar{X}_r,\bar{a}_{j\in \mathcal{B}},H,\bar{\delta})$
			\ENDCASE
			\CASE {$B_3$}
				\STATE $\bar{d} \leftarrow a_{j\in \mathcal{B}}(\bar{X})$
				\STATE $\dbar{d} \leftarrow \dbar{a}_{j\in \mathcal{B}}$
				\STATE $\tau_r \leftarrow $ Compute the Chernoff tau-leap step size using $(\dbar{X}_r,\dbar{a}_{j\in \mathcal{B}},H,\dbar{\delta})$
			\ENDCASE
			\CASE {$B_4$}
				\STATE $\bar{d} \leftarrow a_{j\in \mathcal{B}}(\bar{X})$
				\STATE $\dbar{d} \leftarrow a_{j\in \mathcal{B}}(\dbar{X})$
				\STATE $\tau_r \leftarrow \infty$
			\ENDCASE
		\ENDSWITCH
%			\STATE $A \leftarrow$ Algorithm \ref{alg:auxilco2} with ($\bar{d},\dbar{d}$)
			\STATE $A_1 \leftarrow \min(\bar{d},\dbar{d})$
			\STATE $A_2 \leftarrow \bar{d}-A_1$; $A_3 \leftarrow \dbar{d} -A_1$
			\STATE $H_r \leftarrow \min\{H,t_r{+}\tau_r\}$
%			\WHILE {$t_r < H_r$}
				\STATE $(t_r, \bar{X}_r, \dbar{X}_r, R_{j \mathcal{B}}, P_{j\in \mathcal{B}}) \leftarrow$ Algorithm \ref{alg:auxilco} with $(t_r,H_r,\bar{X}_r,\dbar{X}_r,R_{j\in \mathcal{B}},P_{j\in \mathcal{B}},A)$
%				\ENDWHILE
		\ENDWHILE
			\STATE $\Delta \bar{X} \leftarrow \Delta \bar{X} + (\bar{X}_r{-}\bar{X})$; $\Delta \dbar{X} \leftarrow \Delta \dbar{X} + (\dbar{X}_r{-}\dbar{X})$
		\ENDFOR
		\STATE $\bar{X} \leftarrow \bar{X} + \Delta \hat{X} + \Delta \bar{X}$;  $\dbar{X} \leftarrow \dbar{X} + \Delta \hat{\hat{X}} + \Delta \dbar{X}$
		\STATE $t \leftarrow H$
			\IF{$t<T$}% \AND \NOT out}
	\IF{$\bar H \leq \dbar H$}
		\STATE $\bar{t} \leftarrow$ next grid point in $\seqof{t_i}{i=0}{N}$ larger than $t$ 
\STATE $(\bar{H},\bar{\mathcal{R}}_{\TL},\bar{\mathcal{R}}_{\MNRM},\bar{a})  \leftarrow$ Algorithm \ref{alg:timehor} with ($\bar{X}$,$t$,$\bar{t}$,$T$,$\bar{\delta}$)

	\ENDIF
	\IF {$\bar H \geq \dbar H$}
		\STATE $\dbar{t} \leftarrow$ next grid point in $\seqof{s_j}{j=0}{N'}$ larger than $t$ 
\STATE $(\dbar{H},\dbar{\mathcal{R}}_{\TL},\dbar{\mathcal{R}}_{\MNRM},\dbar{a})  \leftarrow$ Algorithm \ref{alg:timehor} with ($\dbar{X}$,$t$,$\dbar{t}$,$T$,$\dbar{\delta}$)
	\ENDIF
\ENDIF
	\ENDWHILE
\end{small}
\end{algorithmic}
\end{algorithm}

\begin{algorithm}
\caption{\small Compute the next time horizon. Inputs: the current state, $\tilde{X}$, the current time, $t$, the next grid point, $\tilde{t}$, the final time, $T$, the one step exit probability bound, $\tilde{\delta}$, and the propensity functions, $a{=}\seqof{a_j}{j=1}{J}$.
Outputs: the next horizon $H$, the set of reaction channels to which the Tau-leap method should be applied, $\tilde{\mathcal{R}}_{\TL}$, the set of reaction channels to which MNRM should be applied, $\tilde{\mathcal{R}}_{\MNRM}$, and current propensity values $\tilde{a}$.}
\label{alg:timehor}
\begin{algorithmic}[1]
\begin{small}
\STATE $\tilde{a} \leftarrow a(\tilde{X})$
%\STATE $\tilde{t} \leftarrow$ next grid point in $\seqof{t_i}{i=0}{N}$ larger than $t$ 
\STATE $(\tilde{\mathcal{R}}_{\TL},\tilde{\mathcal{R}}_{\MNRM}) \leftarrow$ Algorithm \ref{alg:sel} with $(\bar X,t,\seqof{a_j(\bar X)}{j=1}{J},\tilde \delta,\tilde t,\kappa)$
%split reaction chans given ($\tilde{X}$,$t$,$\tilde{t}$,$T$,$\tilde{\delta}$)
\IF{$\tilde{\mathcal{R}}_{\TL} \neq \emptyset$}
\STATE $\tilde{H} \leftarrow \min\{\tilde{t},t{+}\tau(\tilde{\mathcal{R}}_{\TL}),T\}$
\ELSE
\STATE $\tilde{H} \leftarrow \min\{t{+}\tau(\tilde{\mathcal{R}}_{\TL}),T\}$
\ENDIF
\RETURN $(\tilde H,\tilde{\mathcal{R}}_{\TL},\tilde{\mathcal{R}}_{\MNRM},\tilde a)$
\end{small}
\end{algorithmic}
\end{algorithm}

\begin{algorithm}[h!]
\caption{\small Compute building block 1. This algorithm is part of Algorithm \ref{alg:coupled}.
}
\label{alg:bb1}
\begin{algorithmic}[1]
\begin{small}
		\STATE $t_r \leftarrow t$
		\STATE $\Delta \hat{\hat{X}} \leftarrow 0$; $\Delta \hat{X} \leftarrow 0$
		\WHILE {$t_r < H$}
			\STATE $\bar{\tau}_r \leftarrow $ Compute the Chernoff tau-leap step size using $(\bar{X}{+}\Delta \hat{X},\bar{a}_{j\in B_1},H,\bar{\delta})$
			\STATE $\dbar{\tau}_r \leftarrow $ Compute the Chernoff tau-leap step size using $(\dbar{X}{+}\Delta \hat{\hat{X}},\dbar{a}_{j \in B_1},H,\dbar{\delta})$
			%\STATE $\dbar{T}_r \leftarrow \min(\dbar{T},t_r+\bar{\tau}_r)$
			%\STATE $\dbar{T}_r \leftarrow \min(\dbar{T}_r,t_r+\dbar{\tau}_r)$
			\STATE $H_r \leftarrow \min\{H,t_r{+}\bar{\tau}_r,t_r{+}\dbar{\tau}_r\}$
%			\STATE $A \leftarrow$ Algorithm \ref{alg:auxilco2} with ($\bar{a}_{j \in B_1},\dbar{a}_{j \in B_1}$)
			\STATE $A_1 \leftarrow \min(\bar{a}_{j \in B_1},\dbar{a}_{j \in B_1})$
			\STATE $A_2 \leftarrow \bar{a}_{j \in B_1}-A_1$
			\STATE $A_3 \leftarrow \dbar{a}_{j \in B_1} -A_1$
			\STATE $\Lambda \leftarrow \mathcal{P}(A {\cdot} (H_r{-}t_r))$
			\STATE $\Delta \hat{\hat{X}} \leftarrow \Delta \hat{\hat{X}}+ (\Lambda_1 {+} \Lambda_2)\nu_{j \in B_1}$
			\STATE $\Delta \hat{X} \leftarrow \Delta \hat{X}+ (\Lambda_1 {+} \Lambda_3)\nu_{j \in B_1}$
			\STATE $t_r \leftarrow H_r$
%		\STATE $t \leftarrow \dbar{T}$
%		\STATE do TL-TL step
		\ENDWHILE
\end{small}
\end{algorithmic}
\end{algorithm}

\begin{algorithm}[h!]
\caption{\small The auxiliary function used in algorithm \ref{alg:coupled}. 
Inputs: current time, $t$, current time horizon, $\dbar{T}$, current system state at coarser level and finer level, $\bar{X}$, $\dbar{X}$, respectively, the internal clocks $R$ and $P$, the values, $A$, and the current building block, $B$. Outputs: updated time, $t$, updated system states, $\bar{X}$, $\dbar{X}$, and updated internal clocks, $R_i$, $P_i$, $i{=}1,2,3$. 
}
\label{alg:auxilco}
\begin{algorithmic}[1]
\begin{small}
\STATE $\Delta t_i \leftarrow (P_i - R_i)/A_i$, for $i=1,2,3$
\STATE $\Delta \leftarrow \min_i\{ \Delta t_i \}$
\STATE ${\mu} \leftarrow \text{argmin}_i \{\Delta t_i\}$% and  $\Delta_{\mu}$ the time where the minimum is realized.
					\IF{$t+\Delta > \dbar{T}$}
						\STATE $R \leftarrow R + A {\cdot} (\dbar{T}{-}t)$
						\STATE $t \leftarrow \dbar{T}$
					\ELSE
						\STATE update $\bar{X}$ and $\dbar{X}$ using $\nu_{j \in B}$
						\STATE $R \leftarrow R + A \Delta$
						\STATE $r \leftarrow$ uniform$(0,1)$
						\STATE $P_\mu \leftarrow P_\mu + \log (1/r)$
						\STATE $t \leftarrow t + \Delta$
					\ENDIF
					\RETURN $(t, \bar X, \dbar X, R, P)$
\end{small}
\end{algorithmic}
\end{algorithm}

%\begin{algorithm}[h!]
%\caption{\small Auxiliary function used in algorithm \ref{alg:coupled}. Inputs: the propensity values at the coarse and fine grid, $\bar{a},\dbar{a}$. Output: $A$.
%}
%\label{alg:auxilco2}
%\begin{algorithmic}[1]
%\begin{small}
%			\STATE $A_1 \leftarrow \min(\bar{a},\dbar{a})$
%			\STATE $A_2 \leftarrow \bar{a}-A_1$, $A_3 \leftarrow \dbar{a} -A_1$
%			\RETURN $A$
%\end{small}
%\end{algorithmic}
%\end{algorithm}
%\begin{algorithm}[h!]
%\caption{Auxiliary function used in algorithm \ref{alg:coupled}. Inputs: current state at fine grid, $\dbar{X}$,  next time horizon at coarse grid, $\bar{T}$, and the fine grid,$\seqof{s_j}{j=0}{N'}$. Output: next time horizon at fine grid, $\dbar{T}$, and current propensity values at fine grid, $\dbar{a}$.
%}
%\label{alg:auxilco3}
%\begin{algorithmic}[1]
%			\STATE $\dbar{a} \leftarrow a(\dbar{X})$
%			\STATE $\dbar{\tau} \leftarrow$ Algorithm \ref{alg:Ch} with ($\dbar{X}$,$\dbar{a}$)
%			\STATE $\dbar{T} \leftarrow \min(\bar{T},t+\dbar{\tau})$
%			\STATE $s \leftarrow$ next point in the finer grid $\seqof{s_j}{j=0}{N'}$
%			\STATE $\dbar{T} \leftarrow \min(\dbar{T},s)$
%			\RETURN ($\dbar{T}$,$\dbar{a}$)
%\end{algorithmic}
%\end{algorithm}
%\end{document}


\section{The Multilevel Estimator and Total Error Decomposition}
\label{sec:mlmcerr}
In this section, we first show the multilevel Monte Carlo estimator. We then  analyze and control the computational global error, which is decomposed into three error components: the discretization error, the global exit error, and the Monte Carlo statistical error. Upper bounds for each one of the three components are given. 
Finally, we briefly describe the automatic estimation procedure that allows us to estimate our quantity of interest within a given prescribed relative tolerance, up to a given confidence level.

%we present a control variate for $X(T)$, which in some cases can obtain substantial gains in computational work. 

\subsection{The MLMC Estimator}\label{sec:mlmc}
In this section, we discuss and implement a multilevel Monte Carlo estimator for the mixed Chernoff tau-leap case. 
%The main ingredient of this section is Algorithm \ref{alg:coupled}, which generates coupled hybrid paths at levels $\ell{-}1$ and $\ell$.
%Let us now introduce some notation. 
Consider a hierarchy of nested meshes of the time interval $[0,T]$, indexed by $\ell = 0,1,\ldots,L$.
Let $\Delta t_0$ be the size of the coarsest time mesh that corresponds to the level $\ell{=}0$.
The size of the time mesh at  level $\ell \geq 1$ is given by $\Delta t_{\ell} {=} R^{-\ell} \Delta t_0$, where $R{>}1$ is a given integer constant. 
Let $\{\bar X_{\ell}(t)\}_{t\in[0,T]}$ be a mixed Chernoff tau-leap process with a time mesh of size $\Delta t_{\ell}$ and a one-step exit probability bound $\delta$, and let $g_{\ell}{:=} g({\bar X}_{\ell}(T))$ be our quantity of interest computed with a mesh of size $\Delta t_{\ell}$. We can simulate paths of  $\{\bar X_{\ell}(t)\}_{t\in[0,T]}$ by using  Algorithm \ref{alg:mixed}.
We are interested in estimating $\expt{g_{L}}$, and we can simulate correlated pairs, $(g_{\ell},g_{\ell{-}1})$ for $\ell=1,\ldots,L$, by using Algorithm \ref{alg:coupled}.
Let $A_\ell$ be the event in which the $\ell$-th grid level path, $\bar X_{{\ell}}$, arrives at the final time, $T$, without exiting the state space of $X$. %Let 
%Finally, $g_{\ell}:= g({\bar X}_{\ell}(T))$ was defined in Section \ref{MLMCintro}.

Consider the following telescopic decomposition:
\begin{align*}
\expt{g_L \indicator{A_L}} = \expt{g_0 \indicator{A_0}} {+} \sum_{\ell=1}^L \expt{g_\ell\indicator{A_\ell}-g_{\ell-1}\indicator{A_{\ell-1}}} \COMMA
\end{align*} 
where $\indicator{A}$ is the indicator function of the set $A$. This motivates the definition of our MLMC estimator of $\expt{g(X(T))}$:
\begin{align}\label{MLMCest}
\mathcal{M}_L := \frac{1}{M_0} \sum_{m=1}^{M_0}g_0\indicator{A_0}(\omega_{m,0}) + \sum_{\ell=1}^L \frac{1}{M_{\ell}} \sum_{m=1}^{M_{\ell}} [g_{\ell} \indicator{A_\ell} - g_{\ell-1} \indicator{A_{\ell-1}}](\omega_{m,\ell}) \PERIOD
\end{align}
\subsubsection*{Computational Complexity} 
\label{sec:complex}
A key property of our multilevel estimator is that the computational work is a function of the given relative tolerance, $TOL$, is of the order of $TOL^{-2}$.
 The optimal work is is given by 
\begin{equation*}
w^*_{L} (TOL) = \LP\frac{C_A}{\theta}\sum_{\ell=0}^L \sqrt{\mV{\ell} \psi_{\ell}}\RP^{2} TOL^{-2}\PERIOD
\end{equation*}
From the fact that the sum $\sum_{\ell=0}^{\infty} \sqrt{\mV{\ell} \psi_{\ell}}$ converges, because $\psi_\ell = \Ordo{\psi_{\MNRM}}$, 
%First, observe that the expected computational work per path at level $\ell$, $\psi_{\ell}$, is bounded by a multiple of the expected computational work of the MNRM (see Section \ref{sec:MNR}), \ie, $K \psi_{\text{MNRM}}$. In our numerical experiments, we observe that $K$ around three is enough. Therefore, $\sum_{\ell=0}^{\infty} \sqrt{\mV{\ell} \psi_{\ell}} \leq \sqrt{K \psi_{\text{MNRM}}} \sum_{\ell=0}^{\infty}  \sqrt{\mV{\ell}}$.Since $\mV{\ell}=\Ordo{\Delta t_\ell} \leq C \Delta t_0(1/2)^{\ell}$ for some positive constant $C$, the series $\sum_{\ell =0}^{\infty} \mV{\ell} ^{1/2}$ converges like the geometric series $\sum_{\ell =0}^{\infty} (1/\sqrt{2})^{\ell}$ . 
we conclude that 
$\sup_{L}\{\sum_{\ell=0}^L \sqrt{\mV{\ell} \psi_{\ell}}\}$ is bounded and, therefore, the expected computational complexity of the multilevel mixed Chernoff tau-leap method is $w^*_L(TOL){=}\Ordo{TOL^{-2}}$. 

\subsection{Global Error Decomposition}
In this section, we define the computational global error, $\EE_L$, and show how it can be naturally decomposed into three components: the discretization error, $\EE_{I,L}$, and the exit error, $\EE_{E,L}$, both coming from the tau-leap part of the mixed method, and the Monte Carlo statistical error, $\EE_{S,L}$. %Next, we show how to model and control the global error, $\EE_L$, 
We also give upper bounds for each one of the three components. 

The computational global error, $\EE_L$, is defined as
\begin{equation*}%\label{eq:global error}
\EE_L := \expt{g(X(T))} - \mathcal{M}_L \COMMA
\end{equation*}
and can be decomposed as
\begin{align*}%\label{eq:globalerror}
\expt{g(X(T))} - \mathcal{M}_L &= \expt{g(X(T))(\indicator{A_L}+\indicator{A^c_L})} \pm \expt{g_L \indicator{A_L}} - \mathcal{M}_L\nonumber \\
&= \underbrace{\expt{g(X(T)) \indicator{A^c_L}}}_{=:\EE_{E,L}} 
 + \underbrace{\expt{ \LP g(X(T)) {-} g_L \RP \indicator{A_L}}}_{=:\EE_{I,L}} + \underbrace{\expt{g_L \indicator{A_L}}{-}\mathcal{M}_L}_{=:\EE_{S,L}}\PERIOD 
\end{align*}

We showed in \cite{ourSL} that by choosing adequately the one-step exit probability bound, $\delta$, the exit error, $\EE_{E,L}$, satisfies 
$|\EE_{E,L}| \leq |\expt{g(X(T))}|\, \prob{A^c_L}\leq TOL^2$. 
%observing that we are free to choose $X$ and $\bar X$ independent.

An efficient procedure for accurately estimating $\EE_{I,L}$ in the context of the tau-leap method is described in \cite{ourML}. 
%We adapt this method in Algorithm \ref{alg:weakerror} for estimating the weak error in the hybrid context. A brief description follows. 
For each mixed path,  
$\seqof{\bar X_{\ell}(t_{n,\ell},\bar \omega)}{n=0}{N(\bar{\omega})}$, we define  the sequence of dual weights, $\seqof{\varphi_{n,\ell}(\bar \omega)}{n=1}{N(\bar{\omega})}$, backwards as follows:
\begin{align}\label{eq:phisback}
\varphi_{N(\bar{\omega}),\ell} &:= \nabla g(\bar X_{\ell}(t_{N(\bar{\omega}),\ell},\bar \omega))\\
\varphi_{n,\ell} &:= \LP Id +  {\Delta t_{n,\ell}} \JAC_a^T(\bar X_{\ell}(t_{n,\ell},\bar \omega)) \,\nu^T\RP \varphi_{n+1,\ell}, \quad n=N(\bar{\omega}){-}1,\ldots,1\COMMA\nonumber
\end{align}
where $\Delta t_{n,\ell}{:=}t_{n+1,\ell}{-}t_{n,\ell}$, $\nabla$ is the gradient operator and $\JAC_a(\bar X_{\ell}(t_{n,\ell},\bar \omega)){\equiv}[\partial_i a_j(\bar X_{\ell}(t_{n,\ell},\bar \omega))]_{j,i}$ is the Jacobian matrix of the propensity function, $a_j$, for $j{=}1\ldots J$ and $i{=}1\ldots d$. 
%This procedure computes  $\EE_{I,L}(\bar \omega)$ for every simulated hybrid trajectory, $\seqof{\bar X_{L}(t_n,\bar \omega)}{n=0}{N(\bar{\omega})}$, as a weighted sum of local errors at the mesh times, $\seqof{t_n}{n=0}{{N(\bar{\omega})}}$. 
We then approximate $\EE_{I,L}$ by $\avg{\EE_{I,L}(\bar\omega)}{\cdot}$, where 
\begin{equation*}%\label{eq:EILw}
\EE_{I,L}(\bar \omega) := \sum_{n=1}^{{N(\bar{\omega})}}  \LP\frac {\Delta t_{n,L}}{2} \varphi_{n,L}   \sj \indicator{j\in \mathcal{R}_{\TL}(n)} \nu_j^T \LP a_j(\bar X_{L}(t_{n+1,\ell})){-}a_j(\bar X_{L}(t_{n,\ell}))\RP \RP(\bar \omega)\COMMA
\end{equation*}
$\avg{X}{M}{:=}\frac 1 M \sum_{m=1}^M X(\omega_m)$ and $\svar{X}{M}{:=}\avg{X^2}{M}- \avg{X}{M}^2$ denote
the sample mean and the sample variance of the random  variable, $X$, respectively.
Here $\indicator{j\in \mathcal{R}_{\TL}(n)}{=}1$ if and only if, at time $t_{n,\ell}$, the tau-leap method was used for reaction channel $j$, and we denote by $Id$ the $d\times d$ identity matrix. 

The variance of the statistical error, $\EE_{S,L}$, is given by $\sum_{\ell=0}^L \frac{\mV{\ell}}{M_{\ell}}$,
where $\mV{0}:=\var{g_{0}\indicator{A_0}}$ and 
$\mV{\ell}:=\var{ g_{\ell} \indicator{A_{\ell}} -g_{\ell-1} \indicator{A_{\ell{-}1}}},\,\,\ell \geq 1$. 
In \cite{ourML}, we presented an efficient and accurate method for estimating $\mV{\ell}$, $\ell \geq 1$ using the formula
\begin{align*}
%\label{eq:varhat}
\hat{\mathcal{V}}_\ell&:=\svar{\sum_n \expt{ \varphi_{n+1} \cdot e_{n+1}\SEP\mathcal{F}}(\bar \omega)}{M_{\ell}}+ \avg{\sum_n \var{ \varphi_{n+1} \cdot e_{n+1}\SEP\mathcal{F}}(\bar \omega)}{M_{\ell}}\COMMA
\end{align*}
where $\mathcal{F}$ is a suitable chosen sigma algebra such that $\seqof{\varphi_n(\bar \omega)}{n=1}{N(\bar \omega)}$ is measurable, with $N(\bar \omega)$ being the total number of steps given by Algorithm \ref{alg:coupled}. In this way, the only randomness in $\expt{ \varphi_{n+1} \cdot e_{n+1}\SEP\mathcal{F}}$ and $\var{ \varphi_{n+1} \cdot e_{n+1}\SEP\mathcal{F}}$ comes from the local errors, $\seqof{e_n}{n=1}{N(\bar \omega)}$, defined as $e_{n} := X_{\ell,n} - X_{\ell-1,n}$.
In the aforementioned work, we derived exact and approximate formulas for computing $\expt{ \varphi_{n+1} \cdot e_{n+1}\SEP\mathcal{F}}$ and $\var{ \varphi_{n+1} \cdot e_{n+1}\SEP\mathcal{F}}$.
\begin{rem}[Backward Euler] In \eqref{eq:phisback}, we have that $\varphi_{n,\ell}$ can be computed by a backward Euler formula when too fine time meshes are required for stability, \ie,
$\varphi_{n,\ell} := \LP Id -  {\Delta t_{n,\ell}} \JAC_a^T(\bar X_{\ell}(t_{n,\ell},\bar \omega)) \,\nu^T\RP^{-1} \varphi_{n+1,\ell}$.   
\end{rem}

\subsection{Estimation Procedure}
\label{sec:estim}
In this section, we briefly describe the automatic procedure that estimates $\expt{g(X(T))}$ within a given prescribed relative tolerance, $TOL{>}0$, up to a given confidence level. Up to minor changes, it is the same as the one presented in \cite{ourML}. It is important to remark that the minimal user intervention is required to obtain the parameters needed to simulate the mixed paths, and subsequently, to compute the estimations using \eqref{MLMCest}. 
Once the reaction network is given (stoichiometric matrix $\nu$ and $J$ propensity functions $a_j$), the user only needs to set the required maximum allowed relative global error or tolerance, $TOL$, and the confidence level, $\alpha$. This process has three phases:
\begin{description}
\item[Phase I] Calibration of virtual machine-dependent quantities. In this phase, we estimate the quantities 
$C_\MNRM$, $C_\TL$, $C_s$ and the function $C_P$
that allow us to model the expected computational work, measured in runtime. %This are defined in Section \ref{sec:splitting}.
\item[Phase II] Solution of the work optimization problem:  we obtain 
the total number of levels, $L$, and the sequences, $\seqof{\delta_\ell}{\ell=0}{L}$ and $\seqof{M_\ell}{\ell=0}{L}$, \ie, the one-step exit probability bounds and the required number of simulations at each level. 
In this phase, given a relative tolerance, $TOL{>}0$, we solve the work optimization problem
\begin{align}
\left\{ \!
\begin{array}{l}\label{eq:minworkprob}
\min_{\{\Delta t_0, L,\seqof{M_{\ell},\delta_{\ell}}{\ell=0}{L}\}} \sum_{\ell =0}^{L}\psi_\ell M_{\ell}\\
\mbox{s.t.}\\
\EE_{E,L} +\EE_{I,L}+\EE_{S,L} \leq TOL
\end{array}
\right. \PERIOD
\end{align}
An algorithm to efficiently compute the solution of this optimization problem is given in \cite{ourML}.
Our objective function is the expected total work of the MLMC estimator, $\mathcal{M}_L$, \ie,
$\sum_{\ell =0}^{L}\psi_\ell M_{\ell}\COMMA$
where $L$ is the deepest level, $\psi_0$ is the expected work of a single-level path at level 0, and $\psi_\ell$, for $\ell\geq 1$, is the expected computational work of two coupled paths at levels $\ell{-1}$ and $\ell$. Finally, $M_0$ is the number of single-level paths at level 0, and $M_\ell$, for $\ell\geq 1$, is the number of coupled paths at levels $\ell{-1}$ and $\ell$. 
We now describe the quantities $\seqof{\psi_\ell}{\ell=0}{L}$. 
First, $\psi_0$ is the expected work of a single hybrid path (simulated by Algorithm \ref{alg:mixed}), %The main building block is Algorithm \ref{alg:sel}, which adaptively determines the sets $\mathcal{R}_\TL$ and $\mathcal{R}_\MNRM$ at each decision step. 
%The expected work of a single hybrid path, at level $\ell=0$, is
\begin{align}\label{eq:cost}
\psi_0 &:= C_\MNRM \expt{\NMNRM(\Delta t_0,\delta_0)} + C_\TL \expt{\NTL (\Delta t_0,\delta_0)}\\ &+  \int_{[0,T]}\expt{ \sum_{j \in \mathcal{R}_{\TL}(s)} C_P(a_j(\bar X_0(s))\tau_{Ch}(\bar X_0(s),\delta_0)) ds} \nonumber
\COMMA
\end{align} 
where $\Delta t_0$ is the size of the time mesh at level 0 and $\delta_0$ is the exit probability bound at level 0, and $\mathcal{R}_\TL=\mathcal{R}_\TL(t)$ is the tau-leap set, which depends on time (and also the current state of the process).
The set $\mathcal{R}_\TL$ is determined  at each decision step by Algorithm \ref{alg:sel}.
Therefore, the expected work at level 0 is $\psi_0 M_0$, 
where $M_0$ is the total number of single hybrid paths. 

For $\ell \geq 1$, we use Algorithm \ref{alg:coupled} to generate $M_{\ell}$-coupled paths that couple levels $\ell{-}1$ and $\ell$.
%For that reason, we have to consider two different time meshes of sizes $h_{\ell{-}1}$ and $h_{\ell}$, respectively, but just one delta, $\delta_{\ell}$ (see Remark \ref{rem:whybother}). 
%Given two coupled paths, let $\NMNRMKone(\Delta t_{\ell{-}1},\delta_{\ell-1})$ and $\NMNRMKone(\Delta t_{\ell},\delta_{\ell})$ be the number of exact steps for level ${\ell{-}1}$ (coarse mesh) and $\ell$ (fine mesh), respectively, with associated work $C_1$. We define $\NMNRMKtwo(\Delta t_{\ell{-}1},\delta_{\ell-1})$ and $\NMNRMKtwo(\Delta t_{\ell},\delta_{\ell})$ analogously. 
%Finally, let $\NTLC(h_{\ell{-}1},\delta_{\ell}-1)$, $\NTLC(h_{\ell},\delta_{\ell})$ be the number of tau-leap steps. 
The expected work of a pair of coupled hybrid paths at levels $\ell$ and $\ell-1$ is
\begin{align}\label{eq:costc}
\psi_\ell &:= C_\MNRM \expt{\NMNRMC(\ell)} +C_\TL \expt{\NTLC(\ell)} \\ \nonumber
&+  \int_{[0,T]} \expt{\sum_{j\in \mathcal{R}_{\TL,\ell}(s)} C_P(a_j(\bar X_\ell(s))\tau_{Ch}(\bar X_\ell(s),\delta_\ell)) ds}\\ \nonumber
 &+ \int_{[0,T]} \expt{\sum_{j\in \mathcal{R}_{\TL,\ell-1}(s)}  C_P(a_j(\bar X_{\ell-1}(s))\tau_{Ch}(\bar X_{\ell-1}(s),\delta_{\ell-1})) ds}\COMMA
\end{align} 
where 
\begin{align*}
&\NMNRMC(\ell):=\NMNRM(\Delta t_\ell,\delta_\ell) + \NMNRM(\Delta t_{\ell-1},\delta_{\ell-1})\\
&\NTLC(\ell) :=\NTL(\Delta t_\ell,\delta_\ell)+\NTL (\Delta t_{\ell-1},\delta_{\ell-1})\PERIOD
\end{align*} 
\item[Phase III] Estimation of $\expt{g(X(T))}$. 
\end{description} 

%In the Virus example $g$ is just the projection of a component.
%\subsection{Phase I}
%\subsection{Phase II}
%\subsubsection{Bayesian updating}
%To solve the family of optimization problems \eqref{eq:auxminprob}, it is required that we previously know the sequences  
%$\seqof{\EE_{I,\ell}}{\ell=0}{L}$, $\seqof{\mV{\ell}}{\ell=0}{L}$, and $\seqof{\psi_\ell}{\ell=0}{L}$.
%Notice that  $\mV{0}$ and $\psi_0$ are related to single paths, but for $\ell\geq 1$, 
%$\mV{\ell}$ and $\psi_\ell$ are related to coupled paths.
%We assume for now that we know, for every $\ell$, the value $\delta_\ell$ such that 
%$\EE_{E,\ell}{<}TOL^2$. 
%
% To estimate $\EE_{I,0}$, $\mV{0}$ and $\psi_0$, we just sample single-level paths with parameters $(h_0,\delta_0)$ and use a standard Monte Carlo estimator. The number of samples should be enough for controlling the coefficients of variation of the mentioned estimators.  In principle, we could follow the same strategy to estimate $\EE_{I,\ell}$, $\mV{\ell}$ and $\psi_\ell$, for $\ell{\geq}1$. But notice that for large values of $\ell$, not only is $\psi_\ell$ large, but a large number of sampled coupled hybrid paths are also needed to estimate those quantities with a small coefficient of variation. 
%To overcome this problem, we apply a Bayesian approach that uses the information obtained from sampling at levels $\ell{-}2$ and $\ell{-}1$ to establish prior distributions for $\EE_{I,\ell}$, $\mV{\ell}$ and $\psi_\ell$. In this way, just  by drawing few coupled paths at levels $\ell{-}1$ and $\ell$, we can update our knowledge about those three mentioned quantities.
%This estimation strategy turns out to be very efficient in terms of  work, and also reasonably accurate.
%This is mainly due to the regularity of those quantities as functions of the level $\ell$, for instance, see Figure \ref{fig:dec2-diag}. 
%
%Let us consider with some detail the case of $\mV{\ell}$. 
%There are two sources of prior information to estimate $\mV{\ell}$: i) extrapolating linearly the values 
%$\log{\hmV{\ell-1}}$ and $\log{\hmV{\ell-2}}$, and then computing 
%$\tilde{\mathcal{V}}_{\ell}:= \hmV{\ell-1}^2/\hmV{\ell-2}$, 
%ii) sampling pairs $\LP \bar X_{\ell-2}, \bar X_{\ell-1}\RP$, and then using $\bar X_{\ell-1}$ in the dual formulas \eqref{eq:Ephinu}, \eqref{eq:Vphinu} and \eqref{eq:varhat} to obtain $\hat{\mathcal{V}}_{\ell,\delta_{\ell-1}}$. Of course, this last estimate corresponds to $\delta_{\ell-1}$ and not to $\delta_{\ell}$, but, still, it is useful prior information. 
%
%In our calculations, we use a lognormal model for the prior distribution of $\mV{\ell}$, where its parameters are estimated by combining  $\tilde{\mathcal{V}}_{\ell}$ and $\hat{\mathcal{V}}_{\ell,\delta_{\ell-1}}$.
%Using conjugate analysis (\cite{Fink97, Gelman}), we update this prior knowledge by sampling coupled paths $\LP \bar X_{\ell-1}, \bar X_{\ell}\RP$ and use $\bar X_{\ell-1}$ to compute dual weighted estimates of $\mV{\ell}$.
%Note that we do not use the observations $\Delta g_{\ell}:= \LP g_{\ell}-g_{\ell-1}\RP \indicator{A_{\ell{-}1}} \indicator{A_{\ell}}$ due to its high variance. 
%Figure \ref{fig:arrows} depicts the Bayesian updating process for the estimation of ${\mathcal{V}}_{\ell}$.
%\begin{figure}[h!]
%\begin{displaymath}
%    \xymatrix{ (\bar X_{\ell-2}, \blue{\bar X_{\ell-1}})_{\delta_{\ell-1}}  \ar@{~>}[dr] \ar[dd] &  &\\
%                   & \blue{\hat{\mathcal{V}}_{\ell,\delta_{\ell-1}}} \&\,\, \tilde{\mathcal{V}}_{\ell}\,\,  (\text{prior}) \ar@{.>}[d] &\\
%              (\blue{\bar X_{\ell-1}}, \bar X_{\ell})_{\delta_{\ell}}\ar@{~>}[r]  &  \hat{\mathcal{V}}_{\ell,\delta_{\ell}} }
%\end{displaymath}
%\caption{Scheme for the Bayesian updating of $\hat{\mathcal{V}}_{\ell,\delta}$.}
%\label{fig:arrows}
%\end{figure}
%We proceed analogously with the quantities  $\EE_{I,\ell}$, and $\psi_\ell$ 

%\subsubsection{Some Comments on the Algorithms for Phase II}
%
%In Algorithm \ref{alg:Cal}, we propose an iterative method to obtain an approximate solution to the problem
%\eqref{eq:minworkprob}. Notice that we are assuming that there are at least two levels in the multilevel hierarchy, \ie, $L\geq 1$. 
%%Given $L$, to solve the problem  \eqref{eq:minworkprob}, we need accurate estimations of $\psi_\ell$, $\mV{\ell}$ and $\EE_{I,\ell}$, for $\ell=0,1,\ldots,L$  with minimal amount of  work. 
%%We control that  work as follows: every time we simulate a sample path (single hybrid or coupled hybrid), we decrement a certain quantity, $rt$, that we consider as our computational budget. 
%%If we set $rt{=}\infty$, we can simulate as many paths as needed, until the relevant coefficients of variation are small enough, for instance, the coefficient of variation of $\svar{g(\bar X(T))}{\cdot}$, $\svar{g_{\ell}{-}g_{\ell-1}}{\cdot}$ and $\avg{\EE_{I,L}}{\cdot}$. 
%%Typically, we set $rt{=}\infty$ for $L\leq 2$, but when we go to deepest levels, and consequently $\psi_{\ell}$ gets larger, it is not convenient to set $rt{=}\infty$, and we have to rely on prior knowledge.
%%Our strategy for $L{>}2$ is to control the number of sample coupled paths either by letting the work at level $\ell{+}1$ be at most double the work at level $\ell$, or by solving, when possible, the work optimization problem for $L{+}1$ using prior information on level $\ell{+}1$, and sample at most $M^*_{\ell+1}$ coupled paths.  
%%The strategy of doubling the previous budget is to guarantee that the total work is less than or equal to double the work of the last level. 
%
%To solve the problem \eqref{eq:minworkprob}, we bound the global exit error, 
%$\EE_{E,L}$, by $TOL^2$. More specifically, we choose $\delta_L$ to be sufficiently small such that 
%\begin{equation}\label{eq:lastdelta}
%|\avg{g_{L}}{\cdot}|\,\, \delta_{L}\,\,  \avg{\NTL(\Delta t_L,\delta_L)}{\cdot} < TOL^2\PERIOD
%\end{equation} 
%At this point, it is crucial to observe that, condition \eqref{eq:lastdelta}, when $\ell{<}L$, artificially enforces a dependence on $TOL$ that may result in very small values of $\delta_\ell$. This could affect the expected number of exact steps and tau-leap steps at level $\ell$. In the appendix of \cite{ourSL}, we proved that, when $\delta_\ell$ tends to zero, the expected values of the number of tau-leap steps at level $\ell$ goes to zero, and therefore our hybrid MLMC strategy would converge to the SSA method without the desired reduction in the computational work. 
%In order to avoid the dependence of $\seqof{\delta_\ell}{\ell=0}{L-1}$ on $TOL$, we adopt a different strategy based on the following decomposition:
%\begin{align*}
%\mV{\ell} = \var{g_\ell \indicator {A_{\ell}} - g_{\ell-1} \indicator{A_{\ell-1}}} &= \var{g_\ell - g_{\ell-1} \SEP 
%A_{\ell}\cap A_{\ell-1}}\prob{A_{\ell}\cap A_{\ell-1}}\\
%&+ \var{g_\ell\SEP A_{\ell}\cap A^c_{\ell-1}}\prob{A_{\ell}\cap A^c_{\ell-1}}\\
%&+ \var{g_{\ell-1}\SEP A^c_{\ell}\cap A_{\ell-1}}\prob{A^c_{\ell}\cap A_{\ell-1}}\PERIOD
%\end{align*}   
%We impose that the first term of the right-hand side dominates the other two. This is because the conditional variances appearing in the last two terms are of order $\Ordo{1}$, while the conditional variance appearing in the first term is of order $\Ordo{\Delta t_{\ell}}$, and we make our computations with approximations of $\mV{\ell}$ assuming that $\prob{A_{\ell}\cap A_{\ell-1}}$ is close to one.  We proceed as follows: first, we approximate $\prob{A_{\ell}\cap A_{\ell-1}}$ by $\prob{A_{\ell}} \prob{A_{\ell-1}}$; then, we consider $1{-}\delta_\ell \avg{\NTL(\Delta t_\ell,\delta_\ell)}{\cdot}$ as an approximate upper bound for $\prob{A_{\ell}}$ when  $\delta_\ell \avg{\NTL(\Delta t_\ell,\delta_\ell)}{\cdot}{\ll} 1$. Those considerations lead us to impose 
%\begin{align}\label{eq:deltas}
%&\var{g_\ell - g_{\ell-1} \SEP A_{\ell}\cap A_{\ell-1}}
%\LP 1{-}\delta_\ell \avg{\NTL(\Delta t_\ell,\delta_\ell)}{\cdot}\RP
%\LP 1{-}\delta_{\ell-1} \avg{\NTL(\Delta t_{\ell-1},\delta_{\ell-1})}{\cdot}\RP>\\
%&\var{g_\ell\SEP {A_{\ell}}\cap{A^c_{\ell-1}}}\delta_{\ell-1} \avg{\NTL(\Delta t_{\ell-1},\delta_{\ell-1})}{\cdot} + 
%\var{g_{\ell-1}\SEP{A^c_{\ell}}\cap{A_{\ell-1}}}\delta_{\ell} \avg{\NTL(\Delta t_{\ell},\delta_{\ell})}{\cdot}\PERIOD\nonumber
%\end{align}
%To avoid simultaneous refinements on $\delta_\ell$ and $\delta_{\ell-1}$, based on \eqref{eq:deltas}, we impose to $\delta_\ell$ the following condition:
%\begin{equation*}%\label{eq:deltasapp}
%\hmV{\ell}\LP 1{-}\delta_{\ell} \avg{\NTL(\Delta t_{\ell},\delta_{\ell})}{\cdot}\RP^2 >
% 2\,\, \svar{g}{\cdot}\delta_{\ell} \avg{\NTL(\Delta t_{\ell},\delta_{\ell})}{\cdot} \PERIOD
%\end{equation*}
% 
%Algorithms \ref{alg:pathdualscostSL} and \ref{alg:Cal} provide  $\avg{g_{\ell}}{\cdot}$, $\avg{\NTL}{\cdot}$ and the other required quantities.
%
%%As a general strategy, we proceed as follows:
%%for $L{=}2$, we estimate  $\seqof{\EE_{I,\ell}}{\ell=0}{L}$, $\seqof{\mV{\ell}}{\ell=0}{L}$, and $\seqof{\psi_\ell}{\ell=0}{L}$, stopping by the coefficient of variation rule. 
%%If $w_2{\geq}w_1$, we set the optimal level $L^*{=}1$. Otherwise, based on the samples that we already have at levels $\ell{=}1,2$, we assign prior values to $\EE_{I,3}$, $\mV{3}$ and $\psi_3$. If $\EE_{I,3}{<}TOL_0$, we solve the problem \eqref{eq:greedyKKT} for $L{=}3$, and obtain a auxiliary values $M_{3,aux}^*$ and  estimated optimal cost $w_{3,aux}$. If, $w_{3,aux}{\geq} w_2$, we set the optimal level $L^*{=}2$, otherwise, we start sampling at levels $\ell{=}2,3$ to update our prior knowledge about $\EE_{I,3}$, $\mV{3}$ and $\psi_3$. We use at maximum $M_{3,aux}^*$, however this value could be re-estimated on the fly. 
%%But if $\EE_{I,3}{\geq}TOL_0$, we sample set $M_{3,aux}^* {:=} 2 \, \hat{\psi_2}/\psi_{3,pr} M_2$, where here $\psi_{3,pr}$ stands for our previous knowledge about $\psi_3$, and $M_2$ is the number of coupled sampled paths used at levels $\ell{=}1,2$.
%%We proceed iteratively. 
%
%%This estimations are obtained using the Algorithm \ref{alg:pathdualscost}. To this end, we simulate a certain number of single hybrid paths at level $\ell{=}0$ and coupled hybrid paths at levels $\ell {\geq} 1$ until either the cost reaches the budget or the coefficients of variation of both $\svar{g_{\ell}{-}g_{\ell-1}}{\cdot}$ and $\avg{\EE_{I,L}}{\cdot}$, are  below $cv_0$. 
%
%\begin{rem}[Multilevel estimators used in Algorithm \ref{alg:Cal}]
%\label{rem:algCal}
%Although in algorithm \ref{alg:Cal} we show that the estimations of $\expt{g(X(T))}$ and $\var{g(X(T))}$ are computed using the information of the last level only, in fact we are computing them using a multilevel estimator. We omit the details in the algorithm for the sake of simplicity. In the case of $\expt{g(\bar X(T))}$, we use the standard mutilevel estimator, and, for the case of $\var{g(X(T))}$, we use the following telescopic decomposition:
%\begin{align*}
%\var{g(\bar X_l(T))} = \var{g(\bar X_0(T))} + \sum_{\ell=1}^{l}(\var{g(\bar X_\ell(T))}-\var{g(\bar X_{\ell-1}(T))})\COMMA
%\end{align*}
%where $l>1$ is a fixed level. Using the usual variance estimators for each level, we obtain an unbiased multilevel estimator of the variance of $g(\bar X)$. We refer to \cite{Bierig.Chernov} for details.
%\end{rem}
%\begin{rem}[Coupled paths exiting the lattice, $\latt$]
%Algorithm \ref{alg:coupled} could compute four types of paths. It could happen that no approximate process (the coarse one, $X_{\ell-1}$, or the fine one, $X_\ell$) exits the lattice, which is the most common case. It could also happen that one of the approximate processes exits the lattice. And finally, both approximate processes could exit the lattice. The first case is the most common one and no further explanation is required. We now explain the case that one of process exits the lattice. Suppose the coarse one exited the lattice. In that case, until the fine process reaches time $T$ or exits the lattice, we still simulate the coupled process by simulating only the fine path using the single-level hybrid algorithm, presented in \cite{ourSL}. If the fine path reaches $T$, we have that $\indicator{A_{\ell-1}}=0$, and $\indicator{A_{\ell}}=1$. Vice versa, if the fine process exits and the coarse one reaches $T$, we have  $\indicator{A_{\ell-1}}=1$, and $\indicator{A_{\ell}}=0$.
%\end{rem}
%
%\begin{rem}[Coupling with an exact path]
%\label{rem:unbiased}
%Algorithm \ref{alg:Cal} uses a computational-cost-based stopping critera. That is, the algorithm stops refining the time mesh when the estimated total compuational cost of the multilevel estimator, $\hat W_{\ML}{:=}\sum_{\ell =0}^{l}\hat \psi_\ell M_{\ell}$, at level $l$, is greater than the corresponding cost for level $l{-}1$, and only when condition $\hat \WE {<} TOL{-}TOL^{2}$ is already satisfied. In that case, $L^*{=}l{-}1$. The latter condition is required for obtaining a solution of the optimization problem \eqref{eq:greedyKKT}. In our experiments, we observed that, the cost of two hybrid  coupled paths, $\psi_\ell$, may be greater than the cost of ``hybrid-exact'' coupled paths, that is, a hybrid path at level $l{-}1$ coupled with an exact path at level $l$. That kind of path, used only at the last level, leads to the following unbiased multilevel estimator:
%\begin{align*}
%\mathcal{\tilde M}_L &:= \frac{1}{M_0} \sum_{m=1}^{M_0}g_0\indicator{A_0}(\omega_{m,0}) + \sum_{\ell=1}^{L-1} \frac{1}{M_{\ell}} \sum_{m=1}^{M_{\ell}} [g_{\ell} \indicator{A_\ell} - g_{\ell-1} \indicator{A_{\ell-1}}](\omega_{m,\ell}) \\
%&+ \frac{1}{M_{L}} \sum_{m=1}^{M_L} [g(X(T)) - g_{\ell-1} \indicator{A_{L-1}}](\omega_{m,L})  \PERIOD
%\end{align*}
%Therefore, it is possible to add another stopping criteria to Algorithm \ref{alg:Cal} related to the comparison between the estimated cost of two hybrid coupled paths and the cost of hybrid-exact coupled paths. Please note that the condition $\delta_{L} \avg{N_{\TL,L}}{\cdot}\avg{g_{L}}{\cdot} {\leq} TOL^2$ trivially holds because $\avg{N_{\TL,L}}{\cdot}$ is zero. In the two examples analyzed, there are no significative computational gains in the estimation phase by using that stopping rule and its corresponding estimator. This unbiased estimator is inspired in the work of Anderson and Higham, \cite{Anderson2012}.
%\end{rem}
%
%\subsection{Phase III}
%From Phase II, we found that, to compute our multilevel Monte Carlo estimator, $\mathcal{M}_L$, for a given tolerance, we have to run $M_{0}^*$ single hybrid paths with parameters $(\Delta t_0,\delta_0)$ and  $M_{\ell}^*$ coupled hybrid paths with parameters $(\Delta t_{\ell-1},\Delta t_\ell, \delta_\ell)$ for $\ell=1,2,\ldots,L^*$. But, we will follow a slightly different strategy: we run half of the required simulations and use them to update our estimations of  the sequences  
%$\seqof{\EE_{I,\ell}}{\ell=0}{L^*}$, $\seqof{\mV{\ell}}{\ell=0}{L^*}$, and $\seqof{\psi_\ell}{\ell=0}{L^*}$. Then, we solve problem \ref{eq:auxminprob} again and re-calculate the values of $M_\ell^*$, for all $\ell$. We proceed iteratively until convergence. 
%In this way, we take advantage of the information generated by new simulated paths and update the estimations of the sequences of weak errors, costs, and variances, obtaining more control over the total work of the method.
%
%\subsection{Bayesian Updating}
%To solve the family of optimization problems \eqref{eq:auxminprob}, it is required that we previously know the sequences  
%$\seqof{\EE_{I,\ell}}{\ell=0}{L}$, $\seqof{\mV{\ell}}{\ell=0}{L}$, and $\seqof{\psi_\ell}{\ell=0}{L}$.
%Notice that  $\mV{0}$ and $\psi_0$ are related to single paths, but for $\ell\geq 1$, 
%$\mV{\ell}$ and $\psi_\ell$ are related to coupled paths.
%We assume for now that we know, for every $\ell$, the value $\delta_\ell$ such that 
%$\EE_{E,\ell}{<}TOL^2$. 
%
%To estimate $\EE_{I,0}$, $\mV{0}$ and $\psi_0$, we just sample single-level paths with parameters $(h_0,\delta_0)$ and use a standard Monte Carlo estimator. The number of samples should be enough for controlling the coefficients of variation of the mentioned estimators.  In principle, we could follow the same strategy to estimate $\EE_{I,\ell}$, $\mV{\ell}$ and $\psi_\ell$, for $\ell{\geq}1$. But notice that for large values of $\ell$, not only is $\psi_\ell$ large, but a large number of sampled coupled hybrid paths are also needed to estimate those quantities with a small coefficient of variation. 
%To overcome this problem, we apply a Bayesian approach that uses the information obtained from sampling at levels $\ell{-}2$ and $\ell{-}1$ to establish prior distributions for $\EE_{I,\ell}$, $\mV{\ell}$ and $\psi_\ell$. In this way, just  by drawing few coupled paths at levels $\ell{-}1$ and $\ell$, we can update our knowledge about those three mentioned quantities.
%This estimation strategy turns out to be very efficient in terms of  work, and also reasonably accurate.
%This is mainly due to the regularity of those quantities as functions of the level $\ell$, for instance, see Figure \ref{fig:dec2-diag}. 
%
%Let us consider with some detail the case of $\mV{\ell}$. 
%There are two sources of prior information to estimate $\mV{\ell}$: i) extrapolating linearly the values 
%$\log{\hmV{\ell-1}}$ and $\log{\hmV{\ell-2}}$, and then computing 
%$\tilde{\mathcal{V}}_{\ell}:= \hmV{\ell-1}^2/\hmV{\ell-2}$, 
%ii) sampling pairs $\LP \bar X_{\ell-2}, \bar X_{\ell-1}\RP$, and then using $\bar X_{\ell-1}$ in the dual formulas \eqref{eq:Ephinu}, \eqref{eq:Vphinu} and \eqref{eq:varhat} to obtain $\hat{\mathcal{V}}_{\ell,\delta_{\ell-1}}$. Of course, this last estimate corresponds to $\delta_{\ell-1}$ and not to $\delta_{\ell}$, but, still, it is useful prior information. 
%
%In our calculations, we use a lognormal model for the prior distribution of $\mV{\ell}$, where its parameters are estimated by combining  $\tilde{\mathcal{V}}_{\ell}$ and $\hat{\mathcal{V}}_{\ell,\delta_{\ell-1}}$.
%Using conjugate analysis (\cite{Fink97, Gelman}), we update this prior knowledge by sampling coupled paths $\LP \bar X_{\ell-1}, \bar X_{\ell}\RP$ and use $\bar X_{\ell-1}$ to compute dual weighted estimates of $\mV{\ell}$.
%Note that we do not use the observations $\Delta g_{\ell}:= \LP g_{\ell}-g_{\ell-1}\RP \indicator{A_{\ell{-}1}} \indicator{A_{\ell}}$ due to its high variance. 
%Figure \ref{fig:arrows} depicts the Bayesian updating process for the estimation of ${\mathcal{V}}_{\ell}$.
%\begin{figure}[h!]
%\begin{displaymath}
%    \xymatrix{ (\bar X_{\ell-2}, \blue{\bar X_{\ell-1}})_{\delta_{\ell-1}}  \ar@{~>}[dr] \ar[dd] &  &\\
%                   & \blue{\hat{\mathcal{V}}_{\ell,\delta_{\ell-1}}} \&\,\, \tilde{\mathcal{V}}_{\ell}\,\,  (\text{prior}) \ar@{.>}[d] &\\
%              (\blue{\bar X_{\ell-1}}, \bar X_{\ell})_{\delta_{\ell}}\ar@{~>}[r]  &  \hat{\mathcal{V}}_{\ell,\delta_{\ell}} }
%\end{displaymath}
%\caption{Scheme for the Bayesian updating of $\hat{\mathcal{V}}_{\ell,\delta}$.}
%\label{fig:arrows}
%\end{figure}
%We proceed analogously with the quantities  $\EE_{I,\ell}$, and $\psi_\ell$ 
%\clearpage
\section{A Control Variate Based on a Deterministic Time Change}
\label{sec:cvar}
In this section, we motivate a novel control variate for the random variable $X(T,\omega)$ defined by the random time change representation,
\begin{equation*}
X(T,\omega) = x_0 + \sum_j \nu_j Y_j\LP \int_0^T a_j(X(s))\,ds,\omega\RP\PERIOD
\end{equation*}

First, we replace the independent Poisson processes, $\LP Y_j(s,\omega) \RP_{s\geq 0}$, by the identity function. This defines the deterministic mean field, 
\begin{equation*}
Z(T) = x_0 + \sum_j \nu_j  \int_0^T a_j(Z(s))\,ds\PERIOD
\end{equation*}

Next,  we consider the random variable 
\begin{equation*}
\tilde X(T,\omega) = x_0 + \sum_j \nu_j Y_j\LP \int_0^T a_j(Z(s))\,ds, \omega\RP\COMMA
\end{equation*}
which uses the same realizations of $\LP Y_j(s,\omega) \RP_{s\geq 0}$ that define $X(T,\omega)$. 
In this way, we expect some correlation between  $X(T)$ and $\tilde X(T)$. 
Since $\expt{\tilde X(T)} = Z(T)$ is a computable quantity, we have that $\tilde X(T)$ is a potential control variate for $X(T)$ obtained at almost negligible extra computational cost.

We have that $\tilde X(T,\omega)$ can be considered as a deterministic time change approximation of $X(T,\omega)$.

To implement this idea, we first consider the sequence $Z_{k}$, defined as a forward Euler discretization of the mean field over a suitable mesh,
 $\{t_0{=}0,  t_1,\ldots, t_K{=}T\}$, $\Delta t_k :=  t_{k+1}{-} t_k$, $k{=}0,1,\ldots,K{-}1$; that is,
$$\left\{ 
\begin{array}{ll}
Z_{k+1} = Z_{k} + \sum_j \nu_j a_j(Z_{k}) \Delta t_{k}, & k{=}0,\ldots,K{-}1\\
Z_0 = x_0  &
\end{array}
\right. .
$$
The sequence  $Z_{k}$ allow us to define another sequence, $\hat \Lambda_{j,k}$, by 
$$\left\{ 
\begin{array}{ll}
\hat \Lambda_{j,k+1} = \hat \Lambda_{j,k} + a_j(Z_k)\Delta t_{k}, & k{=}1,\ldots,K{-}1\\
\hat \Lambda_{j,0} = 0  &
\end{array}
\right. ,
$$
where $\hat \Lambda_{j,K}$ approximates $\int_0^T a_j(Z(s))\,ds$.



Then, for each realization of $\bar{X}(T,\omega)$, which is an approximation of ${X}(T,\omega)$, we compute the control variate:
\begin{equation}\label{eq:XK}
\hat X_K = x_0 + \sum_j \nu_j Y_j\LP \hat \Lambda_{j,K}\RP\COMMA
\end{equation}
which is the corresponding approximation of  $\tilde{X}(T,\omega)$ and  has the computable expectation 
$$\mu_K := \expt{\hat X_K} = x_0 + \sum_j \nu_j  \hat \Lambda_{j,K}. $$


Now, we consider the random sequence, $\{\bar X_n(\omega)\}_{n=0}^{N(\omega)}$, generated in this case by the mixed algorithm. Here, $\bar X(\omega)_{N(\omega)}$ is an approximation of $X(T,\omega)$.
The sequence of mixed random times, $\{ \bar \Lambda_{j,n}(\omega) \}$, is defined by
$$\left\{ 
\begin{array}{ll}
\bar \Lambda_{j,n+1} = \bar \Lambda_{j,n} + a_j(\bar X_n(\omega))\Delta s_{n}, & n{=}0,\ldots,N(\omega){-}1\\
\bar \Lambda_{j,0} = 0 &
\end{array}
\right. \COMMA
$$
over the mesh $\{s_0{=}0,  s_1,\ldots, s_{N(\omega)}{=}T\}$, $\Delta s_n :=  s_{n+1}{-} s_n$, $n{=}0,1,\ldots,N(\omega){-}1$.

%$$Z(T) = x_0 + \sum_j \nu_j  \underbrace{\int_0^T a_j(Z(s))\,ds}_{\hat \Lambda_{j,T}} \COMMA$$

At this point, it is crucial to observe that we can keep track of the values $Y_j(\bar \Lambda_{j,n},\omega)$, since at each step of the approximation algorithm, we are sampling the increments of the processes, $Y_j$. From now on, we omit $\omega$ in our notation.

The values $Y_j\LP \hat \Lambda_{j,K}\RP$, required in \eqref{eq:XK}, can by obtained by sampling the process $Y_j$ as follows.
For each realization of $\bar{X}$, we have two scenarios:
\begin{enumerate}
\item for some $n$, $\bar \Lambda_{j,n} <\hat \Lambda_{j,K}< \bar \Lambda_{j,n+1}$.  Since $Y_j\LP \bar \Lambda_{j,n} \RP$ and $Y_j\LP \bar \Lambda_{j,n+1} \RP$ are known, we sample a Poissonian bridge (binomial), \ie,
$$Y_j(\hat \Lambda_{j,K})\sim Y_j\LP \bar \Lambda_{j,n} \RP + \text{binomial}\LP Y_j\LP \bar \Lambda_{j,n+1} \RP-Y_j\LP \bar \Lambda_{j,n} \RP, \frac{\hat \Lambda_{j,K} - \bar \Lambda_{j,n}}{\bar \Lambda_{j,n+1} - \bar \Lambda_{j,n}}\RP. $$ 
\item $\bar \Lambda_{j,K} >\hat \Lambda_{j,N}$. Since we know the value $Y_j\LP \bar \Lambda_{j,N} \RP$,
we just have to sample a Poisson random variate as follows:
$$Y_j(\hat \Lambda_{j,K})\sim Y_j\LP \bar \Lambda_{j,N} \RP + \text{Poisson}(a_j(\bar X_N)(\hat \Lambda_{j,K}- \bar \Lambda_{j,N})).$$
\end{enumerate}


%\begin{itemize}
%\item 
%\item $\bar \Lambda_{j,n+1} = \bar \Lambda_{j,n} + a_j(\bar X_n)\Delta s_{n}$, $\bar \Lambda_{j,0} = 0$  $j{=}1,...,J$
%\item over a mesh: $\{s_0{=}0,  s_1,\ldots, s_N{=}T\}$, $\Delta s_n :=  s_{n+1}{-} s_n$, $n{=}0,1,\ldots,N{-}1$
%\end{itemize}

Finally, using the aforementioned control variate, we can estimate $\expt{g(\bar{X}(T))}$ with 
$$\frac{1}{M} \sum_{m=1}^{M} g(\bar{X}_N(\omega_m)) - \beta \frac{1}{M}\sum_{m=1}^{M} (g(\hat{X}_K(\omega_m))-g(\mu_K)) \COMMA$$
for any linear functional, $g$. 
For polynomial observables, $g$, this estimator can be easily extended by Taylor expansions in such way that we can estimate $\expt{g(\tilde X(T))}$ by powers, $g\LP \LP\expt{\tilde X(T)}\RP^k\RP$.
%construct a control process
%$$\hat X(T) = x_0 + \sum_j \nu_j Y_j\LP \int_0^T a_j(Z(s))\,ds\RP \PERIOD$$

\begin{rem}[Reducing the variance at the coarsest level]
The main application of the deterministic time change control variate, $\tilde X(T)$, in this work is at the coarsest level of our multilevel hierarchy. 
Consider the trivial decomposition 
$$g(\bar X_0(T)) = g(\tilde X(T)) + \LP g(\bar X_0(T))-  g(\tilde X(T)) \RP.$$ Therefore, 
$$\expt{g(\bar X_0(T))} = \expt{g(\tilde X(T))} + \expt{g(\bar X_0(T))-  g(\tilde X(T))}.$$ 
Since we can compute exactly $ \expt{g(\tilde X(T))}$, we just have to estimate $\expt{g(\bar X_0(T))- g(\tilde X(T))}$ instead of $\expt{g(\bar X_0(T))}$ in our multilevel scheme. The computational gain lies in the fact that $\var{g(\bar X_0(T))- g(\tilde X(T))}$ could be substantially lower than $\var{g(\bar X_0(T))}$.
\end{rem}
\begin{rem}[Computational Cost]
An advantage of this control variate is that the computational cost is almost negligible because we only need to store two scalars, $\bar \Lambda_{j,n}$ and $\bar \Lambda_{j,n+1}$, for each reaction, $j$. These values are determined at each step by $a_j(\bar X_n)$, which is a quantity that is already computed at each time step of the mixed algorithm. Also, for each realization of the control variate, at most one Poisson random variate is needed for each reaction channel.
\end{rem}
\begin{rem}[Empirical Time Change]
We can also compute the final times, $\hat \Lambda_{j,K}$, using a sample average of mixed paths instead of the mean field. We found no significant improvements when using that approach, which requires a lot more computational work. We conjecture that, for settings in which the mean field is not representative, this approach is the only reasonable option.
\end{rem}

\section{Numerical Examples}
\label{sec:examples}

%\newcommand{\MLVirus}{MLM_Virus_2014_4_27_12_22_29}
\newcommand{\MLVirus}{MLM_Virus_2014_4_28_10_14_33}
\newcommand{\Stiff}{MLM_Simple_stiff_2014_5_25_11_7_22}

In this section, we present two examples to illustrate the performance of our proposed method, and we compare the results with the hybrid MLMC approach given in \cite{ourML}. For benchmarking purposes, we use Gillespie's Stochastic Simulation Algorithm (SSA) instead of the Modificed Next Reaction Method (MNRM) because the former is widely used in the literature.

%\subsection*{Hybrid vs. Mixed}

\subsection*{Intracellular Virus Kinetics}

%The classical radioactive decay model provides a simple and important example for the application of our  method. 
This model, first developed in \cite{srivastava2002stochastic}, has four species and six reactions,
\begin{itemize}
\item $E \xrightarrow{1} E {+} G$, the viral template (E) forms a viral genome (G),
\item $G \xrightarrow{0.025} E$, the genome generates a new template,
\item $E \xrightarrow{1000} E {+} S$, a viral structural protein (S) is generated,
\item $G{+}S \xrightarrow{7.5\times10^{-6}} V$, the virus (V) is produced, 
\item $E\xrightarrow{0.25}\emptyset$, $S\xrightarrow{2}\emptyset$ degradation reactions.
\end{itemize}
Its stoichiometric matrix and its propensity functions, $a_j:\zset_+ \rightarrow \rset$, are given by 
\begin{align*}
\nu =\left( 
 \begin{array}{ccccc}     1    & 0 &    0 & 0\\
     -1    & 0 &    1 & 0 \\
     0   & 1 &    0 & 0 \\
     -1 & -1 &0 &1 \\
    0    & 0 &   -1&  0 \\
     0   & -1 &  0 & 0 
 \end{array} 
 \right) ^{tr}\mbox{    and     }  a(X) =\left( \begin{array}{c} E \\ 0.025 \, G \\ 1000 \, E \\ 7.5{\times} 10^{-6} G \, S \\ 0.25 \, E \\ 2 \, S \end{array} \right)\COMMA
\end{align*}
respectively. 

In this model, $X(t)=(G(t),S(t),E(t),V(t))$, and $g(X(t))=V(t)$, the number of viruses produced. The initial condition is $X_0{=}(0, 0, 10, 0)$ and the final time is $T{=}20$.
%This particular setting can be found in \cite{Anderson2012}.
%\begin{align*}
%X \xrightarrow{c} \emptyset\PERIOD
%\end{align*}
%Its stoichiometric matrix, $\nu \in \rset$, and the propensity function, $a:\zset_+ \rightarrow \rset$, are given by 
%\begin{align*}
%\nu = -1 \mbox{   and   }  a(X) = c\,X \PERIOD
%\end{align*}
%Here, we choose $c = 1$ and $g(x)=x$. In this particularly simple example, we have that $\expt{g(X(T))|X(t) = X_0} = X_0 \exp(-c(T{-}t))$.
%Consider the initial condition $X_0{=}10^5$, and the final time, $T{=}0.5$. In this case, the process starts relatively far from the boundary, \ie, it is a tau-leap dominated setting. 
This example is interesting because i) it shows a clear separation of time scales, ii) our previous hybrid Chernoff method has no compuational work gain with respect to an exact method, and iii) in \cite{Anderson2012} the authors take an alternative approach, not using the multilevel aspect of their paper.

We now analyze an ensemble of 10 independent runs of the phase II algorithm (see Section \ref{sec:estim}), using different relative tolerances. % (its summary is given in Table \ref{tab:vir}). 
In Figure \ref{fig:dec2-worktime}, we show the total predicted work (runtime) for the multilevel mixed method and for the SSA method, versus the estimated error bound. We also show the estimated asymptotic work of the multilevel mixed method. We remark that the computational work of the multilevel \emph{hybrid} method is the same as the work of the SSA. 
%To predict the work of the single-level method, we used the corresponding calibration algorithm for the single-level, given in \cite{ourSL}. 
%In the right panel, we show, for different tolerances, the actual work (runtime), using a 12 core Intel GLNXA64 architecture and MATLAB version R2012b. 
%In Table \red{PUTREF}, we summarize an ensemble run of the calibration algorithm, where
%$\hat W_{\ML}:=\sum_{\ell=0}^L M_\ell \hat{\psi}_{\ell}$ is the estimated computational work of the multilevel estimator and $W_{\ssa}:=M_{\ssa}\, \,C_1\,\avg{\NSSAP}{\cdot}$ is the corresponding estimated computational work of the SSA. Both quantities can be computed using Algorithm \ref{alg:Cal}. 
%$W_{\ML}$ is the average actual computational work of the multilevel estimator (the sum of all the seconds taken to compute the estimation) and $W_{\ssa}$ is the corresponding average actual work of the SSA. We compare those values with the corresponding estimations, $\hat W_{\ML}$, and $\hat W_{\ssa}$. The difference between the estimated work and the actual one gets smaller as we refine the tolerance. %For this setting, our method is asymptotically more than ten times faster than the SSA.

\begin{figure}[h!]
\centering
\begin{minipage}{0.49\textwidth}
\includegraphics[scale=0.33]{Figures/\MLVirus/Virus_T20_work_err_ens}
\end{minipage}
\hfill
\begin{minipage}{0.49\textwidth}
%\includegraphics[scale=0.31]{Figures/\MLVirus/Virus_T20_comp_time_err_ens}
\begin{tabular}{lcc}
\label{tab:vir}
$TOL$ & $L^*$ & $\frac{\hat{W}_{ML}}{\hat{W}_{\ssa}}$ \\ \noalign{\smallskip} \hline\noalign{\smallskip} 
1.00e-01 & 1.0  & 0.02 $\pm$0.001   \\ 
5.00e-02 & 1.0  & 0.02 $\pm$0.001   \\ 
2.50e-02 & 1.2 $\pm$0.261 & 0.02 $\pm$0.001 \\ 
1.25e-02 & 2.2 $\pm$0.261 & 0.03 $\pm$0.002   \\ 
6.25e-03 & 3.4 $\pm$0.320 & 0.04 $\pm$0.004  \\ 
3.13e-03 & 4.6 $\pm$0.320 & 0.04 $\pm$0.002   \\ 
1.56e-03 & 5.8 $\pm$0.261 & 0.06 $\pm$0.008  \\ 
7.81e-04 & 7.4 $\pm$0.433 & 0.07 $\pm$0.006  \\ 
3.91e-04 & 8.6 $\pm$0.320 & 0.06 $\pm$0.007   \\ 
\noalign{\smallskip}\hline 
\end{tabular}
\end{minipage}
\caption{Left: Predicted work (runtime) versus the estimated error bound, with $95\%$ confidence intervals. The multilevel mixed method is preferred over the SSA and the multilevel hybrid  method for all the tolerances. 
Right: Details for the ensemble run of the phase II algorithm.
Here, $\hat{W}_{ML}= \sum_{\ell =0}^{L^*}\hat \psi_\ell M_{\ell}$ and  $\hat{W}_{\ssa} = M_{\ssa}\, \,C_\ssa \,\avg{\NSSAP}{\cdot}$. 
As an example, the fourth row of the table tells us that, for a tolerance $TOL{=}1.25\cdot 10^{-2}$, 2.2 levels are needed on average. 
The work of the multilevel hybrid method is, on average, $3\%$ of the work of the SSA and the multilevel hybrid method. Confidence intervals at 95\% are also provided.
%Right: Actual computational work (runtime) versus the estimated error bound. Notice that the computational complexity has order $\Ordo{TOL^{-2}}$.
}
\label{fig:dec2-worktime}
\end{figure}
%\begin{table}[h!]
%\centering
%\begin{tabular}{lclllcl}
%\input{Figures/\MLVirus/Virus_T20_data1}
%\end{tabular}
%
%\bigskip
%\caption{Details for the ensemble run of phase II.
%%Here $C_{\ML}:= \sum_{\ell=0}^L rt_{C,\ell} M_\ell$ and  $C_{\ssa} := M_{\ssa}\, \,C_1\,\avg{\NSSAP}{\cdot}$.
%As an example, the fourth row of the table tells us that, for a tolerance $TOL{=}1.25\cdot 10^{-2}$, 2.4 levels are needed on average. 
%The work of the multilevel hybrid method is, on average, $3\%$ of the work of the SSA and the multilevel hybrid method. Confidence intervals at 95\% are also provided.
%}
%\label{tab:vir}
%\end{table}

%\begin{figure}[h!]
%\centering
%\begin{minipage}{0.49\textwidth}
%\includegraphics[scale=0.33]{Figures/\DECtwo/Virus_T20_time_est_err}
%\end{minipage}
%\caption{Predicted and actual runtime versus the estimated error bound for one run of the Algorithm \ref{alg:Cal}. Observe how close the predictions are to the actual work (runtime).}
%\label{fig:dec2-predactual}
%\end{figure}
In Figure \ref{fig:dec2-diag}, we can observe how the estimated weak error, $\WEH{\ell}$, and the estimated variance of the difference of the functional between two consecutive levels, $\hmV{\ell}$, decrease linearly as we refine the time mesh, which corresponds to a tau-leap dominated regime. This linear relationship for the variance starts at level $1$, as expected. When the MNRM dominated regime is reached, both quickly converge to zero as expected.
The estimated total path work, $\hat \psi_{\ell}$, increases as we refine the time mesh. Observe that it increases linearly for the coarser grids, until it reaches a plateau, which corresponds to the pure MNRM case where the computational cost is independent of the grid size. %This is because the work needed for generating Poisson random variables becomes smaller as we refine the time mesh.
In the lower right panel, we show the total computational work, only for the cases in which $\WEH{\ell} < TOL{-}TOL^2$. %The minimum is attained at level $L^*=13$.
\newcommand{\TOLDEC}{3.91e-04}
\begin{figure}[h!]
\centering
\begin{minipage}{0.49\textwidth}
\includegraphics[scale=0.31]{Figures/\MLVirus/Virus_T20_h_WE_\TOLDEC_ens}
\end{minipage}
%\hfill
\begin{minipage}{0.49\textwidth}
\includegraphics[scale=0.31]{Figures/\MLVirus/Virus_T20_h_var_\TOLDEC_ens}
\end{minipage}
\begin{minipage}{0.49\textwidth}
\includegraphics[scale=0.31]{Figures/\MLVirus/Virus_T20_h_psi_\TOLDEC_ens}
\end{minipage}
%\hfill
\begin{minipage}{0.49\textwidth}
\includegraphics[scale=0.31]{Figures/\MLVirus/Virus_T20_l_Work_\TOLDEC_ens}
\end{minipage}
\caption{Upper left: estimated weak error, $\WEH{\ell}$, as a function of the time mesh size, $h$. Upper right: estimated variance of the difference between two consecutive levels, $\hmV{\ell}$, as a function of $h$. Lower left: estimated path work, $\hat \psi_\ell$, as a function of $h$. Lower right: estimated total computational work, $\sum_{l=0}^L \hat \psi_l M_l$, as a function of the level, $L$.
}
\label{fig:dec2-diag}
\end{figure}

\begin{figure}[h!]
\centering
\begin{minipage}{0.49\textwidth}
\includegraphics[scale=0.31]{Figures/\MLVirus/Virus_T20_TOL_Stat_ens}
\end{minipage}
%\hfill
\begin{minipage}{0.49\textwidth}
\includegraphics[scale=0.31]{Figures/\MLVirus/Virus_T20_VarPsi_\TOLDEC_ens}
\end{minipage}
\caption{Left: Percentage of the statistical error over the total error. As we mentioned in Section \ref{sec:complex}, it is well above $0.5$ for all the tolerances. Right: $\sqrt{\hmV{\ell}\hat{\psi}_\ell}$, as a function of $\ell$, for the smallest tolerance, which  decreases as the level increases. Observe that the contribution of level 0 is less than 50\% of the sum of the other levels.
}
\label{fig:statdec2}
\end{figure}

In Figure \ref{fig:dec2-out}, we show the main outputs of the phase II algorithm, $\delta_\ell$ and $M_\ell$ for $\ell=0,...,L^*$, for the smallest considered tolerance. In this example, $L^*$ is 8 or 9, depending on the run. We observe that the number of realizations decreases slower than linearly, from levels $1$ to $L^*{-}1$, until it drops, due to the change to a MNRM dominated regime. % to $M_{L^*}{=}100$.

\begin{figure}[h!]
\centering
\begin{minipage}{0.49\textwidth}
\includegraphics[scale=0.31]{Figures/\MLVirus/Virus_T20_l_delta_\TOLDEC_ens}
\end{minipage}
%\hfill
\begin{minipage}{0.49\textwidth}
\includegraphics[scale=0.31]{Figures/\MLVirus/Virus_T20_l_M_\TOLDEC_ens}
\end{minipage}
\caption{The one-step exit probability bound, $\delta_\ell$, and $M_\ell$ for $\ell{=}0,1,...,L^*$, for the smallest tolerance.
%Upper left: Lower left. Upper right: Lower right:
}
\label{fig:dec2-out}
\end{figure}

In Figure \ref{fig:effvir}, 
%we show the performance of the formula \eqref{eq:varhatestimated}, implemented in Algorithm \ref{alg:varggl}, used to estimate the strong error, $\mV{\ell}$. The quotient of $\hmV{\ell}$ over a standard Monte Carlo estimate of $\mV{\ell}$ is almost 1 for the first ten levels. At levels 11 and 12, we obtain  0.99 and  0.91, respectively. Both quantities are estimated using a coefficient of variation less than 5\%, but there is a remarkable difference in terms of computational work in favor of our dual-weighted estimator. In the same figure, we also 
we show $TOL$ versus the actual computational error. It can be seen that the prescribed tolerance is achieved with the required confidence of 95\%, since $C_A{=}1.96$, for all the tolerances.
The QQ-plot in the right part of Figure \ref{fig:effvir} was obtained as follows: i) for the range of tolerances specified in the first column of Table \ref{tab:vir}, we ran the phase II algorithm 5 times, ii) for each output of the calibration algorithm, we sampled the multilevel estimator $\mathcal{M}_L$, defined in \ref{MLMCest}, 100 times. This plot reaffirms our assumption about the Gaussian distribution of the statistical error. %Observe that the estimates are concentrated around the theoretical value $X_0 \exp(-c(T{-}t)) = 10^5 \exp(-0.5)\approx 6.0653e+04$.
\begin{figure}[h!]
\centering
%\begin{minipage}{0.49\textwidth}
%\includegraphics[scale=0.31]{Figures/Virus_T20_strong_errors}
%\end{minipage}
%\hfill
\begin{minipage}{0.49\textwidth}
\includegraphics[scale=0.31]{Figures/\MLVirus/Virus_T20_TOL_total_error}
\end{minipage}
\begin{minipage}{0.49\textwidth}
\includegraphics[scale=0.31]{Figures/\MLVirus/Virus_T20_qqplop}
\end{minipage}
\caption{Left: $TOL$ versus the actual computational error. The numbers above the straight line show the percentage of runs that had errors larger than the required tolerance. We observe that in all cases, the computational error follows the imposed tolerance with the expected confidence of 95\%. Right: quantile-quantile plot based on realizations of $\mathcal{M}_L$.
}
\label{fig:effvir}
\end{figure}

\begin{rem}
In the simulations, we observe that, as we refine $TOL$, the optimal number of levels approximately increases logarithmically, which is a desirable feature. We fit the model $L^*=a\log(TOL^{-1})+b$, obtaining $a{=}1.47$ and $b{=}3.56$.
\end{rem}

\begin{rem}[Pareto rule]
\label{rem:pareto}
Using the cost-based rule (see remark \ref{rem:pare}), we estimate the threshold for the Pareto rule, obtaining $\nu = 0.95419$. It turns out that, for this example, $\hat W_{MixPareto}/ \hat W_{Mix}$ ranges from $0.6$ to $0.75$ (for most $TOL$s). This shows that it is possible to  increase the computational work gains further in some examples. %in which the threshold does not vary along the paths.
\end{rem}

\begin{rem}
The savings in computational work when generating Poisson random variables heavily depend on MATLAB's performance capabilities. In fact, we would expect better results from our method if we were to implement our algorithms in more performance-oriented languages or if we were to sample Poisson random variables in batches.
\end{rem}

\subsection*{A Simple Stiff System}
This model, adapted from \cite{gillespie2005}, has three species and a mixture of fast and slow reaction channels,
$$X_1 \xrightleftharpoons[c_2]{c_1} X_2 \xrightarrow{c_3} X_3 \xrightarrow{c_4} \emptyset, \,\,\, c_2 \gg c_3 > c_4 \PERIOD$$

Its stoichiometric matrix and propensity functions, $a_j:\zset_+ \rightarrow \rset$, are given by 
\begin{align*}
\nu =\left( 
 \begin{array}{ccc}     -1  &   1 &    0 \\
     1  & -1 &   0\\
     0 &   -1  &   1\\
     0  &   0  &  -1
 \end{array} 
 \right) ^{tr}\mbox{    and     }  a(X) =\left( \begin{array}{c} c_1 X_1 \\ 
 c_2 X_2 \\ c_3 X_2 \\ c_4 X_3 \end{array} \right)\COMMA
\end{align*}
respectively, where $g(X(t))=X_3(t)$.
In this model, successive firings of the reaction $X_2 \rightarrow X_3$ are separated by many reversible firings between $X_1$ and $X_2$, which takes a lot of computational work in a standard SSA run. In \cite{gillespie2013}, Gillespie et al. claim that this inefficiency cannot be addressed using ordinary tau-leaping because of the stiffness of the system. We show here that we have substantial gains using our mixed method, which also controls the global error. In this example, we also show the performance of the control variate idea, presented in Section \ref{sec:cvar}.
We analyze 10 independent runs of the phase II algorithm (see Section \ref{sec:estim}), using different relative tolerances. % (its summary is given in Table \ref{tab:vir}). 
In Figure \ref{fig:stiff-worktime}, we show the total predicted work (runtime) for the multilevel mixed method with and without a control variate at level 0 and for the SSA method versus the estimated error bound. We also show the estimated asymptotic work of the multilevel mixed method. % For practical tolerances, gains are of the order of two magnitudes. 
Observe that, for practical tolerances, the computational work gains with respect to the SSA method, when using the control variate, are of a factor of 500 times. Without using the control variate, computational gains are also substantial. %Also the estimated asymptotic work when using the control variate is half when not using it. %It is useful to mention that, for this particular quantity of interest, which is an integer value, asymptotic tolerances has little use in practice.

\begin{figure}[h!]
\centering
\begin{minipage}{0.49\textwidth}
\includegraphics[scale=0.34]{Figures/\Stiff/Simple_stiff_T100_comparison_cv0}
\end{minipage}
\hfill
\begin{minipage}{0.49\textwidth}
%\includegraphics[scale=0.31]{Figures/\MLVirus/Virus_T20_comp_time_err_ens}
\begin{small}
\begin{tabular}{lccc}
$TOL$ & $L^*$ & $\frac{\hat{W}_{ML_{cv}}}{\hat{W}_{\ssa}}$ & $\frac{\hat{W}_{ML}}{\hat{W}_{\ssa}}$ \\ \noalign{\smallskip} \hline\noalign{\smallskip} 
3.13e-03 & 1.0  & 0.002 $\pm$0.0004 & 0.03 $\pm$0.001\\ 
1.56e-03 & 1.0  & 0.003 $\pm$0.0004 & 0.04 $\pm$0.001 \\ 
7.81e-04 & 1.0  & 0.003 $\pm$0.0010 & 0.04 $\pm$0.002\\ 
3.91e-04 & 1.0  & 0.004 $\pm$0.0004 & 0.06 $\pm$0.003\\ 
1.95e-04 & 2.0  & 0.013 $\pm$0.0015 & 0.09 $\pm$0.008\\ 
9.77e-05 & 3.0  & 0.027 $\pm$0.0040 & 0.13 $\pm$0.016\\ 
4.88e-05 & 4.0  & 0.065 $\pm$0.0146 & 0.19 $\pm$0.025 \\ 
2.44e-05 & 6.0  & 0.100 $\pm$0.0136 & 0.21 $\pm$0.020\\ 
1.22e-05 & 6.0  & 0.109 $\pm$0.0299 & 0.22 $\pm$0.029\\ 
6.10e-06 & 6.0  & 0.108 $\pm$0.0168 & 0.19 $\pm$0.020\\ 
\noalign{\smallskip}\hline 
\end{tabular}
\end{small}
\end{minipage}
\caption{Left: Predicted work (runtime) versus the estimated error bound, with $95\%$ confidence intervals, for the simple stiff model with and without using the control variate at level 0, as described in Section \ref{sec:cvar}. 
Right: Details of the ensemble run of the phase II algorithm using the control variate (third column) and without using the control variate (fourth column).
As an example, the fifth row of the table tells us that, for a tolerance $TOL{=}1.95\cdot 10^{-4}$, 2 levels are needed on average. 
The work of the multilevel mixed method using the control variate at level 0 is, on average, $1\%$ of the work of the SSA. When not using the control variate, it is $9\%$.
Confidence intervals at 95\% are also provided.
%Right: Actual computational work (runtime) versus the estimated error bound. Notice that the computational complexity has order $\Ordo{TOL^{-2}}$.
}
\label{fig:stiff-worktime}
\end{figure}

%\green{
%\begin{enumerate}
%\item{Prpoensities, jumps, observable, final time}
%\item{Why we put each Figure, what is important here}
%\item{Particular conclusions, comparison with Single Level}
%\end{enumerate}
%}
%\input{num_examples}

\section{Conclusions}
\label{sec:conclusions}

In this work, we addressed the problem of approximating the quantity of interest $\expt{g(X(T))}$, where $X$ is a non-homogeneous Poisson process that  describes a stochastic reactions network, and $g$ is a given suitable observable of $X$, within a given prescribed relative tolerance, $TOL{>}0$, up to a given confidence level at near-optimal computational work.

We developed an automatic, adaptive reaction-splitting multilevel Monte Carlo method, based on our Chernoff tau-leap mehthod \cite{ourSL,ourML}. 
Its computational complexity is $\Ordo{TOL^{-2}}$. This method can be therefore seen as a variance reduction of the SSA, which has the same complexity. %It has important advantages over a pure tau-leap method in a multilevel context specially for stiff problems. 
In our numerical examples, we obtained substantial gains with respect to SSA and, for systems in which the set of reaction channels can be adaptively partitioned into ``high'' and ``low'' activity, over our previous multilevel hybrid Chernoff tau-leap method.

We also presented a novel control variate for $g(X(T))$, which adds negligible computational cost when simulating a path of $X(T)$, and it may lead to additional dramatic cost reductions.

%Our next step is to extend this method to implicit tau-leap schemes and to introduce a spatial component to the problem setting.

%For reaching this optimality, we derived novel formulas based on dual-weighted residual estimations for computing the variance of the difference of the observables between two consecutive levels in coupled hybrid paths.
%Unlike the multilevel algorithms proposed in \cite{Anderson2012}, we do not need to distinguish between biased and unbiased methods. When our hybrid algorithm chooses exact paths at the bottom level, we automatically have an unbiased algorithm.
%Future extensions may involve better hybridization techniques as well as implicit and higher order versions of the hybrid MLMC.
%We developed a novel multilevel Monte Carlo approach for the hybrid Chernoff tau-leap algorithm presented in \cite{ourSL}. This approach has two main features. First, we developed a multilevel Monte Carlo estimator that has a high probability of providing estimates that lies within a certain prescribed tolerance TOL from the true value. Second, it provides all the parameter values needed for the simulation stage, \ie, number of levels,
%one-step exit probability bound for each level, initial time mesh and the number of simulations at each level, such that approximately optimizes the computational work.
%
%The bias of our estimator is $\expt{g(X(T))} - \LP \expt{g_{L} \indicator{A_{L}}} + \Ordo{L\,TOL^{\orden}} \RP$. 
%In principle, one could reduce the bias to $\Ordo{L\,TOL^{\orden}}$, following the ideas of Anderson and Higham in $\cite{Anderson2012}$, \ie, using at the last level, $L$, a coupling scheme between hybrid and exact paths. We found that, in the cases where the jump activity is high, coupling with exact paths could be unaffordable. 
%
%In that sense, for the cases in which we can afford the exact coupling, we also have a virtually unbiased method.
%
%Our approach for the optimization of the total computational work needs to compute the total cost of the coupled path, the variance of the difference between the quantity of interest computed with a mesh of level $\ell$ and of level $\ell{-}1$, and the discretization error for an extra level. %This computations depends on simulations done at a level one more that the optimal. 
%We may avoid this extra level using appropriate  models for the mentioned quantities, but we left this for future work.
%
%\red{TODO}
%
%In this work, we addressed the problem of accurately estimating the expected value of an observable of a Markov pure jump process at a given final time within a certain prescribed tolerance with high probability. 
%
%In this work, we proposed a hybrid algorithm that, at each step, adaptively chooses to adopt the SSA method when the work of the tau-leap step becomes high. As a consequence, the expected work of a hybrid path remains bounded by the expected work of an SSA path and potentially can be much smaller.
%
%The global exit error is related to the fact that, at any time, a tau-leap path can attain a non-physical value. %This is true because the increments of the tau-leap method are linear combination of the stoichiometric vectors with coefficients given  by independent Poisson random variables whose rates are proportional to the time step size. 
%Pre-leap checks are common techniques for dealing with this problem by controlling the time step size. Here, we presented a novel non-asymptotic Chernoff-type hard bound to control large deviations of linear combinations of independent Poisson random variables. This bound allows us not only to obtain a pre-leap check for the tau-leap method, which does not change the distribution of the increments nor does it require any type of assumption regarding the reactions that can occur, but also to estimate and control the global exit error. To the best of our knowledge, there is no previous attempt in the literature to estimate and control this type of error at the path level.
%
%Another important contribution of this work is a calibration algorithm that can determine if it is suitable to use the hybrid algorithm for a given problem and also that can provide the associated simulation parameters. In the hybrid case, the calibration algorithm provides the one-step exit probability bound, the time mesh and the number of hybrid paths that are needed for computing the mentioned expected value with low computational work.
%
%It is worth mentioning that, by simulating hybrid paths, we obtained accurate estimates of the average number of steps required by the SSA method to reach the final time. This is especially relevant in problems where the process visits regions of the state space where the total propensity is very high. 
%
%The numerical results that we obtain from different models show that the hybrid method proposed here is suitable for addressing problems in which one or more species has few individuals while the total propensity is high. In these types of problems (e.g., the gene transcription and translation model), the reaction-rate ODEs not only do not provide accurate approximation for the average behavior of the process, but also the cost of the exact methods is also high.   
%%The hybrid decision rule is based on a fine-tune modelling of the computational times associated to the exact and tau-leap methods. 
%Moreover, we observed that generating Poisson random variables makes the computational work of a tau-leap step much higher than the work of an SSA step. 
%This last argument, together with the advantages already discussed in terms of the time discretization error and the global exit error, adds more evidence in favor of  avoiding the tau-leap whenever possible.
%
%Our next step is to extend this hybrid algorithm to the Multi Level Monte Carlo setting \cite{GilesMLMC, Anderson2012}. We aim to obtain substantial  computational work gains with respect to the traditional exact methods (SSA or MNRM) and the single-level hybrid Chernoff tau-leap, and show that the computational complexity of this multilevel extension is of order $\Ordo{TOL^{-2}}$.


%\green{
%\begin{enumerate}
%\item{What we did.}
%\item{Why matters.}
%\item{We are better than ... because....}
%\end{enumerate}}
%\input{conclusions}

\section*{Acknowledgments}
The research reported here was supported by King Abdullah University of Science and Technology (KAUST).
The authors are members of the SRI Center for
Uncertainty Quantification in Computational Science and Engineering at KAUST. 
\appendix
%\input{CalibrationAlgorithms}
%\newpage
%*flatex input: [main_mixed.bbl]
\begin{thebibliography}{10}

\bibitem{poivariate}
J.~Ahrens and U.~Dieter.
\newblock Computer methods for sampling from gamma, beta, {P}oisson and
  bionomial distributions.
\newblock {\em Computing}, 12:223--246, 1974.

\bibitem{Anderson2012}
D.~{Anderson} and D.~{Higham}.
\newblock Multilevel {M}onte {C}arlo for continuous {M}arkov chains, with
  applications in biochemical kinetics.
\newblock {\em SIAM Multiscal Model. Simul.}, 10(1), 2012.

\bibitem{anderson2007modified}
D.~F. Anderson.
\newblock A modified next reaction method for simulating chemical systems with
  time dependent propensities and delays.
\newblock {\em The Journal of Chemical Physics}, 127(21), 2007.

\bibitem{KMC}
C.~C. Battaile.
\newblock The kinetic monte carlo method: Foundation, implementation, and
  application.
\newblock {\em Computer Methods in Applied Mechanics and Engineering}, Volume
  197(41-42):3386–3398, 2008.

\bibitem{gillespie2005}
Y.~Cao, D.~T. Gillespie, and L.~R. Petzold.
\newblock The slow-scale stochastic simulation algorithm.
\newblock {\em The Journal of Chemical Physics}, 122(1):--, 2005.

\bibitem{Cao2006}
Y.~Cao, D.~T. Gillespie, and L.~R. Petzold.
\newblock {Efficient step size selection for the tau-leaping simulation
  method.}
\newblock {\em The Journal of Chemical Physics}, 124(4):044109, 2006.

\bibitem{CaoPetzold2006}
Y.~Cao and L.~Petzold.
\newblock Accuracy limitations and the measurement of errors in the stochastic
  simulation of chemically reacting systems.
\newblock {\em Journal of Computational Physics}, 212(1):6--24, 2006.

\bibitem{kurtzmp}
S.~N. Ethier and T.~G. Kurtz.
\newblock {\em Markov Processes: Characterization and Convergence (Wiley Series
  in Probability and Statistics)}.
\newblock Wiley-Interscience, 2nd edition, 9 2005.

\bibitem{GibBruck}
M.~A. Gibson and J.~Bruck.
\newblock Efficient exact stochastic simulation of chemical systems with many
  species and many channels.
\newblock {\em The Journal of Physical Chemistry A}, 104(9):1876--1889, 2000.

\bibitem{gillespie_ssa}
D.~T. Gillespie.
\newblock A general method for numerically simulating the stochastic time
  evolution of coupled chemical reactions.
\newblock {\em Journal of Computational Physics}, 22:403--434, 1976.

\bibitem{gillespie_tau_leap}
D.~T. {Gillespie}.
\newblock {Approximate accelerated stochastic simulation of chemically reacting
  systems}.
\newblock {\em Journal of Chemical Physics}, 115:1716--1733, July 2001.

\bibitem{gillespie2013}
D.~T. Gillespie, A.~Hellander, and L.~R. Petzold.
\newblock Perspective: Stochastic algorithms for chemical kinetics.
\newblock {\em The Journal of Chemical Physics}, 138(17):--, 2013.

\bibitem{Harris2006}
L.~Harris and P.~Clancy.
\newblock A “partitioned leaping” approach for multiscale modeling of
  chemical reaction dynamics.
\newblock {\em J. Chem. Phys}, Volume 125, 2006.

\bibitem{Haseltine2002}
E.~Haseltine and J.~Rawlings.
\newblock Approximate simulation of coupled fast and slow reactions for
  stochastic chemical kinetics.
\newblock {\em J. Chem. Phys}, 117(15), 2002.

\bibitem{li}
T.~Li.
\newblock Analysis of explicit tau-leaping schemes for simulating chemically
  reacting systems.
\newblock {\em Multiscale Model. Simul.}, 6(2):417--436 (electronic), 2007.

\bibitem{ourSL}
A.~Moraes, R.~Tempone, and P.~Vilanova.
\newblock Hybrid {C}hernoff tau-leap.
\newblock {\em To appear in SIAM Multiscale Modeling and Simulation}, 2014.

\bibitem{ourML}
A.~Moraes, R.~Tempone, and P.~Vilanova.
\newblock Multilevel hybrid {C}hernoff tau-leap.
\newblock {\em arXiv:1403.2943}, 2014.

\bibitem{Plyasunov2005}
S.~Plyasunov.
\newblock Averaging methods for stochastic dynamics of complex reaction
  networks: description of multi-scale couplings.
\newblock {\em arXiv:physics/0510054v1}, 2005.

\bibitem{Puchalka2004}
J.~Puchalka and A.~Kierzek.
\newblock Bridging the gap between stochastic and deterministic regimes in the
  kinetic simulations of the biochemical reaction networks.
\newblock {\em Biophysical Society Biophysical Journal}, 86(3):1357--1372,
  2004.

\bibitem{srivastava2002stochastic}
R.~Srivastava, L.~You, J.~Summers, and J.~Yin.
\newblock Stochastic vs. deterministic modeling of intracellular viral
  kinetics.
\newblock {\em Journal of Theoretical Biology}, 218(3):309--321, 2002.

\end{thebibliography}

% flatex input end: [main_mixed.bbl]
%FLATEX-REM:\bibliographystyle{abbrv}
%FLATEX-REM:\bibliography{references}

%\clearpage
%\section*{To do list}
%\begin{enumerate}
%%\item Other splitting methods in the intro.
%%\item run the GTT problem with the mixed method and compare gains with the hybrid.
%\end{enumerate}

\end{document}

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: t
%%% End: 

% flatex input end: [main_mixed.tex]

1406.1989.extracted_bib
{poivariate J.~Ahrens and U.~Dieter. \newblock Computer methods for sampling from gamma, beta, {Poisson and bionomial distributions. \newblock {\em Computing, 12:223--246, 1974.
1406.1989.extracted_bib
{Anderson2012 D.~{Anderson and D.~{Higham. \newblock Multilevel {Monte {Carlo for continuous {Markov chains, with applications in biochemical kinetics. \newblock {\em SIAM Multiscal Model. Simul., 10(1), 2012.
1406.1989.extracted_bib
{anderson2007modified D.~F. Anderson. \newblock A modified next reaction method for simulating chemical systems with time dependent propensities and delays. \newblock {\em The Journal of Chemical Physics, 127(21), 2007.
1406.1989.extracted_bib
{KMC C.~C. Battaile. \newblock The kinetic monte carlo method: Foundation, implementation, and application. \newblock {\em Computer Methods in Applied Mechanics and Engineering, Volume 197(41-42):3386–3398, 2008.
1406.1989.extracted_bib
{gillespie2005 Y.~Cao, D.~T. Gillespie, and L.~R. Petzold. \newblock The slow-scale stochastic simulation algorithm. \newblock {\em The Journal of Chemical Physics, 122(1):--, 2005.
1406.1989.extracted_bib
{Cao2006 Y.~Cao, D.~T. Gillespie, and L.~R. Petzold. \newblock {Efficient step size selection for the tau-leaping simulation method. \newblock {\em The Journal of Chemical Physics, 124(4):044109, 2006.
1406.1989.extracted_bib
{CaoPetzold2006 Y.~Cao and L.~Petzold. \newblock Accuracy limitations and the measurement of errors in the stochastic simulation of chemically reacting systems. \newblock {\em Journal of Computational Physics, 212(1):6--24, 2006.
1406.1989.extracted_bib
{kurtzmp S.~N. Ethier and T.~G. Kurtz. \newblock {\em Markov Processes: Characterization and Convergence (Wiley Series in Probability and Statistics). \newblock Wiley-Interscience, 2nd edition, 9 2005.
1406.1989.extracted_bib
{GibBruck M.~A. Gibson and J.~Bruck. \newblock Efficient exact stochastic simulation of chemical systems with many species and many channels. \newblock {\em The Journal of Physical Chemistry A, 104(9):1876--1889, 2000.
1406.1989.extracted_bib
{gillespie_ssa D.~T. Gillespie. \newblock A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. \newblock {\em Journal of Computational Physics, 22:403--434, 1976.
1406.1989.extracted_bib
{gillespie_tau_leap D.~T. {Gillespie. \newblock {Approximate accelerated stochastic simulation of chemically reacting systems. \newblock {\em Journal of Chemical Physics, 115:1716--1733, July 2001.
1406.1989.extracted_bib
{gillespie2013 D.~T. Gillespie, A.~Hellander, and L.~R. Petzold. \newblock Perspective: Stochastic algorithms for chemical kinetics. \newblock {\em The Journal of Chemical Physics, 138(17):--, 2013.
1406.1989.extracted_bib
{Harris2006 L.~Harris and P.~Clancy. \newblock A “partitioned leaping” approach for multiscale modeling of chemical reaction dynamics. \newblock {\em J. Chem. Phys, Volume 125, 2006.
1406.1989.extracted_bib
{Haseltine2002 E.~Haseltine and J.~Rawlings. \newblock Approximate simulation of coupled fast and slow reactions for stochastic chemical kinetics. \newblock {\em J. Chem. Phys, 117(15), 2002.
1406.1989.extracted_bib
{li T.~Li. \newblock Analysis of explicit tau-leaping schemes for simulating chemically reacting systems. \newblock {\em Multiscale Model. Simul., 6(2):417--436 (electronic), 2007.
1406.1989.extracted_bib
{ourSL A.~Moraes, R.~Tempone, and P.~Vilanova. \newblock Hybrid {Chernoff tau-leap. \newblock {\em To appear in SIAM Multiscale Modeling and Simulation, 2014.
1406.1989.extracted_bib
{ourML A.~Moraes, R.~Tempone, and P.~Vilanova. \newblock Multilevel hybrid {Chernoff tau-leap. \newblock {\em arXiv:1403.2943, 2014.
1406.1989.extracted_bib
{Plyasunov2005 S.~Plyasunov. \newblock Averaging methods for stochastic dynamics of complex reaction networks: description of multi-scale couplings. \newblock {\em arXiv:physics/0510054v1, 2005.
1406.1989.extracted_bib
{Puchalka2004 J.~Puchalka and A.~Kierzek. \newblock Bridging the gap between stochastic and deterministic regimes in the kinetic simulations of the biochemical reaction networks. \newblock {\em Biophysical Society Biophysical Journal, 86(3):1357--1372, 2004.
1406.1989.extracted_bib
{srivastava2002stochastic R.~Srivastava, L.~You, J.~Summers, and J.~Yin. \newblock Stochastic vs. deterministic modeling of intracellular viral kinetics. \newblock {\em Journal of Theoretical Biology, 218(3):309--321, 2002.