macpan_ms.tex

\documentclass[12pt]{article}\usepackage[]{graphicx}\usepackage[]{color}
% maxwidth is the original width if it is less than linewidth
% otherwise use linewidth (to make sure the graphics do not exceed the margin)

\input{McMasterReport_preamble.tex}

\title{A compartmental model for epidemic parameter estimation and forecasting, with applications to SARS-CoV-2}

\author{Michael Li, Jonathan Dushoff, David J.\,D.\ Earn,\\
  Irena Papst, Benjamin M.\ Bolker\\
  McMaster University}

\begin{document}

\linenumbers
\maketitle

\begin{abstract}
Compartmental epidemiological models are widely used to understand, manage, and forecast the SARS-CoV-2 (COVID-19) pandemic. 
We introduce a new compartmental modeling framework that shares many characteristics with existing models, but includes a number of new and noteworthy features.
In particular, it includes a flexible structure based on the \emph{flow matrix} (the \emph{per capita} rates of transitions between compartments) that allows it to be used interchangeably for discrete or continuous time and for deterministic or discrete-state stochastic models; the capacity to set starting conditions based on the expected distribution of states during an exponential phase of the epidemic; automatic computation of $\R_0$ and the mean and dispersion of the generation interval for specified parameters; explicit structures incorporating the intensity of testing and delays between test administration and reporting; time-varying parameters based on breakpoints, spline bases, or external covariates such as cellphone-based mobility indices; and the ability to calibrate model parameters to multiple data streams such as case reports, hospitalization and ICU admission rates.
We demonstrate the model by calibrating it to multiple COVID-19 time series (positive tests, negative tests, hospitalizations, and deaths) for each of the Canadian provinces from 2020-02-27 to 2020-08-30.
We estimate epidemiological parameters, including the effective reproduction number $\R_t$ over the course of these early six months of the pandemic.
\end{abstract}

\david{Target journal?  Ideas: PLoS Comp Biol}
\ben{PLoS Comp Biol seems to be shooting very high for this.
  I'd be fine with PeerJ/PLoS ONE, but would be willing to try harder,
  e.g. some low-tier-but-not-terrible epi or modeling venue}

%%\vfill
%%\subsubsection*{Our group's COVID-19 research}
%% A brief summary, including publications by our group to date, is available at\\
%%\url{https://mac-theobio.github.io/covid-19/}.
%%\david{I updated Park et al 2020, which is now published in JRSI.  JD,
%%  please update that page with anything else you've published.}

%%\vfill

\tableofcontents

%%\vfill

%%\subsubsection*{\underline{Current date range for analysis}: \ \ date_range[1] -- date_range[2]}

%%\newpage


\section{Introduction}

SARS-CoV-2, the etiological agent of coronavirus disease 2019 (COVID-19), has been circulating in Canada since at least January 2020 \cite{onpr_200125}.
Response to the worldwide pandemic \cite{Li+20,Fauc+20} has been guided to a substantial extent by mathematical modelling \cite{Flax+20}.

%% Virological testing and contact tracing have been conducted to varying
%% degrees in different jurisdictions\david{REFS}.  While testing is
%% imperative for surveillance, its value for control is less clear.
%% Agent-based simulations indicate that the combination of testing and
%% contact tracing can be very effective in mitigating spread of
%% SARS-CoV-2 \cite{Ng+20}.  Emerging research is demonstrating that
%% saliva-based tests, as opposed to naso-pharyngeal swabs for RT-qPCR,
%% may be more effective, cheaper, and viable for the public to
%% self-administer \cite{Will+20,Wyll+20}, potentially making very
%% large-scale, frequent testing possible.

Here, we present a compartmental framework that was developed over the course of the pandemic. 
It incorporates the standard epidemiological compartments required to model COVID-19 as well as compartments that track health-care utilization and COVID-induced
mortality. 
Case reporting can be modeled either as a time-delayed convolution of incidence or by enabling a factorial expansion of the model that accounts for the testing status of individuals \cite{Fris+20}. 
The model also allows for time-dependent variation in any rate parameter --- particularly the transmission rate --- and allows this variation to be indexed by external covariates such as cellphone-based mobility metrics, or to follow smooth (spline) curves over time. 
The same model structure can be run as an ordinary differential equation, a discrete-time deterministic model, or a discrete-time stochastic model.
Finally, the model can be used to
(1) simulate specific scenarios for planning purposes;
(2) calibrate parameters to match multiple input time series such as hospital admissions or occupancy, cases, or deaths; or
(3) forecast future epidemic dynamics on the basis of past calibration.

\section{Methods}

\subsection*{Compartmental structure}

The epidemiological structure of the model is based on a susceptible-exposed-infectious-removed (SEIR) model with additional compartments reflecting the biology of COVID-19 and the structure of the health-care system. 
The COVID-19-specific compartmental structure of the model resembles many other COVID-19 models \cite{childs2021impact, tuite2020mathematical, kainChopping2021} in separating infectious individuals into sub-compartments reflective of the epidemiology of COVID-19; it additionally adds compartments for hospitalized individuals in acute care or intensive care. 
All symptomatic individuals are presumed to have undergone a period of pre-symptomatic infectiousness (\code{p}). 
Infections can be asymptomatic (\code{a}), mildly or moderately symptomatic (\code{m}) or severely symptomatic (\code{s}); all individuals with severe symptoms go to the hospital (acute care or ICU), or die before reaching hospital (e.g., in long-term care facilities).
Some fraction \david{add parameter symbol in brackets} of individuals who go to the ICU die. Recovered individuals (\code{R}) are assumed to be immune.  
The model includes additional compartments that facilitate book-keeping: cumulative hospital admissions (\code{X}), individuals in acute care after discharge from ICU (\code{H2}), and cumulative deaths (\code{D}) (\cref{fig:flowchart}).
(The version of the model discussed here was developed and used before vaccines were available and introduction of variants of concern and reinfections, we have since expanded the model to include other relevant compartments.)

\david{It is weird that we do not have a table listing all the parameters of the model.
I guess that's why the symbol for the proportion in ICU who die was not given above.}

This version of the model assumes homogeneous mixing --- all classes of infectives contribute additively to the force of infection (the \emph{per capita} infection rate of susceptibles).
We assumed hospitalized individuals do not contribute to the force of infection. \mike{re-write this point nicely? it is a strong assumption.}  
We did add one feature to the model to account for heterogeneity in susceptibility in the population, which we typically imagine is driven by heterogeneity in exposure (e.g., front-line and essential workers will be infected earlier), but could also be influenced by genetic, immunological, or other factors.
A \emph{phenomenological heterogeneity} parameter $\zeta$, modifies the force of infection by a factor of $\left(S(t)/N\right)^\zeta$. Since $0 < S/N < 1$, a positive value of $\zeta$ will make the force of infection decrease as the remaining fraction susceptible decreases, capturing the fact that the most susceptible individuals tend to be infected first, after which
the average level of susceptibility in the remaining population decreases.\david{revised that last sentence, which I kept misreading} While modelers have considered the dynamics of epidemiological models incorporating incidence functions of this general form \cite{WilsWorc+45,Liu+87}, previous attempts to model phenomenological heterogeneity have used transmission proportional to an exponential function of prevalence (i.e., $S I/N \times e^{-(I/N)^n}$ where $n$ is a shape parameter) rather than using a power law \cite{Will+06,Gran+09}. \ben{JD please check that this makes sense/expand where necessary; check Dwyer? Williams et al 2006: ``To allow for heterogeneity in sexual behaviour and to fit the observed asymptotic prevalence of infection, the transmission parameter takes the value $\lambda_0$ at the start of the epidemic and declines exponentially at rate $\alpha$ times the prevalence of infection.'' Granich et al. say ``To allow for heterogeneity in sexual behaviour and for the observed steady state prevalence of HIV, we let the transmission decrease with the prevalence, P. If n=1, the decrease is exponential; if $n=\infty$, the decrease is a step function. Both have been used in previous models (5, 29).'' Granich et al. also cite \cite{williamsHIV2006} (their ref. 29), but I don't see a transmission model in there anywhere ...}

%\begin{figure}[ht!]
\begin{figure}
% \includegraphics[width=\maxwidth]{figure/flow.chart.crop.pdf}
\includegraphics[width=\maxwidth]{pix/Macpan_Base_Epi-model.png}
\caption{Flow chart for basic compartmental mechanistic transmission model. 
Compartments: \code{S} (susceptible), \code{E} (exposed), \code{Ia} (asymptomatic infection), \code{Ip} (presymptomatic infection), \code{Im} (mild/moderately symptomatic infection), \code{Is} (severely symptomatic infection), \code{H} (hospitalized [acute care]), \code{ICUs}/``ICU Survive'' (ICU with prognosis of survival), \code{ICUd}/``ICU Die'' (ICU with prognosis of death), \code{H2}/``Hospitalized Recovery'' (acute care after ICU stay). 
%Compartments denoted by rectangles are accumulators, primarily used in the condensation step to compute incidences: \code{R} (recovered), \code{D} (dead), \code{X} (accumulator for cumulative hospital admissions).
}
\label{fig:flowchart}
\end{figure}

%% \begin{figure}[ht!]
%% \definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}
%% \includegraphics[width=\maxwidth]{figure/ratematrix.pdf}

%% \caption{Flow matrix for basic compartmental mechanistic transmission
%%   model.  Gray indicates non-zero flows: compartments as in
%%   \cref{fig:flowchart}.  \ben{do we want this? It's equivalent to
%%     the flowchart, could be used to show differences in magnitude (via shading)
%%   if we sorted out the details} \jd{I think we don't.}
%%   }
%% \label{fig:flowmatrix}
%% \end{figure}

Internally, the model is defined by a \emph{flow matrix} $\mathbf{M}$, the elements of which specify the \emph{per capita} rates at which individuals move from one compartment to another.
For the basic model, the only element of the flow matrix that needs to be recomputed at each time step is the incidence (flow from $S$ to $E$); all other rates are piecewise constant (they are adjusted instantaneously when policies change, for example).
%%
This set-up allows for considerable flexibility:
\begin{enumerate}
\item For a numerical differential equation solver, we need the time derivatives of each compartment. The absolute rates are computed by columnwise multiplication by the state vector $\mathbf s$ (${\mathbf F}_{ij} = \mathbf{M}_{ij} \mathbf{s}_j$).
The gradient is the difference between the total flows into (column sums of $\mathbf F$) and out of (row sums of $\mathbf F$) each compartment.
\item For a discrete-time model, we compute the vectors of outflows and inflows as above, but simply compute the new states as (original state + (inflow - outflow)$\Delta t$).
\item We can also use a discrete-time model with a \emph{hazard correction} 
  \david{Is there a published paper that uses this terminology?  or at least a paper
    we can cite that uses the idea?  I believe that using this to make a discrete-time
    model is discussed by Hoppensteadt in an old book that I could find.}
  to the flows to address the possibility that a state will go negative due to an overly large outflow. Instead of the total \emph{per capita} outflow of a compartment being equal to the sum of \emph{per capita} flow to each other compartment ($\textsub{f}{tot} = \sum_i {\mathbf M}_{ij}$), we let the total outflow be 
\begin{equation}\label{eq:hazard}
\textsub{f'}{tot} = 1-\exp(-\textsub{f}{tot} \Delta t) \,; 
\end{equation}
the individual flows are then adjusted by a factor of $\textsub{f'}{tot}/\textsub{f}{tot}$.  \david{Notation is not ideal here.  $\textsub{f}{tot}$ depends on $j$; perhaps $f_j$ would be better?  Why $f$?  Maybe
  $M_j$?  Also $f'$ conveys ``derivative''.  Perhaps $\tilde f$ would be better?}
This adjustment accounts for the effects of depletion during the course of a time step.
\item Finally, if we choose to run a fully stochastic simulation, the flow matrix $\mathbf M$ is what we need to sample the flows between compartments as Euler-multinomial deviates \cite{breto+09}, which take the hazard correction \eqref{eq:hazard} and use it to compute probabilities for draws from binomial or multinomial random deviates.
  \david{Maybe explicitly write the formula for the probabilities
  in terms of the hazard correction.}
\end{enumerate}

While the flow matrix description is convenient for most epidemiological dynamics, there are a few epidemiological processes that are more naturally captured by absolute rather than \emph{per capita} rates, for example (1) when intensities of public health interventions such as numbers of tests or vaccines administered are reported by public health agencies; or (2) in models incorporating births or immigration.
We have typically handled the former case by taking observed vaccination or testing rates and dividing them by current compartment sizes in order to set the relevant entries in the flow matrix; we have not yet tried to build models including inflows from outside the system (the latter case).\david{The phrasing ``not yet tried'' feels odd in a paper.  Not sure how to express this.  Maybe ``the version of the software described here does not include inflows from outside the system (the latter case).''}

We typically compute the model dynamics deterministically, as a discrete-time model with a hazard correction \eqref{eq:hazard} (option 3 above).  Once the trajectories are computed, we reduce the full state vector to a more convenient, collapsed state vector in a step we call \emph{condensation}, for example by summing all of the infectious compartments to a single $I$ state vector, or collapsing the different acute-care ($H$, $H2$) or ICU ($\textrm{ICU}_\textrm{s}$, $\textrm{ICU}_\textrm{d}$)
\david{We should be consistent wrt font used for compartments.  We
were using \macro{texttt} earlier.}
compartments.
As well as allowing us to visualize model results more conveniently, condensation also allows us to compare the simulated state vector to available data streams.
In addition to summing compartments, we can also compute incidences as time-lagged differences of accumulator compartments (for example, differencing accumulated deaths $D$ to derive a mortality rate) or perform more complicated operations such as convolution\david{maybe cite Goldstein et al 2009, PNAS}. 
Our main use of convolution is to convert incidence---the force of infection (FOI) multiplied by the number of susceptibles---to a case-reporting (CR) time series: 
\begin{linenomath*}
\begin{equation}\label{eq:CR}
\textrm{CR}(t) = \sum_i \phi(i) (\textrm{FOI}(t-i) S(t-i))\,,
\end{equation}
\end{linenomath*}
where we typically set $\phi(i)$ to be a Gamma distribution with
moments chosen to match empirical estimates of case-reporting delays.
\mike{Do we have/need a formula for FOI?}
\ben{how did we pick these values? 
Probably don't have a formal ref but some statement of what we used to guess that this was reasonable would be good ...}.
\ben{done, but still not clear that there are a lot of such estimates in the literature \ldots}\david{I added ``empirical'', which can cover unpublished estimates.}
When computing case reports from incidence we also assume a case-report proportion \code{c\_prop} to account for the fact that the majority of COVID infections are never reported \cite{Doug+2020}; this value is usually calibrated from data.

After condensation, the model also allows us to add observation error, which we typically simulate from a negative binomial distribution with a variable-specific dispersion parameter \cite{linden2011using}.

\subsection*{Expansion to accomodate testing}

At the cost of additional complexity, we can add explicit testing compartments to the model.  In the simpler version of the model, we assume that a specified fraction of infections are reported as cases (or we calibrate this fraction from joint data on cases and hospitalizations), and impose a distributed delay between infection and reporting via convolution \eqref{eq:CR}.
Here, we instead expand the susceptible and all infected compartments factorially to include the possibilities that individuals in those epidemiological classes have one of four testing statuses:
\begin{itemize}
\item untested;
\item tested and awaiting negative results;
\item tested and awaiting positive results; or 
\item tested positive.
\end{itemize}
After receiving negative results (whether true negatives, from individuals in $S$, or false negatives, from individuals in one of the infectious compartments), individuals cycle back to the corresponding “untested” sub-compartment, since they can be tested again. After receiving positive results, they remain in the positively tested sub-compartment; depending on the model parameters, their transmission may be reduced due to self-isolation (controlled by the parameter \code{iso\_t}).  We assume here that people waiting for test results do not isolate.
Gharouni \emph{et al.}\/ \cite{Ghar+22} thoroughly analyze the epidemiological consequences of this structure in a simpler framework that expands a basic SIR model rather than the COVID-specific compartmental model (\cref{fig:flowchart}) as the foundation.
Individuals may progress between epidemiological compartments (e.g., becoming infectious or recovering) while awaiting the results of tests. In general, progressing individuals move to the same testing subcompartment, e.g., from \code{Im}-negative-waiting to \code{R}-negative-waiting.

We define a weight vector $w$ across epidemic compartments that gives weight $w_\textrm{a}$ to asymptomatic classes (\code{S}, \code{E}, \code{Ia}, \code{Ip}, \code{R}) and 1 to symptomatic classes (i.e., all other classes).
If $W$ is the weighted sum of compartment occupancies $X_i$, i.e., 
\begin{equation}
W = \sum_i w_i X_i/N \,,
\end{equation}
then for a daily \emph{per capita} testing rate $\rho$ we might expect that the corresponding \emph{per capita} testing rate in compartment $i$ would be ${\mathcal T}_i = \rho w_i/W$. 
\david{We should try to keep notation consistent with \cite{Ghar+22}.  I changed
$T_i$ to ${\mathcal T}_i$.  Not sure that's sufficient.}
However, under some extreme conditions (if testing is so extreme that few untested symptomatic people are left) this formulation can allow the \emph{per capita} testing rate to explode. To address this problem we add a maximum daily \emph{per capita} testing rate $\tau$ to the model such that
\begin{equation}
{\mathcal T}_i = \frac{\rho \tau w_i}{\tau W + \rho} \,;
\end{equation}
see Appendix A.5 of \cite{Ghar+22} for details.  The testing flows out of each untested subcompartment are divided into flows to ``positive waiting'' and ``negative waiting'' compartments according to the infection status of the relevant compartment and specificity/sensitivity parameters (specificity/sensitivity refers to the probability that a truely negative/positive individual tests negative/positive).  \david{phrasing could probably be improved, but I felt we should define specificity and sensitivity} If tests are assumed to be perfectly specific and sensitive, then all individuals from non-infectious compartments (\code{S}, \code{E}, \code{R}) enter the ``negative waiting'' compartment and those from infectious compartments enter the ``positive waiting'' compartment.

We assume that all individuals admitted to the hospital for COVID-19 are immediately tested. (These tests are not included in the accounting of test distribution above, but in the COVID setting they represent a small fraction of the overall number of tests administered.)

\begin{figure}
  \includegraphics{figure/testFlow.1.Rout.pdf}
  \caption{Testing flow. Every epidemiological compartment is subdivided into the four subcompartments shown.
    Black arrows represent flows between subcompartments due to testing processes (test administration, reporting of tests); blue arrows represent progression between epidemiological compartments.
    Dashed arrows represent the accumulation of negative and positive test reports, which can be compared against data.}
  \label{fig:testing_flow}
\end{figure}

\subsection*{Model parameterization}

The model allows for time-varying, piecewise-constant changes in any parameter. 
Our early analyses focused on changes in time-varying effective reproductive number (\Rt) due to behaviour change and non-pharmaceutical interventions, which we model by changing the transmission rate.
The transmission rate $\beta(t)$ is taken to be a time-varying function of the form $\beta_0 \beta_1(t)$ where $\beta_0$ is the baseline value for transmission from symptomatic individuals.
Individuals in different symptomatic classes (presymptomatic, asymptomatic, mild, severe) may have their transmission modified by a specified multiplier (given as model parameters); we assume that hospital transmission is negligible.
\david{Can we find a reference to support most transmission occuring in the community / outside hospitals?}
The time-varying (relative) transmission $\beta_1(t)$ can incorporate a variety of different effects, one at a time or in combination:
\begin{itemize}
\item abrupt (piecewise) changes on specified dates when control measures are known to have been implemented;
\item proportional to a power of observed mobility or some other exogenous proxy for contact behaviour:
\begin{equation}
  \beta_1(t) \propto \log(M(t))^{\pmob} \,, 
\end{equation}
where $M(t)$ is the relative mobility index at time $t$, for some power $\pmob>0$;
\item according to an arbitrary spline curve, i.e., a linear combination of components of a B-spline basis.
\end{itemize}
All of these sub-models for temporal change in the transmission rate can be subsumed under a single log-linear model:
\begin{equation}\label{eq:betamodel}
\log \beta(t) = \log \beta_0 + \XX\boldsymbol{c}
\,,
\end{equation}
where $\XX$ is a model matrix 
\david{not ideal that we have $X_i$ meaning something else entirely above}
that can contain any combination of covariates and $\boldsymbol{c}$ is a vector of covariates determining differences in the log of the relative transmission rate.
\david{Only $\boldsymbol{c}$ can vary in time, right?  We should therefore
write $\boldsymbol{c}(t)$ in \cref{eq:betamodel}.}
In particular, piecewise breaks correspond to indicator variables for which period an observation falls in; the mobility model corresponds to a column containing the log of relative mobility; and the spline model corresponds to a set of columns containing the basis vectors of a B-spline basis with specified knots.

In practice, when we use the model without incorporating covariates we adopt a simpler strategy of providing a list of the breakpoints and the parameters that change at those breakpoints --- a special case of the more general log-linear model. In the models illustrated below, we pick breakpoints to denote several time periods $\{P_1, P_2, \ldots\}$ and construct $\XX$ to allow the effects of mobility, and the baseline contact rate, to change at breakpoints.
Specifically, for each breakpoint we define a logistic transition curve 
\begin{equation}
S(j) = \frac{1}{1+e^{[t-t_\textrm{brk}(j)]/s}} \,,
\end{equation}
where $s$ (set to 3 days) defines the speed of transition. Because the baseline transmission rate is included in the model (as parameter \code{beta0}), we do not need to include an intercept in $\XX$; the first column is $\log(M(t))$, the dependence of (log) transmission rate on (log) mobility during the first period. For $j \ge 1$, subsequent columns $2j$ and $2j+1$ of $\XX$ are defined as $S(j)$ and $S(j) \log M(t)$, respectively. The parameters associated with these columns denote the change in baseline contact rate and the change in the effect of mobility between periods $j$ and $j+1$.
  
\subsection*{Derived parameters and parameter setting}

It is useful to be able to compute several quantities derived from the parameters of a model, in particular:
\begin{enumerate}
\item the dominant eigenvector of the system in the exponential growth phase
  (as described below, we use this eigenvector to set sensible initial conditions for all state variables, including those that cannot be observed);
\item the basic reproduction number \Rzero and the intrinsic growth rate $r$;
\item the mean and coefficient of variation of the generation interval. 
\end{enumerate}
In principle we could compute these values directly from the flow matrix, by constructing the Jacobian matrix and the next-generation matrix and performing the appropriate eigenvector/eigenvalue calculations (as discussed by \cite{VandWatm02} for differential equations and \cite{Casw00} for discrete-time systems). However, we found it simpler to derive these values by simulation.

To compute the eigenvector, we run a simulation where we set the outflow from the susceptible compartment to zero (while maintaining the \emph{inflow} from $S$ to $E$), which mimics the dynamics near the disease-free equilibrium where susceptible depletion is negligible. After running the simulation for a long time (100 steps by default), the state vector is close to the eigenvector.  
We use this value to set starting conditions when starting near the beginning of the epidemic. 
It is easy to specify a scalar value (say, 1\% of the population) to indicate the initial size of the epidemic; by distributing these individuals according to the calculated eigenvector, we reduce numerical instability at the beginning of our simulations. 

To compute the other summaries (\Rzero, $r$, and moments of the generation interval) we rely on the fact that the software saves the force of infection at each time step. 
We simulate the progression of a \emph{single} individual through the infection process, i.e., setting the population size to 1 and setting the initial state to $E=1$ with all other compartments empty. 
Simulating deterministically then generates a time course of the probability that an individual is in any given box at a particular “age of infection”, and thence the expected force of infection generated by a single infectious individual through time. 
This vector $K(t)$ is the same as the transmission kernel in a renewal equation \cite{Cham+18}. 
We can easily compute 
\begin{equation}\label{eq:Rzero}
\Rzero = \sum_t K(t) \,, 
\end{equation}
mean generation interval  
\begin{equation}\label{eq:Gbar}
\bar G = \sum_t \frac{t\,K(t)}{\Rzero} \,, 
\end{equation}
and generation interval coefficient of variation,
\begin{equation}\label{eq:GCV}
\textrm{CV}(G) = \frac{\sqrt{\frac{1}{\Rzero}\sum_t K(t) (t- \bar G)^2}}{\bar G} \,.
\end{equation}
The growth rate can be computed by numerically solving for $r$ in the Euler-Lotka equation,
\begin{equation}\label{eq:Euler-Lotka}
\sum_t  K(t) e^{-r t} = 1 \,.
\end{equation}

\jd{Suggest dropping the detailed equations for Gbar and CV, but instead adding $g(t)$ and saying that we also calculate its mean and CV.}
\jd{Suggested easy ref is ChamDush15: http://dx.doi.org/10.1098/rspb.2015.2026 (instead of current Cham18)}
\ben{JD: please go ahead and implement this change}

In addition to its use in summarizing a given set of parameters, the computation of $\Rzero$, $r$, and the generation interval is useful as an initial step in calibrating the model. 
Estimates of these summary statistics are more broadly available \cite{park2020reconciling}, and more epidemiologically relevant, than the more specific mechanistic parameters describing the relative infectiousness and duration of each of the different infectious compartments (although this detailed information is still important for determining the effectiveness of interventions like contact tracing). 
We typically start with mechanistic parameters gathered from the literature and pre-calibrate them to specified target values of $r$ (which is easy to estimate from the observed initial growth rate of the epidemic in a region) and the mean generation interval by adjusting $\beta_0$ (baseline transmission) and simultaneously scaling the values of all of the epidemiological transition rates ($\sigma$, \textsub{\gamma}{s}, \textsub{\gamma}{m}, \textsub{\gamma}{a};
see \cref{tab:litparm})
by a single factor until the target values are achieved.

\subsubsection*{Calibration}

Once we have run a deterministic model simulation for a particular set of parameters (including a starting number infected, distributed across non-susceptible classes according to the exponential-phase eigenvector computed as described above) and we have some set of time series data to calibrate against (for example, case reports and hospital admissions), we can calculate a log-likelihood (e.g., \cite{Bolk08}). 
We assume that every observation is independently negative binomially distributed, with a series-specific estimated dispersion parameter (i.e., the variability in cases, hospital admissions, \emph{etc.}, will differ). 
We use standard nonlinear optimization algorithms built into R, such as Nelder-Mead, to find the maximum likelihood estimates; when we have had difficulty with numerical instability, we have performed an initial fit with differential evolution \cite{Mull+11} followed by a final fit with Nelder-Mead. \mike{Do we need a line about process/dynamical error here? Say we have it or we are not including it.} 

A link function can be added for any parameter in the model to constrain it to a sensible domain; the user specifies this by adding an appropriate prefix to the name of the parameter in the list of starting values for parameters to be calibrated. 
For example, specifying \code{log\_beta0 = -1} would specify that the baseline transmission parameter $\beta_0$ should be calibrated on the log scale, ensuring that the value of $\beta_0$ is always positive, and using a starting value of -1 on the log scale (i.e., initial $\beta_0 = \exp(-1)$). 
Specifying \code{logit\_nonhosp\_mort = -0.2} would specify that the value of \code{nonhosp\_mort} (the fraction of mortality that occurs outside hospitals) should be fitted on a logit, or log-odds, scale, ensuring that it is bounded between 0 and 1, and using a starting value of $\beta_0 = \textrm{logit}^{-1}(-{0.2}) = 0.45$.

The model includes a general framework for adding a prior probability distribution for any parameter, using any distribution available in R. 
For example, \code{dbeta(nonhosp\_mort, 2, 2)} would specify a $\textrm{Beta}(2,2)$ prior for \code{nonhosp\_mort}.
We do not need to adopt a fully Bayesian framework to make use of priors; instead, we can think of them as convenient regularizing factors to keep the model-fitting process numerically stable. If we do want to be Bayesian, then the fitting procedure described above will return maximum \emph{a posteriori} (MAP) parameter values, not a sample from the full posterior distribution as is standard with frameworks that use Markov chain Monte Carlo.
\david{A standard ref is needed for MCMC.  \cite{Bolk08} or something else?}

For the province of Ontario, Canada, in 2020
we calibrated to deaths, and new confirmations, all of which were available publicly. 
In the expanded model, we can include time series of both positive and total tests ( in our calibration. 
Reported new confirmations (daily positive tests) are the most reliable and voluminous source of epidemic information.
Unfortunately, they are also subject to many inherent biases, including substantial variation over time in \term{testing intensity} (\ie tests \emph{per capita} per day).

We simultaneously estimate the temporal pattern of the transmission rate [\cref{eq:betamodel}] (two mobility intercepts and three slopes), and several basic model parameters (\cref{tab:estparm}).
Other model parameters are taken from the literature (\cref{tab:litparm}).

\begin{table}
  \centering
  \input{litparm_table.tex}
  \label{tab:litparm}
\end{table}

We do not include hospitalization and ICU occupancy in our calibration because our model assumes all severe cases go through ICU, whereas during the early pandemic period
that we focus on here, many severe cases actually occurred entirely in Long Term Care Facilities (LTCFs); furthermore, capacity limitations in the health care system lead to ceilings on hospitalization that are
not represented in our model.
\jd{I removed earlier comments, and simplified by removing stuff about why we didn't want to use an inappropriate model.}

\subsubsection*{Forecasting}

Our calibrations yield values for the parameters of our deterministic model (\cref{tab:estparm}).
Using the calibrated parameters and the estimated covariance matrix of the sampling distribution of the parameters, we draw many (typically 1000) sets of parameters from a multivariate normal distribution \cite{Bolk08,krinskyThree1991a} and feed them into the deterministic model to generate an ensemble of forecasts. We extend the end date of the calibration window (i.e., the last observed data point for calibration) to create a forecast window. In the forecast window, the user can either investigate scenarios by inputting covariates (i.e., relative mobility and/or testing rates, depending on the model) or assume that the last set of inputs remains constant through the forecasting window to obtain a
\emph{status quo} forecast.

\subsection{Data sources}

\subsubsection*{Public COVID-19 data for Ontario, Canada}

Daily reported data on COVID-19 testing and outcomes (confirmed cases, deaths, hospital and ICU occupancies) are publicly available for Ontario on the official provincal website:
\url{https://data.ontario.ca/dataset/f4f86e54-872d-43f8-8a86-3892fd3cb5e6/resource/ed270bb8-340b-41f9-a7c6-e8ef587e6d11/download/covidtesting.csv}.
However, this dataset does not include testing counts before 15 April 2020.  Since the beginning of the COVID-19 pandemic (and until the present at the time of writing), one of us (ML) has maintained a public web site containing Canadian COVID-19 data at the provincial level.  Data are frequently downloaded from a variety of sources and cleaned.  See \url{https://wzmli.github.io/COVID19-Canada/}.

\hypertarget{Mobility data}{}
\subsubsection*{Mobility data}

We use mobility data from Apple\footnote{\url{https://raw.githubusercontent.com/ActiveConclusion/COVID19_mobility/master/apple_reports/applemobilitytrends.csv}} and Google\footnote{\url{https://www.gstatic.com/covid19/mobility/Global_Mobility_Report.csv}}.  From these data we derive a \term{relative mobility index} using the ``driving'' index from Apple and the ``retail and recreation'' and ``workplaces'' indices from Google; we compute a 7-day moving average of these indices, rescale all of them to have a baseline (pre-pandemic) value of 1.0, and average the three indices (equally weighted), to obtain an overall index for each day.

\david{The ``important dates'' section is in \texttt{cuts.tex}.  Do we really not use those dates at all?}
\mike{We didn't use it because mobility is handling most of the changes in transmission.}
\mike{We replaced it with the important dates when we gave up on mobility later on. I do prefer mobility for this paper for a more generalized approach to incorprating different exogenous data (I guess same arugement applies to important dates as well).}

\begin{figure}[ht!]
\includegraphics[width=\maxwidth]{figure/ontario_mobility.Rout.pdf}

\caption{Ontario mobility. \ben{need more description/standalone: 7-day average, etc.}
}
\label{fig:Ont_mobility}
\end{figure}

\FloatBarrier

\hypertarget{sec:Results}{}
\section{Results}

We fitted two different models to the Ontario data demonstrate the capability of this modeling framework. The first model is our simplest model, following the schematic illustrated in Figure \ref{fig:flowchart}. This model calibrates to new confirmation and death time series and includes mobility with two additional breaks and phenomenological heterogeneity. We will refer to this as the ``base model''.
The second model extends the base model by incorporating the testing structure illustrated in Figure \ref{fig:testing_flow}. We will refer to the second model as the ``testify model''. 
We did not explicitly include break dates (i.e., time-varying transmission rates) to match different public health measures and restrictions in place within the fitting window, but instead included them in the mobility index. The transmission rate is a function of time-varying intercepts and slopes on the mobility index, which is a sensible proxy for changes in behavior in response to public health restrictions. Two change points were included to allow the transmission rate to have a different intercept and slope with respect to the mobility index, capturing the adaptiveness of the mobility patterns.
Both models were used to forecast 30 days ahead, under the assumption that the mobility index would remain constant at the last observed value for the forecast period.
\david{That was true in 2020 but is no longer true.
  We should rephase.  We could perhaps say ``was a good proxy\dots\ during the focal time period'' or something like that, but we then need to say this stopped being true eventually and some reference(s) to support this would be very helpful.} 
\mike{Maybe add it in the discussion?}
\jd{What isn't true? I feel the text has left these comments behind?}

\subsection{Base model}

Figure \ref{fig:Ont_calibration_base_forecast} shows the fit to new confirmations and death time series from Febrary 24 2020 to August 30 2020. This time frame captures the first wave of the COVID outbreak and the beginning of the second wave in Ontario. The model capture to fit the new confirmations data well within the fitting window but not for the death time series. The model has adaqualt flexibility to capture the effective transmission parameter with mobility data and changes in mobility effects for each mobility break on Apr 1st and Aug 7, and mortality parameters were assumed to be constant at the optimized level resulting in a overestimation of deaths. 

The 30 day forecast is expected to increase; when over-laying actual observed time series for that time-frame, the observed time series are lower than predicted but within 95\% prediction interval for the new confirmations and missed all the death observations. The parameters for mortality were assumed to be constant throughout the fitting and forecast window. 

Figure \ref{fig:Ont_calibration_base_rt} shows the R0t and Ret from the base model. \mike{thought on the interval band on Rts?}

\begin{figure}[ht!]
\includegraphics[width=\maxwidth]{figure/ontario_base_forecast_plot.Rout.pdf}

\caption{Ontario base forecast. \ben{need to explain what's going on here; at the very least the caption should say (even if it's kind of obvious) that gray points represent calibration points, red are forecast points.  Why is the death forecast so far off?  Should we just leave it out (probably)? (WZ ...)}}
\jd{More comfortable leaving it out after we figure it out. It's bizarre. I could be convinced that it's water under the bridge, but it's worth a look.}
\label{fig:Ont_calibration_base_forecast}
\end{figure}

\begin{figure}[ht!]
\includegraphics[width=\maxwidth]{figure/ontario_base_rt_plot.Rout.pdf}

\caption{Ontario base reproductive numbers (Reff and Rt).}
\label{fig:Ont_calibration_base_rt}
\end{figure}

\begin{table}
\centering
\input{combo_table.tex}
\caption{
  \label{tab:estparm}
  Parameter estimates for base model calibration. \ben{JD/DE: should some parameters (change in transmission etc.) be back-transformed? Are parameter names comprehensible?  Should we add confidence intervals/std errs? We should explain that dispersion model is converging to Poisson; do we need to justify/explain the very high value of $\zeta$ (phenom het)?}\david{We definitely need to comment on
    the value of $\zeta$; it seems absurdly high.}
  \david{CIs would be nice.}
}
\end{table}

\FloatBarrier

\subsection{Testify model}

The base model is relatively simple but does not account for testing practices. Testing strategies changed frequently over the course of the pandemic due to availability of testing facilities and shifting regulations on eligibility for testing.

Figure \ref{fig:Ont_calibration_testify} shows the results from calibrating the model to positive tests, with the testing flow incorporated and with observed testing rates fed into the model as a known series.
The variation in the predicted line is driven by high-frequency changes in testing rate, including day-of-week effects.
It also shows the daily testing in Ontario in the calibration time frame. The vertical line is the cut-off between the calibration window (before September 1st 2020) and 30 days beyond it for the forecast window. The testing intensity over the forecast window is assumed to be constant at the last observed testing intensity on August 30th 2020 (since the actual increase in testing that occurred during the forecast window was not known at the time of the forecast). The forecast under predicts the number of new positive cases in the forecast time-frame \ref{fig:Ont_calibration_testify}. The actual testing intensitive did increase during the forecast period. 

\david{I suggest using open circles during the forecast window, as we do for observed cases and hospitalizations.}

\begin{figure}[ht!]
\includegraphics[width=\maxwidth]{figure/ontario_testing_plot.Rout.pdf}

\caption{Calibration of positive tests and testing intensity}
\label{fig:Ont_calibration_testify}
\end{figure}

\begin{figure}[ht!]
\includegraphics[width=\maxwidth]{figure/ontario_testify_forecast_plot.Rout.pdf}

\caption{Ontario calibration with testing flows. \mike{do we want to use the raw testing intensity or smooth it with 7dayMA?}}
\label{fig:Ont_calibration_testify_forecast}
\end{figure}

\mike{We (I) probably should make another forecast using the testing intensity of the forecast window. Thoughts?}
\david{Yes it would be good to see how using the correct testing intensity affects the forecast.  This needs to be expressed carefully to avoid confusion.  We need to get across that this emphasizes how much better the forecast would be with advance knowledge of how testing intensity would change over the forecast period.}

\FloatBarrier

\section{Discussion}

We present a compartmental framework developed in real time during the SARS-CoV-2 pandemic, first to make simple projections, and later to fit epidemic data, combine different data streams, and make forecasts for a variety of public-health agencies and governmental advisory bodies. Valuable lessons were learned in this process that should be consolidated to provide a better foundation for future pandemics. The modeling approaches used here, and the challenges faced, have already been incorporated into newer, more flexible tools, and development is continuing.

\jd{Can somebody brief me on how these fits relate to what we actually did in real time?}
We show two example fits for the province of Ontario early in the pandemic, from models that vary in complexity.  Our simple base model is typical of compartmental epidemic models used to study infectious disease epidemics.  The testify model is a practical expansion of the base model that allows us to take advantage of data on time-varying testing practices and policies.

\mike{Talk about the assumptions of the forecast period for each model. Freezing vs extending/predicting mobility and testing intensity.}
\jd{Less interested in these details. What do we \emph{learn} from the fits and the from the process? What's the contrast between the two methods.}

\jd{Limitations should be short (or absent? I don't care about the limitations of these particular fits), and the Conclusion (or just the end of the Discussion should talk more about what we've learned, what we're doing about, and what should continue to be done about it.}

\subsection{Limitations}

Our model assumes homogeneous mixing of the population and excludes population-level heterogenetity such as age, spatial/geographical contact structures.
Calibrating to multiple time-series are challenging; it has to respect the biological process/relationship as well as the flexibility to change in various times. 
\david{The following sentence needs to be rewritten.  I don't understand what was intended to be expressed.}
We did not include hospitalization in the calibrations due to reporting from capacity issues which became a serious issue with the admissions and occpuancy underreporting the severity of the pandemic.
In addition, long-term care facilities (LTCFs), where many elderly people in Canada have died without going to hospital, have not been modelled explicitly.  
While we do not anticipate any qualitative differences in results, explicitly creating compartments and calibrating to data for LTCFs would likely allow us to match ICU occupancy and forecast
pressure on ICUs more accurately but will not resolve the issues on hospital intensity.
\david{the last phrase needs to be rephrased.   What is ``hospital intensity''?}
\mike{How busy the hospitals get. I.e. addmission only shows how many the hosp admit but not how many people how actually showed up that needs treatment.}

When forecasting into the future, many assumptions have to be made for both exogenous input parameters. For example, future mobility patterns are assumed to be the same as the last observed patterns for both models, and similarly for testing intensity in the testify model. As observed in the base model, a fixed mortality ratio is insufficient to capture the mortality patterns. Even if they are captured accurately, there is a possibility that they may differ in the forecast window.
 
\subsection{Conclusions}

\david{This conclusion paragraph leaves a lot to be desired.  I think we should discuss what we think the key take-home messages should be, and then write a new conclusion paragraph.  We should also state that we have continue to develop the software during the pandemic, and will describe further innovations in subsequent papers.}

SARS-CoV-2 continues to be a global burden after 3 years. 
Here we demonstrated a simple modeling framework that can flexibly fit epidemic data and easily allow additional complexities to account for features during the epidemic. 
We learned two things about fitting epidemic data. 
First, it is important to have ways to incorporate time-varying parameterizations fitting changes in the observed data or behavioural changes. Time-varying transmission rates played an important role in our calibration. 
Second, allowing for testing intensity to disaggregate the transmission rates through time (i.e. does increase in positive cases/confirmation due to an increase in transmission or testing capacity). 
The flexibility of incorporating various time-vary parameterizations allows users to make more realistic features and scenario exploration. 

\clearpage

\bibliographystyle{vancouver}
\bibliography{McMasterReport}

\end{document}