text:
\section{Introduction}

   As of July, 2022, the coronavirus disease 2019 (COVID-19) pandemic, caused by the  severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), has taken more than six million lives worldwide. In response to the pandemic, scientists from all disciplines have made a concerted effort to address the ever growing analytic challenges of prediction, intervention, and control. The mathematical tools employed range from the purely deterministic \ac{ODE}/\ac{PDE} models to describe population dynamics (at the macro or ecological scale) to fully stochastic agent-based models (at the micro scale); from physics-inspired mechanistic models, both stochastic and deterministic, to purely statistical approaches such as ensemble models.  However, despite the longstanding history of the discipline of mathematical epidemiology and the enormous recent efforts, the pandemic has laid bare crucial gaps in the state-of-the-art methodology.

   While the macro models are simple to interpret and easy to calibrate, the micro agent-based models provide more flexibility to model elaborate what-if scenarios. In a similar vein, the mechanistic models provide insights into the underlying biology and epidemiology of the disease but are often outperformed by purely statistical methods when it comes to accuracy in prediction and forecasting. While there is no obvious way to completely bypass such trade-offs between these often diametrically opposed modeling approaches, the \acf{DSA} method  \cite{choi2019modeling,KhudaBukhsh2020InterfaceFocus,DiLauro2022nonMarkov,Wascher2021IDSA,OSU_whitepaper,Harley2022Ebola,KhudaBukhsh2021Prison,Caleb2020Stone,KhudaBukhsh2022Israel}, a survival analytic statistical method derived from dynamical systems, holds some promise. The present paper is about a likelihood-free means of performing \ac{DSA} in the context of non-Markovian models of infectious disease epidemiology. More specifically, we develop a \ac{DSA} method for a non-Markovian epidemic model with vaccination based on the \ac{ABC} framework.

   {Why non-Markov models?} The standard compartmental Markovian models, which employ \acp{CTMC} to keep track of counts of individuals in different compartments (\textit{e.g.}, individuals with different immunological statuses), assume the infectious period and the contact interval \cite{kenah2011biostat,DiLauro2022nonMarkov} are exponentially distributed and are thus characterized by a constant hazard function. This simplistic assumption almost always misrepresents the underlying biology of the disease and hence, is untenable. See \cite{DiLauro2022nonMarkov,vanKampen1998nonMarkov} for a detailed discussion on this point. Also, see \cite[Table 1 and Figure 1]{kenah2011biostat} for a numerical illustration of the bias in the estimates of model parameters if a Markovian model is wrongly assumed when the underlying model is non-Markovian.

   A popular and analytically convenient approach to building more realistic non-Markovian models is to make use of the general theory of  measure-valued processes that keep track of not only the population counts but also additional covariates such as  individual's  age. Here  the word ``age'' is used as an  umbrella term  to refer to the physical age of the individual, the age of infection, time since vaccination etc. While the measure-valued processes do require more mathematical sophistication than their \ac{CTMC} counterparts, as we discuss in \Cref{sec:nonMarkovModel}, the age-stratified population densities can be described by a comparatively simple system of \acp{PDE} in the limit of a large population \cite{DiLauro2022nonMarkov,Franco2021Renewal,Hyman2007Infection_age,Sherborne2018nonMarkov}. The crux of the \ac{DSA} method is to interpret this  mean-field limiting system of \acp{PDE} as describing probability distributions of transfer times (the time required to move from one compartment to another). We shall make this notion clear in \Cref{subsec:DSA} and discuss the statistical benefits of the \ac{DSA} method in \Cref{sec:discussion}.

   As an illustration of the \ac{ABC}-based \ac{DSA} method, we apply it to the COVID-19 epidemic in the state of Ohio, USA. The method is shown to fit the real case count data well and capture nontrivial trends. Detailed numerical results are provided in \Cref{sec:num_results}. In addition to the introduction of \ac{ABC}-\ac{DSA} methodology, our second contribution in this paper is a solution method for the mean-field limiting system of \acp{PDE} with nonlocal boundary conditions. For the sake of algorithmic implementation, we also present it in the form of a pseudocode. An implementation in the Julia programming language is  made available for the wider community.

  The rest of the paper is structured as follows: The stochastic non-Markovian epidemic model is described in \Cref{sec:nonMarkovModel}. The mean-field limit in the form of a system of \acp{PDE} and the \ac{DSA} methodology are also described in \Cref{sec:nonMarkovModel}. In \Cref{sec:PDE_solution}, we describe the main technical details of how we solve the limiting mean-field \ac{PDE} system. We describe the statistical approach to parameter inference in \Cref{sec:prm}, followed by numerical results in \Cref{sec:num_results}. Finally, we conclude with a brief discussion in \Cref{sec:discussion}. For the sake of completeness, additional mathematical details, numerical results, and figures are provided in the Appendix.

  \section{Non-markovian mass-action model}\label{sec:nonMarkovModel}

  Our stochastic model is adapted from \cite{DiLauro2022nonMarkov}. As shown in \Cref{fig:sir}, the model has four compartments: \textbf{S}usceptible, \textbf{V}accinated, \textbf{I}nfectious, and \textbf{R}emoved. Individuals are in exactly one of the four compartments. Upon vaccination, susceptible individuals move to the \textbf{V} compartment. We assume only susceptible individuals are vaccinated. Both (unvaccinated) susceptible and vaccinated individuals can get infected, in which case they move to the \textbf{I} compartment. Finally, infected individuals either recover or are removed. In either case, they move to the terminal compartment \textbf{R}.  In addition to the counts of individuals in the four compartments, we keep track of the age distribution. Here, the word ``age'' refers to the physical age for the susceptible individuals, time since vaccination for the vaccinated individuals, time since infection or age of infection for the infected individuals, and finally, time since recovery or removal for the removed individuals. The age of the removed individuals are of no interest because the removed individuals do not contribute to the dynamics. It is possible to include other important covariates, such as sex, comorbidity, in the model. However, for the sake of simplicity, we only keep track of  age. Suppose we have $n$ individuals.

      \begin{figure}
        \centering
        \begin{tikzpicture}[scale=1.0, every node/.style={transform shape}]
        \begin{scope}
            \node(box1)[draw,rectangle, minimum size=1cm,line width=1pt] at (0,0){Susceptible (\textbf{S})};
            \node(box2)[draw,rectangle, minimum size=1cm,line width=1pt] at (6,-2){Infectious (\textbf{I})};
            \node(box3)[draw,rectangle, minimum size=1cm, line width=1pt] at (12.0,0){Removed (\textbf{R})};
            \node(box4)[draw,rectangle, minimum size=1cm, line width=1pt] at (6,2){Vaccinated (\textbf{V})};

            \draw  [-latex, line width=1pt](box1) to node [minimum size=0.3cm, xshift=-0.5cm, yshift=-0.30cm] {$\Upsilon(t)$} (box2);
            \draw  [-latex, line width=1pt](box1) to node [minimum size=0.3cm, xshift=-0.5cm, yshift=0.450cm] {$\lambda(t,s)$} (box4);
            \draw  [-latex, line width=1pt](box2) to node [minimum size=0.3cm, xshift=0.0cm, yshift=-0.40cm] {$\gamma(s)$} (box3);
            \draw  [-latex, line width=1pt](box4) to node [minimum size=0.3cm, xshift=0.0cm, yshift=0.40cm] {$\alpha(s)$} (box3);
            \draw  [-latex, line width=1pt](box4) to node [minimum size=0.3cm, xshift=-1.25cm, yshift=0.30cm] {$(1-\alpha(s))\Upsilon(t)$} (box2);

            \end{scope}

        \end{tikzpicture}
        \caption{\label{fig:sir}%
          Diagrammatic representation of the non-Markovian compartmental model. At time $t$, a susceptible individual with age $s$ moves either to a vaccinated compartment   with instantaneous rate $\lambda(t,s)$, or to the infectious compartment with instantaneous rate $\Upsilon(t) \coloneqq \int_0^\infty\beta(s)y_{I}(t,s)\differential{s}$. The vaccinated individuals are also subject to infection. An infected individual at age $s$ can be removed with instantaneous rate $\gamma(s)$. The vaccine efficacy $\alpha$ is assumed age-dependent.
          }
      \end{figure}

  The age distribution of individuals in different compartments is  described in terms of finite, point measure-valued processes whose atoms are individual ages. See \cite{DiLauro2022nonMarkov}. Such an approach has been previously used to model age-stratified \ac{BD} processes \cite{Tran2008limit,Tran2009traits}, spatially stratified populations \cite{Fournier2004microscopic}, delays in \acp{CRN} \cite{KhudaBukhsh2020Delay}. The instantaneous rates of jump are assumed to depend on the individual ages. In particular, a susceptible individual of age $s$ has an instantaneous rate $\lambda(t,s)$ of getting vaccinated at time $t$. Moreover, a susceptible individual is also subject to an infection pressure exerted by the infectious individuals. We denote the hazard function for the probability distribution of the contact interval \cite{DiLauro2022nonMarkov,kenah2011biostat} by $\beta$. Therefore, an infectious individual of age of infection $s$ makes an infectious contact with a susceptible or vaccinated individual at an instantaneous rate $\beta(s)$. The hazard function of the probability law of the infectious period is denoted by $\gamma$. The vaccine efficacy is also assumed to depend on the age of the vaccinated individual (time since vaccination). At age $s$, a vaccinated individual moves directly to the \textbf{R} compartment at rate unity with probability $\alpha(s)$, while with probability $(1-\alpha(s))$, they are subject to the same infection pressure as the susceptible individuals and can get infected at the same instantaneous rate, which is the population sum total, scaled by $1/n$, of the $\beta$'s evaluated at the individual ages of infection.

  In the limit of a large population ($n\to \infty$) and under suitable constraints on the hazard functions $\beta, \gamma, \lambda$, the measure-valued processes, when appropriately scaled by $n^{-1}$, can be shown to converge to deterministic continuous measure-valued functions \cite{DiLauro2022nonMarkov,Tran2008limit,Fournier2004microscopic,Tran2009traits,KhudaBukhsh2020Delay}. The densities (with respect to the Lebesgue measure) of those limiting measure-valued functions can be then described in terms of a system of \acp{PDE}.  Let us define
  \begin{itemize}
      \item $y_{S}(t,s)$: The density of susceptible individuals with age $s$ at time $t$;
      \item $y_{V}(t,s)$: The density of vaccinated individuals with age (of vaccination) or time since vaccination $s$ at time $t$;
      \item $y_{I}(t,s)$: The density of infectious individuals with age (of infection) or time since infection $s$ at time $t$;
      \item $y_{R}(t)$: The proportion of removed individuals.
  \end{itemize}
  These quantities are taken over rectangular domains that share a $t$-axis but, in general, may have different length $s$-axes. In particular, the quantities $\ys$ and $\lambda$ are defined over a common domain $R_S$. The quantities $\yv$ and $\alpha$ are defined over a common domain $R_V$, and the quantities $\yi$, $\gamma$, and $\beta$ are defined over a common domain $R_I$. When these distinctions are not needed, we subsume these s-axes into the common interval $\left[0,\infty\right)$.

  Analogous to \cite{DiLauro2022nonMarkov}, the limiting system can be described as
  \begin{align}
      (\partial_t  + \partial_s)\, y_{S}(t,s) &{} = - \left( \lambda(t,s) + \int_0^\infty \beta(u) y_{I}(t,u) \differential{u}  \right) y_{S}(t,s), \label{eq:ys} \\
      (\partial_t  + \partial_s)\, y_{V}(t,s) &{} = - \left(\alpha(s) + (1-\alpha(s))\int_0^\infty \beta(u) y_{I}(t,u) \differential{u} \right) y_{V}(t,s), \nonumber\\% \label{eq:yv} \\
      (\partial_t  + \partial_s)\, y_{I}(t,s) &{} = - \gamma(s) y_{I}(t,s) \nonumber\\% , \label{eq:yi} \\
      \frac{\differential{}}{\differential{t}} y_{R}(t) &{} = \int_0^\infty \left(\alpha(s)y_V(t,s) + \gamma(s)y_{I}(t,s)  \right)\differential{s} , \label{eq:yr}
  \end{align}
  with initial and boundary conditions
  \begin{align}
      y_S(t,0) &{} =0 \mbox{, for all } t\geq 0, \label{eq:bdys}\\ y_S(0,s) &{}  = f_S(s) ,  \nonumber\\
      y_V(t,0) &{} = \int_0^\infty \lambda(t,s)y_{S}(t,s) \differential{s}\label{eq:bdyv}\\  y_V(0,s) &{} =0 \mbox{, for all } s\geq 0,  \nonumber\\
      y_I(t, 0) &{} = \int_0^\infty y_{S}(t,s) \int_0^\infty \beta(u) y_{I}(t,u) \differential{u} \differential{s} + \int_0^\infty (1-\alpha(s))y_{V}(t,s) \int_0^\infty \beta(u) y_{I}(t,u) \differential{u} \differential{s} , \label{eq:bdyi}\\
      y_I(0,s) &{}  = \rho f_I(s) ,, \nonumber\\
      y_R(0) &{} = 0.
  \end{align}
  We assume there are no vaccinated or removed individuals initially.
  The nonnegative functions $f_{S}$ and $f_{I}$ describe the initial age distributions of the susceptible and the infected individuals,  and satisfy
  \begin{align*}
      \int_0^\infty f_{S}(s) \differential{s}  = 1 , \quad
      \int_0^\infty f_{I}(s) \differential{s}  = 1 .
  \end{align*}
  The parameter $\rho$ is the initial proportion of infected individuals in the population. It is worthwhile to point out that the system is mass conserved, \textit{i.e.},
  \begin{align*}
      \int_0^\infty \left( y_{S}(t, s) + y_{V}(t,s) + y_{I}(t,s) \right) \differential{s}  = 1 + \rho - y_R(t), \quad \text{ for all } t \ge 0 .
  \end{align*}
  A consequence of the above conservation law is that the equation \eqref{eq:yr} is in fact redundant.

    \subsection{Continuity constraints on boundary conditions}\label{subsec:cnstbcs}
    When we seek solutions of \eqref{eq:ys}-\eqref{eq:bdyi} in the space of continuous functions, the explicit and implicit boundary data must be assigned continuously at the origin. This criterion is not satisfied freely but imposes additional constraints on the equation coefficients and the initial data. We obtain these constraints by equating the expressions for the explicit and implicit data at the origin and derive
    \begin{align}
        f_S(0) &= 0, \label{eq:cnstbcs} \\
        \lambda(0,s) & =  0 ,  \text{ for all } s \ge 0,  \nonumber \\
        f_I(0) &= \int^\infty_0 \beta(u)f_I(u)\differential{u}. \nonumber
    \end{align}
    The last equality generally reduces the combined degrees of freedom for $\beta$ and $f_I$ by one. We note also that the above constraint on $\lambda$ could be relaxed to holding over the support of $f_S$ only.

  \subsection{\ac{DSA} perspective of the mean-field \acp{PDE}}\label{subsec:DSA}
  Before describing how we solve the limiting mean-field \ac{PDE} system in the next section, let us briefly discuss how \ac{DSA} interprets the limiting \acp{PDE} as probabilistic quantities.

  The \ac{DSA} method \cite{KhudaBukhsh2020InterfaceFocus,DiLauro2022nonMarkov,Wascher2021IDSA,OSU_whitepaper,Harley2022Ebola,KhudaBukhsh2021Prison,Caleb2020Stone,KhudaBukhsh2022Israel} combines classical dynamical systems theory and survival analysis to interpret the mean-field limits of scaled population counts or densities as characterizing probability laws of transfer times between compartments. In the Markovian model, this interpretation boils down to treating the mean-field \acp{ODE} as satisfying a time inhomogeneous Chapman--Kolmogorov equation for the marginal probability law of a Markov chain on the state space $\{\textbf{S}, \textbf{V}, \textbf{I}, \textbf{R} \}$ describing the time-evolving status of a single individual embedded in an infinitely large population. See \cite[Section 3.4]{DiLauro2022nonMarkov} for the standard \ac{SIR} model example. In the non-Markovian case, we  construct a Markov process on the state space  $\{\textbf{S}, \textbf{V}, \textbf{I}, \textbf{R} \}\times [0, \infty)$ to keep track of the time-evolving status as well as the age information of an individual embedded in an infinitely large population. As such, \ac{DSA} interprets the mean-field limiting \acp{PDE} as describing the transition kernel for the Markov process. The transition kernel could be used to simulate individual trajectories.

  Note that the individual-based Markov process (or chain in the Markovian case) is entirely characterized by the mean-field limiting \acp{PDE} (or \acp{ODE} in the Markovian case) describing population-level densities (or counts).  It is precisely in this sense that \ac{DSA} turns an ecological model into an agent-based model! Moreover, this agent-based description gives us the following probability measures
  \begin{align}
      \mu_S(t, A) = \frac{\int_A y_S(t,s) \differential{s} }{ 1+ \rho}, \quad  \mu_{V}(t, A) = \frac{\int_A y_V(t,s) \differential{s} }{ 1+ \rho}, \quad  \mu_I(t, A) = \frac{\int_A y_I(t,s) \differential{s} }{ 1+ \rho} ,
  \end{align}
  for Borel subsets $A$ of $[0, \infty)$. Here, the quantity $\mu_S(t, A)$ describes the probability of a randomly chosen individual with age in the set $A \subset [0, \infty)$ to be in the \textbf{S} compartment at time $t$. The other probability measures are interpreted in a similar fashion. Based on a random sample of infection, vaccination and/or recovery times, which are allowed to be censored, truncated or even aggregated over time or individuals, the above probability measures (and the transition kernel) can be used to write an individual-based product-form likelihood function, called the \ac{DSA}-likelihood. See \cite[Section 3.3]{DiLauro2022nonMarkov} for an explicit example of the \ac{DSA}-likelihood. The \ac{DSA}-likelihood can be used for likelihood-based approaches to parameter inference, such as \ac{MLE}, or Bayesian methods employing \ac{MCMC} techniques. In this paper, we focus on an alternative likelihood-free approach based on the \ac{ABC} method, which we describe in \Cref{sec:num_results}. In the next section, we describe how we solve the mean-field limiting \ac{PDE} system.

  \section{Solving the limiting \ac{PDE} system}\label{sec:PDE_solution}

  For the remainder of the paper, we will make a simplifying assumption that the initial data $f_S$ and $f_I$ are compactly supported. This assumption holds in almost all practical cases of interest. For example individuals are of bounded age, and infection lasts for at most a bounded length of time. Moreover, by a suitable enlargement of the rectangular domains of $\ys$, $\yv$, and $\yi$, respectively $R_S$, $R_V$, and $R_I$, we may also assume without loss of generality that the solutions stay compactly supported in their domains for the time horizon simulated. This follows from the method of characteristics used in \Cref{subsec:solmthd} and the forcing terms on the right-hand side of \eqref{eq:ys} which determine exponential solutions. For concreteness, let us write $R_S=\left[0,T\right]\times\left[0,L_S\right]$, $R_V=\left[0,T\right]\times\left[0,L_V\right]$, and $R_I=\left[0,T\right]\times\left[0,L_I\right]$.
  \subsection{Solution method}\label{subsec:solmthd}
  The governing equations, \eqref{eq:ys}-\eqref{eq:bdyi}, constitute a quasi-linear \ac{PDE} system of nonlocal conservation laws. Since the associated directions of their derivative operators are constant and do not intersect, the solutions into the domain interior are naturally suited for construction by the method of characteristics, which recasts the first order \acp{PDE} into equivalent flow equations. For more details on this construction, we refer readers to standard references, such as \cite[Ch 3]{Evans10} and \cite[Ch 7]{DiBenedetto10}.

  We note that the implicit and nonlocal boundary conditions, \eqref{eq:bdyv}-\eqref{eq:bdyi}, are not standard. In order to handle the boundary conditions,  we exploit a special feature of this system that remarkably depends on the nonlocal character of the implicit boundary terms. Specifically, we  differentiate the implicit boundary conditions in $t$ and recover a governing flow for them, which only depends on the solution and not on the solution's derivatives. For example, consider differentiating \eqref{eq:bdyv}:
  \begin{align*}
  \partial_t\yv(t,0) &= \int^{L_S}_0\left(\partial_t\lambda\ys + \lambda\partial_t\ys\right)ds = \int^{L_S}_0\left( \partial_t\lambda\ys+\lambda\left(-\partial_s\ys-\lambda\ys -\ys\int^{L_I}_0\beta(u)\yi(t,u)du \right)\right)ds\\
  &= \left.-\left[\lambda\ys\right](t,s)\big|^{L_S}_{s=0}\right. + \int^{L_S}_0\left[\left(\partial_t\lambda+\partial_s\lambda\right)\ys\right](t,s) ds - \int^{L_S}_0\left[\lambda^2\ys\right](t,s)ds \\
  &\;\;\;\;- \left(\int^{L_S}_0\left[\lambda\ys\right](t,s)ds \right)\left(\int^{L_I}_0\left[\beta\yi\right](t,s)ds\right).
  \end{align*}
  In the second equality, we have substituted for $\partial_t\ys$ using the $\ys$ equation of \eqref{eq:ys}, and in the third equality, we performed an integration by parts. Continuing we may use the compactly supported initial data to conclude that the boundary term at $s=L_S$ is vanishing and substitute \eqref{eq:bdys} into the boundary term at $s=0$. For the particular case of $\yv$, we see both these terms are vanishing. A similar argument works for $\partial_t\yi(t,0)$. Below we present the resulting evolution equations for the solution's boundary values. Of course, $\ys(t,0)\equiv 0$ by \eqref{eq:bdys}, so we need only consider the remaining two:
  \begin{align}
  \partial_t\yv(t,0) &= \int^{L_S}_0\ys(\dcl)\lambda  ds - \int^{L_S}_0 \lambda^2\ys ds - \left(\int^{L_S}_0 \lambda\ys ds\right)\left(\int^{L_I}_0 \beta\yi ds\right),\label{eq:bdflow_yv}\\
  \partial_t\yi(t,0) &=
  \left(\int^{L_I}_0\beta\yi ds\right)\left[-\int^{L_S}_0\lambda\ys ds - \left(\int^{L_S}_0\ys ds\right)\left(\int^{L_I}_0\beta\yi ds\right)-\int^{L_V}_0\yv(\dt +\ds)\alpha ds\right.\nonumber\\
		&\left. \;\;\;\;\;\; +\left[\big(1-\alpha(t,0)\big)\yv(t,0)\right] -\int^{L_V}_0 \left(1-\alpha\right)\alpha\yv ds- \left(\int^{L_V}_0\left(1-\alpha\right)^2\yv ds\right)\left(\int^{L_I}_0\beta\yi ds\right)\right]\nonumber\\
		&\;\;+\left(\int^{L_S}_0\ys ds+\int^{L_V}_0(1-\alpha)\yv ds\right)\left(\int^{L_I}_0\yi(\dt+\ds)\beta ds + \beta(0)\yi(t,0) - \int^{L_I}_0\beta\gamma\yi ds  \right).\nonumber
  \end{align}
  Note that in \eqref{eq:bdflow_yv}, it is understood that all integrands are evaluated at $(t,s)$ and integrated in $s$. In our particular application, we also have $\left(\dcl\right)\alpha=\partial_s\alpha$ and similarly for the hazard function $\beta$. We have chosen to present \eqref{eq:bdflow_yv} in generality to highlight the underlying structure of \eqref{eq:ys}-\eqref{eq:bdyi}. Note also, the validity of the equations \eqref{eq:bdflow_yv} depends on the compact support of the initial data. Otherwise, these flow equations would need to contain additional boundary terms at the right-end points of the s-axes.

  Our solution method of \eqref{eq:ys}-\eqref{eq:bdyi} now consists of coupling \eqref{eq:ys} in the interior domain together with the alternative flow formulation \eqref{eq:bdflow_yv} of the implicit boundary conditions at the $\left[s=0\right]$ boundary. The explicit boundary data assigned at $\left[t=0\right]$ in \eqref{eq:bdys}-\eqref{eq:bdyi} becomes this system's initial data. For the benefit of readers mostly interested in applications, we defer a rigorous analysis of this alternative formulation and its equivalence with the original one to Appendix~\ref{app:subsecflow}. Also, in this investigation, we assume equation coefficients in classical function spaces and classical notions of solvability. Future work may investigate  deeper questions of equivalence between the original and flow formulations when equation coefficients and the solution belong to suitably identified Sobolev spaces. For more information on the distinctions between classical solutions and weak solutions in Sobolev spaces, we refer again to \cite{Evans10,DiBenedetto10}. We present a numerical demonstration of this method approximating the implicit boundary conditions, \eqref{eq:bdyv}-\eqref{eq:bdyi}, in \Cref{appsubsec:bdconv}.

  \subsection{Numerical implementation}
  For numerically approximating \eqref{eq:ys} and \eqref{eq:bdflow_yv}, we take a semi-discrete approach where we discretize the space axis, $s$, into a number of uniform intervals whose intersection points are called nodes. The number of nodes determines the mesh resolution. We then approximate each component of the true solution -- $\ys$, $\yv$, and $\yi$ -- by its own nodal basis expansion
  \[
  y(t,s) \approx \sum^n_{i=1} y_i(t)\phi_i(s).
  \]
  The $\left\{\phi_i(s)\right\}_i$ functions are piecewise linear, $C^0$ basis splines determined by a given mesh resolution. Each basis spline is associated with an $s$-axis node, $s_i$, and is defined by taking the value 1 at that node and 0 at all others. At each $t$-level, we thus approximate the true solution with a linear, interpolant spline whose value at $(t,s_i)$ is given by $y_i(t)$. Integrals of the true solution and its derived quantities are approximated by taking integrals of their corresponding linear, interpolant splines by trapezoidal rule.

  The $\left\{y_i(t)\right\}_i$ evolve according to the flows of \eqref{eq:ys} and  \eqref{eq:bdflow_yv}. We initialize at the $t=0$ level by projecting the initial data onto the interpolant basis. This is accomplished by sampling those functions at the nodes. We then use an explicit Euler scheme to extend the solution up to the next time step, $t+\Delta t$ along the boundary axis where the implicit data is prescribed according to the \acp{ODE} of \eqref{eq:bdflow_yv}. This determines the first value $y_1(t+\Delta t)$ at the node $s=0$.  For the remaining nodal values, each of their nodes originate along a characteristic from either the previous $t$-level or within that part of the boundary axis that was previously approximated in computing $y_1$.  We now evolve each of those originating, solution values by the method of characteristics along the \acp{ODE} that correspond to \eqref{eq:ys} by a second explicit Euler scheme to obtain the remaining $\left\{y_i(t+\Delta t)\right\}^n_{i=2}$. By iterating this process, we construct a numerical approximation over the desired time interval for simulation.

  In practice, we use an adaptive Euler step to improve solution accuracy. We specify tolerances atol and rtol, which are respectively the absolute and the relative tolerances set by the user. The solver approximates the solution using a single Euler step of size $\Delta t$ and two consecutive smaller Euler steps of size $\Delta t/2$. The solver then checks if at each node the difference between the two Euler approximations satisfies either the atol or rtol thresholds. It does this similarly for the values of the nonlocal integrals. It accepts the step if all quantities meet this criterion, or else it halves the step size and repeats this process. Note that in evaluating the absolute error, we take the magnitude difference in the Euler solutions and then scale by the length of their $s$-axes. This is done because, in general, $\ys$, $\yv$, and $\yi$ are defined on $s$-axes of different lengths, while probability densities vary inversely with the length of their domain. By this adjustment, we normalize the errors in $\ys$, $\yv$, and $\yi$ to a common scale.

  Below we also provide a pseudocode of the algorithm described in this section.
  \begin{pseudocode} Numerical solver for \eqref{eq:ys} and \eqref{eq:bdflow_yv}
  \begin{enumerate}
      \item Initialize the solution at $t=0$ by storing values of the initial data at the nodes,
      \item Propagate from $y_1(t)$ to $y_1(t+\Delta t)$ by taking an explicit Euler step of size $\Delta t$ along \eqref{eq:bdflow_yv},
      \item Propagate from $\left\{y_i(t)\right\}^n_{i=1}$ and $y_1(t+\Delta t)$ to $\left\{y_i(t+\Delta t)\right\}^n_{i=2}$ by taking Euler steps of size at most $\sqrt{2}\Delta t$ along the characteristics of \eqref{eq:ys},
      \item Repeat Steps 2-3 for two consecutive Euler steps with $\Delta t = \Delta t/2$,
      \item If the y-nodal values and integral quantities computed by the single and two-half Euler steps satisfy either of the absolute and relative error thresholds, atol and rtol, then accept the Step 4 solution at $t+\Delta t$ by storing its nodals and continue with $\Delta t = 2\Delta t$. Otherwise, repeat steps 2-4 with $\Delta t=\Delta t/2$.
  \end{enumerate}
  \end{pseudocode}

  \subsubsection*{Code availability}
  We provide an implementation of the algorithms described in this paper in the Julia programming language along with jupyter notebooks and html at \cite{github}. Additional information on Julia may be found at \cite{Julia2017}.

  \section{Model parameters}\label{sec:prm}

  For the vaccine efficacy, $\alpha$, we begin with a function that linearly increases from 0 to a maximum of $\alpha_{\mathrm{eff}}$ over a period of $\alpha_L$ days. After $\alpha_L$, it was taken to be constant. For the contact interval and infectious period hazard functions, $\beta$ and $\gamma$, we choose Weibull distributions as the Weibull distribution has been shown to be a flexible choice in modelling epidemics \cite{DiLauro2022nonMarkov}. Therefore, we have
    \begin{align*}
      \alpha(s) &{} \coloneqq \frac{\min (s, \alpha_L)}{\alpha_L} \times \alpha_{\mathrm{eff}} , \quad  s\ge 0 , \\
      \beta(s) &{} \coloneqq \frac{\beta_\alpha}{\beta_\theta} \left(\frac{s}{\beta_\theta} \right)^{\beta_\alpha-1} , \quad  s \ge 0 , \\
      \gamma(s) &{} \coloneqq \frac{\gamma_\alpha}{\gamma_\theta} \left(\frac{s}{\gamma_\theta} \right)^{\gamma_\alpha-1} , \quad s \ge 0 .
    \end{align*}

  \noindent Further, we smoothed $\alpha$ by a moving integral average to ensure its derivative was everywhere classical. We took the density for the initial infected population, $f_I$, to be approximately uniform over a period of fourteen days subject to the further requirement that it should be compactly supported and continuous. For this purpose, $f_I$ was constant up to day 13.25 and then linearly decayed to 0 by day 13.75. As with $\alpha$, we  smoothed $f_I$ by an integral moving average.

  We inferred the parameters in \eqref{eq:ys}-\eqref{eq:bdyi} using COVID-19 epidemic data for the US state of Ohio, \cite{ODH21}, during the time period Nov 15, 2020- Jan 15, 2021. We directly estimated the hazard rate for vaccination for this same time period from the \ac{ODH} and \ac{CDC} data, \cite{ODH21,CDC21}, using b-splines. We also directly estimated the density for the initial susceptible population, $f_S$, consistent with US census information on age distributions in Ohio, \cite{USCEN21}.
  \subsection{Empirical estimates}\label{subsec:prm_data}
  \subsubsection{The distribution for the initial susceptible population, $f_S$}
  \Cref{fig:fs} shows an empirical density for age distributions in the state of Ohio partitioned into ten year age groups. This data suggested a simple piecewise linear characterization of $f_S$, which is also shown. Additionally, we normalize $f_S$ so that it integrates to one.
    \begin{figure}[h!]
      \centering
      \includegraphics[scale=1]{figs/fs.pdf}
      \caption{\label{fig:fs}%
      An empirical density for the age distributions in Ohio partitioned into ten year brackets. From this we extracted a probability density curve, $f_S$. Later we additionally make this density decay to 0 at the origin consistent with \eqref{eq:cnstbcs}.
      }
  \end{figure}
  \subsubsection{The hazard rate for time to vaccination $\lambda$}
  \Cref{fig:mst_vax} shows the empirical, cumulative vaccination doses given in Ohio from Nov 15th, 2020 - Jan 15th, 2021 (top panel) broken down by age group. These age-breakdowns were estimated from the vaccination counts provided by the \ac{ODH} and \ac{CDC}, \cite{ODH21, CDC21}. We observe that vaccine rollout did not begin until Dec 15th. It also shows the derived empirical estimations of the hazard rate for vaccination, $\lambda$, from this data (bottom two panels).  In order to do this, we measured an empirical survival function for vaccination using the cumulative vaccine doses curves and the age demography in Ohio and then applied to these values the negative $\ln$-transformation. These data points were then fit by least squares with $C^2$ cubic splines, and their pointwise derivatives corresponded to the desired vaccination hazard rate. For more information on computing with splines, we refer to \cite{Schumaker15}. These several curves, parameterized in $t$, were combined into a single $\lambda(t,s)$ function by using smoothed linear transitions in the age-variable $s$ across age-groups.
  \begin{figure}
      \centering
      \includegraphics[scale=0.9]{figs/mst_vax.pdf}
      \caption{\label{fig:mst_vax}%
      \textit{(Top)} Cumulative vaccine doses administered in Ohio. \textit{(Bottom left)} From the cumulative vaccine doses and the Ohio population's age distribution, we log transformed the empirical vaccine survival function. The empirical data is shown as dots. These points were then fit by least-squares cubic splines joined together with $C^2$ smoothness which are also shown. The slopes of these curves represent the vaccine hazard rate. \textit{(Bottom right)} The nonparametric hazard rates, computed as the derivatives of the spline curves in the adjacent panel, are shown.}

  \end{figure}
  \subsection{Estimation by \acl{ABC}}\label{subsec:prm_abc}
  The remaining model parameters were not explicitly given by data, and we used instead an \ac{ABC} scheme to estimate their values. The goal of Bayesian inference is to obtain samples from the posterior distribution of the parameters, which is proportional to a prior distribution on the parameters times the likelihood of the data given the parameters. However, in many cases, it is difficult or expensive to compute or even sample from the posterior. In particular, many methods for sampling from the posterior require computing the likelihood, which may itself be intractable.

  \ac{ABC} provides a method to approximate posterior samples without computing the likelihood. The simplest way to perform \ac{ABC} is the rejection sampling method, which goes as follows: First, one proposes a vector of parameters from the joint prior. Next, one simulates data (or some summary statistic thereof) from the model using this vector of parameters. Finally, one compares this simulated data to the observed data using some distance metric and either accepts the proposed vector of parameters as a sample from an approximation of the posterior or rejects it depending on how close it is to the observed data. What is sufficiently close for acceptance may be determined by an absolute threshold or by retaining some proportion of the best samples. The intuition is that a vector of parameters with higher likelihood will more often produce simulated data that is close to the empirical data, and so \ac{ABC} approximates the usual acceptance procedure based on the likelihood ratio that is often used in MCMC. For a more thorough discussion of \ac{ABC}, we direct the reader to \cite{kypraios2017tutorial,ABC}.

  We proposed the vector of parameters $(\beta_{\alpha}, \gamma_{\alpha}, \gamma_{\theta}, \rho, \alpha_{\textrm{eff}})$ by drawing each quantity independently from a uniform prior with the bounds given in Table \ref{tab:prms}. We then accepted it if the associated infectious period was less than three weeks with 99.9\% probability and if the mean contact interval was less than the mean infectious period. The parameter $\alpha_L$ was fixed to be two weeks. The parameter $\beta_{\theta}$ was not proposed, and hence, does not have a prior, since its value was determined from the other parameters. This followed from the last constraint in \eqref{eq:cnstbcs}. In fact, for the particular pairing of a Weibull distribution for the contact interval and a uniform distribution for the initial infected, it follows from a little algebra that \eqref{eq:cnstbcs} implies $\beta_\theta$ is the length of the support of $f_I$ and independent of the other parameters.

  We then used the \ac{PDE} model to generate a predicted incidence trajectory for the Nov 15, 2020 - Jan 15, 2021 time period and compared the predicted trajectory to the empirical trajectory using the root mean square error (\ac{RMSE}). Note that the prevalence at time $t$ is given by $n \int_0^\infty y_{I}(t,s) \differential{s}$, where $n$ is the total population size.  We generated 5000 \ac{ABC} sample trajectories in this fashion and retained the $10\%$ of sample trajectories with the lowest \ac{RMSE}. The results are shown in  \Cref{fig:abcvar}.

  \section{Numerical results}\label{sec:num_results}

  Here we demonstrate how the \ac{PDE}-\ac{DSA} model can be used in conjunction with with the \ac{ABC} method to infer model parameters from epidemic data. As discussed in \Cref{subsec:prm_abc}, the \ac{ABC} method does not require an explicit likelihood but instead uses a computable error between synthetic and true data values to assess the quality of proposed parameter values. Aggregate population-level counts of infection are prototypical data that would be used with \ac{DSA}, and so we apply this approach to daily reported incidence supplied by the ODH during the period of  Nov 15, 2020, - Jan 15, 2021. This period encompassed the epidemic wave promoted by the rise of the COVID-19 $\alpha$-variant \cite{CDC22} as well as the start of vaccine roll-out in Ohio.

  The primary parameters of interest were for the contact interval characterized by $\beta$ and the infectious period characterized by $\gamma$, along with estimates of the initial amount of infection $\rho$. The vaccine efficacy parameter $\alpha_{\mathrm{eff}}$, while also relevant, would not be identifiable due to still low rate of vaccine administration  during  the time window  of interest. \Cref{fig:abcvar} shows the model predictions  together with the \ac{ODH} reported daily incidence.  The best fit obtained across all 5000 \ac{ABC} samples captures the nontrivial empirical trends.  \Cref{tab:prms} lists the model parameters with their best fit values and the \ac{ABC} posterior credible intervals. The posterior distributions correspond to the parameter values that produced trajectories within the shaded bands of \Cref{fig:abcvar}. In \Cref{fig:mst_meanR0}, we see that the mean time before transmission was approximately estimated to be between $\left(12.25,13.0\right)$ days while the mean time to recovery was approximately estimated to be between $\left(12.25,15.0\right)$ days. These estimates are in line with estimates from  similar studies reported  in the literature  \cite{DiLauro2022nonMarkov}.

\begin{figure}[h!]
      \centering
      \includegraphics{figs/abcvar.pdf}
      \caption{The \ac{PDE} model fit to the \ac{ODH} daily incidence data as calibrated by \ac{ABC}. A seven day moving average of the \ac{ODH} data is shown as a dashed orange line. The solid blue line corresponds to the model's best fit observed across 5000 \ac{ABC} samples, while the blue band represents the best 10\% of all \ac{ABC} sampled trajectories as measured by the \ac{RMSE} between the model prediction and empirical incidence.}
      \label{fig:abcvar}
\end{figure}

\noindent In the second panel of \Cref{fig:mst_meanR0}, we give the posterior distribution for the basic reproduction number $R_0$. This quantity is computed using the formula $R_0=\int^\infty_0 S_{\gamma}(s)\beta(s)\differential{s}$ as in \cite{KhudaBukhsh2020InterfaceFocus}, where $S_\gamma(s)$ is the survival function of the probability distribution characterized by the hazard rate $\gamma$. Note that this formula ignores vaccination, and therefore, is an upper bound on the true basic reproduction number. The model estimated a mean of $R_0=1.28$ with a 95\% credible interval of $\left(0.89,1.94\right)$, which is consistent with other studies of basic reproduction numbers associated with the COVID-19.
  \begin{table}[h!]
      \centering
      \begin{tabular}{|c|c|c|c|c|c|}
           \hline
           Parameter & Unit & Description &
           \begin{tabular}{@{}c@{}}
                Value/  \\
                Prior
           \end{tabular} & \begin{tabular}{@{}c@{}}Best \\Fit\end{tabular} &
           \begin{tabular}{@{}c@{}}
                \ac{ABC} Posterior \\
                95\% Credible Interval
           \end{tabular}\\
           \hline
           \multicolumn{6}{|c|}{\textbf{Contact interval}}\\
           \hline
           $\beta_\alpha$ & -- & Weibull shape parameter & $\left(5,20\right)^*$ & 5.24 & $\left(5.06,9.55\right)$\\
           \hline
           $\beta_\theta$ & day & Weibull scale parameter & -- & 13.5038 & $\left( 13.5036,13.5076\right)$\\
           \hline
           \multicolumn{6}{|c|}{\textbf{Infectious period}}
           \\
           \hline
           $\gamma_\alpha$ & -- & Weibull shape parameter & $\left(2,8\right)^\ast$ & 6.05 & $\left(4.98,7.94\right)$\\
           \hline
           $\gamma_\theta$ & day & Weibull scale parameter & $\left(1,18\right)^*$ & 13.6 & $\left(13.43,15.03\right)$\\
           \hline
           \multicolumn{6}{|c|}{\textbf{Starting infection}}
           \\
           \hline
           $\rho$ & \% &
           \begin{tabular}{@{}c@{}}
                mass of the initial \\
                infected population \\
                relative the size of the \\
                susceptible population
           \end{tabular}
           & $\left(0.1,5\right)$ & 4.48 & $\left(0.79,4.95\right)$
           \\
           \hline
           \multicolumn{6}{|c|}{\textbf{Vaccine}}
           \\
           \hline
           $\alpha_L$ & day & \begin{tabular}{@{}c@{}}
                time for vaccine to achieve  \\
                full efficacy
           \end{tabular} & 14 & -- & --\\
           \hline
           $\alpha_{\mathrm{eff}}$ & \% & vaccine blocking efficacy & $\left(70,100\right)$ & 97.0 & $\left(70.7,99.2\right)$\\
           \hline
      \end{tabular}
      \caption{Model parameters used to fit the \ac{ODH} daily case data from Nov 15, 2020 - Jan 15, 2021. We additionally report their assigned \ac{ABC} priors, best fit values, and \ac{ABC} posterior credible intervals from the numerical simulations presented in \Cref{sec:num_results}. The \ac{ABC} posteriors were obtained by conditioning 5000 \ac{ABC} samples to the best 10\% of observed \ac{RMSE} between the \ac{ODH} daily incidence and model prediction. The priors for $\beta$ and $\gamma$ are listed with * since, in addition to these uniform ranges, we further conditioned on parameters where the infectious period lasted no more than 21 days with 99.9\% probability and where the mean contact interval was less than the mean infectious period. Moreover, a prior for $\beta_\theta$ is not given since \eqref{eq:cnstbcs} determines $\beta_\theta=\beta_\theta(\beta_\alpha,f_I)$. The narrow posterior of $\beta_\theta$ was expected as reasoned in \Cref{subsec:prm_abc}.}
      \label{tab:prms}
  \end{table}
  \begin{figure}[t!]
      \centering
      \includegraphics[scale=0.9]{figs/mst_meanR0.pdf}
      \caption{\textit{(Left)} Posterior distribution of the mean contact interval compared with the posterior distribution for the mean infectious period. \textit{(Right)} Posterior distribution for $R_0$. }
      \label{fig:mst_meanR0}
  \end{figure}

  \section{Discussion}\label{sec:discussion}

The COVID-19 pandemic has spurred the development of a cottage industry of ever-more elaborate mathematical models  of epidemic dynamics that have been applied to   empirical data across  the world with varying degree of success \cite{Koelle22,Driscoll21, Holmdahl20,Jewel20,Barda20}.  Since most of the current and historic  COVID-19 data have been   available in aggregate, the  models tend to focus on  aggregate behavior which may sometimes lead to erroneous insights \cite{klaus2022assortative}.  The \ac{DSA} method  discussed here offers a different modeling approach  and in particular  accounts more fully that a typical compartmental model  for  the heterogeneity of individual  behaviors.

Indeed, even since its introduction in \cite{KhudaBukhsh2020InterfaceFocus} on the eve of an outbreak of the global  COVID-19 pandemic, the \ac{DSA} approach has been shown to be  a viable  way of analyzing epidemic data  in order to predict epidemic progression, evaluate the long term  effects  of public health policies (testing, vaccination, lock downs, etc.) and  individual-level  decisions (masking, social distancing, vaccine hesitancy) \cite{KhudaBukhsh2022Israel, Harley2022Ebola,KhudaBukhsh2021Prison,Wascher2021IDSA,OSU_whitepaper}.
As it was argued in  \cite{KhudaBukhsh2020InterfaceFocus}, the \ac{DSA} method has several advantages over  more traditional  approaches  to  analyzing data in SIR-type epidemic curves. Perhaps the most important one is  that the   method does not require the knowledge of the full curve trajectory or the size of the total susceptible population. Indeed, such quantities are rarely known in practice and often require ad-hoc adjustments leading to severely biased analysis \cite{kenah2011biostat}.

In this work we  discuss a relatively simple  yet powerful non-Markovian extension of the original \ac{DSA} method formulation  introduced in
\cite{KhudaBukhsh2020InterfaceFocus}.  This  extension allows one to take into account the  additional heterogeneity of the transmission patterns due  to both the changes  in infectiousness  over the individuals' infectious periods  and the changes in  immunity over the individuals' periods of vaccine-derived protection.
However, the practical price to pay for these modeling improvements is that a more  elaborate numerical scheme is required to evaluate the \ac{DSA} model along with its increased computational cost to fit empirical data.   Here   we have chosen  to apply a likelihood-free \ac{ABC} approach,  which allows us to avoid the large numerical overhead usually associated with  \ac{DSA} likelihood-based methods \cite{DiLauro2022nonMarkov}. This is accomplished by pregenerating parameter samples and then running simulations in parallel using  modern multi-core capabilities. This capacity of ABC for parallel processing  helps also mitigate the computational burden of a  non-Markovian \ac{DSA}, namely evaluating its nonlocal and implicit  boundary conditions and associated flow formulation of  \eqref{eq:bdflow_yv}. The analytic properties of \eqref{eq:ys} and \eqref{eq:bdflow_yv} are of intrinsic interest in their own right, as  a deeper understanding of solution regularity could help identify numerical schemes of higher and optimal convergence order. We hope to further pursue these lines of investigation in the future. As seen from the numerical examples in \Cref{sec:num_results}, using \ac{ABC} we were able to fit the highly irregular model of Ohio COVID-19 epidemic  with proper accounting  for various sources of uncertainty in all of the relevant  model parameters.

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}
The Atacama Cosmology Telescope ACT is a six-meter telescope located in the Atacama desert intended to measure with high (arcminute) angular resolution the Cosmic Microwave Background (CMB) radiation. Advanced ACTPol, a third generation receiver for ACT is the second upgrade to its cryogenic camera, which fields three two kilo-pixel arrays of time-domain-multiplexed (linear) polarization-sensitive dichroic horn-coupled superconducting Transition Edge Sensors  (TESes)  \cite{henderson_advanced_2016, henderson_highly-multiplexed_2018}. The three arrays in Advanced ACTPol are named PA4, PA5 and PA6 and feature 64, 55 and 55 rows (and 32 columns) of detectors at 150-220, 90-150 and 90-150 GHz respectively.

The Advanced ACTPol array is read out using time domain multiplexing. Each column is read in the time domain by visiting its rows successively. Rows are sampled by the Multi Channel Electronics \cite{battistelli_functional_2008} (MCE), which  digitally (anti-alias) filters and decimates the read out feedback current on each TES to limit the bandwidth of the stored data stream. The user has access to this anti-aliased down-sampled version of the feedback current.

Aliasing arises when a digitally sampled continuous-time signal has
frequency content above the Nyquist frequency ($f_{Ny}$, half of the sampling frequency $f_s$).
Higher frequencies than $f_{Ny}$  fold to lower frequencies and
their power adds to the low frequency content \cite{oppenheim_signals_1997, gallardo_optimizing_2019}. This effect, inherent to the process of sampling adds the out-of-band power to the in-band power spectrum. In the case of a detector system, this out-of-band power adds to the inherent detector noise  increasing the apparent noise power spectral density by a fraction (see section \ref{sec:fom}), which depends on the details of the out-of-band noise. It is common (in the design stage) to add a hardware RL filter to mitigate the out-of-band power in a TES bolometer system, however there is a trade-off between removing frequency content (bandwidth) from the detector by filtering and the ability to bias it (stability). The way the hardware anti-aliasing filter operates with a particular type of sensor often is studied experimentally as the interaction is complicated in practice (it depends on the excess detector noise, readout noise, RL filter fall-off, etc). Figure \ref{fig:average_spectrum} shows the average noise power spectral density for varying bias points in the transition from 40 to $90\% \, \rm {R_n}$ (where $\rm{R_n}$ is the normal resistance of the device) on PA4 read out using only 4 rows, which gives a sampling frequency $16 $ times faster than the nominal science observation mode (see table \ref{tab:Aliasing_mceAquisitionParams} for details) allowing us to probe the high frequency part of the noise in the detection system.

In photon background-limited CMB experiments like Advanced ACTPol, the noise budget of the system is dominated by photon noise; the mapping speed of the overall experiment depends on the total noise budget. Percent level changes in this budget  decrease mapping speed \cite{hill_bolocalc:_2018}. Because aliasing occurs at the end of the detection chain, it is of interest to characterize how it impacts the overall noise budget.

In this work we present measurements of the aliased noise characteristics of the Advanced ACTPol receiver as it was deployed in the field. Measurements are done by re-configuring the software readout electronics to increase the sampling frequency at which detectors are read out at the expense of reading only a fraction of the array. In section \ref{sec:fom} we describe the figure of merit used in this study to measure the aliasing noise contribution. In section \ref{sec:Measurements} we describe the measurement strategy. In section \ref{sec:Results} we discuss our results and in section \ref{sec:Conclusion} we conclude and interpret our results.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{average_HighFreqSpectrum_AR4_detFreq_150.pdf}
    \caption{Average noise spectrum for $150\, \rm{GHz}$ detectors in PA4 acquired by reading out only 4 rows which gives a sampling frequency $16\times$ faster than nominal science observation mode (see table \ref{tab:Aliasing_mceAquisitionParams}). Detectors shown have been flagged as nominally operating using the aliasing fraction cut described in Section \ref{sec:Measurements}.}
    \label{fig:average_spectrum}
\end{figure}

\section{Figure of Merit}
\label{sec:fom}

We quantify aliased noise by evaluating the fraction defined as \begin{equation}
\label{eq:af}
    AliasFraction = \frac{\bar P_{slow} [dac^2/Hz]}{\bar P_{fast}[dac^2/Hz]},
\end{equation}
where $\bar P$ denotes an estimator for the mean power spectral density (in units of $dac^2/Hz$) within the 10 to 60 Hz band at the nominal observing sampling frequency (denoted slow) and a higher sampling frequency (denoted fast). In this study we use the median as an estimator for the mean for being a robust estimator in the presence of line contamination. The sampling frequencies (before decimation) that are accessed are $7.8\, \rm{kHz}$ ($9.1\,\rm{kHz}$) and $125 \, \rm {kHz}$ ($125\, \rm{kHz}$) for PA4 (PA5 and PA6).

This figure of merit has the benefit of being agnostic to modeling systematic errors (as it is purely observational) that can arise from modeling the detector and readout noise in the system.

\section{Measurements}
\label{sec:Measurements}

We reconfigure the Advanced ACTPol readout system \cite{henderson_readout_2016} (MCE) to sample at a faster than nominal science (CMB) observation sampling speeds. We achieve this by varying the number of rows the MCE visits in each column. By having less rows to read, the multiplexer spends less time reading out one column and therefore the overall sampling rate increases. The sampling rate in the MCE is given by
\begin{equation}
    \label{eq:fs}
    f_{s} = \frac{50\, \rm{MHz}}{num\_rows \times row\_len},
\end{equation}
where $num\_rows$ is the number of rows to be read and $row\_len$ is the number of $50\,\rm MHz$ cycles spent per row.
This digitally sampled signal is then digitally filtered via a Butterworth  filter  and a time-stream with sampling frequency
\begin{equation}
    \label{eq:fs-readout}
    f_{s\_readout} = \frac{50\,\rm{MHz}}{num\_rows \times row\_len \times data\_rate}
\end{equation}
is given to the user, where $data\_rate$ is the decimation factor at which data are subsampled. In nominal science observation mode, the (decimated) sampling frequency returned to the user is 300.5 Hz for PA4 (64 rows) and 395 Hz for PA5 and PA6 (55 rows). In this study we read out down to 4 rows which allows to obtain decimated sampling frequencies as high as 4.8kHz for PA4 and 5.4kHz for PA5 and PA6.

Data is acquired with reflective covers on the receiver window entrance. The MCE is configured to acquire data  at different sampling frequencies with the parameter $num\_{rows}$ varying between 4 and the total number of rows in the array. Table \ref{tab:Aliasing_mceAquisitionParams} shows the parameters used in the MCE and their corresponding sampling frequencies as well as the Butterworth filter frequencies at 0.99 gain. We acquire one minute long time streams in each acquisition.
\begin{table}
    \centering
    \begin{tabular}{c|c|c|c|c|c|c}
        array & num\_rows & row\_len & data\_rate & $f_s$[kHz] & $f_{s\_readout}$[Hz] & $f_{c}[Hz]$\\
        \hline \hline
         PA4      & 64, 4 & 100 & 26 & 7.8, 125 &  300.5, 4807 &  64.4, 1031.5\\
         PA5, PA6 & 55, 4 & 100 & 23 & 9.1, 125 &  395.3, 5435 & 75.0, 1031.5
    \end{tabular}
    \caption{MCE parameters used to vary the sampling frequency. Sampling frequencies obtained (defined according to eqs. \ref{eq:fs} and \ref{eq:fs-readout}) are also shown. Digital anti-alias Butterworth filter cutoff frequency at $99\%$ gain is shown as $f_c$.}
    \label{tab:Aliasing_mceAquisitionParams}
\end{table}
Data was visually inspected to discard time intervals where glitches are present. We find two kinds of glitches in the data: discontinuities in the data and short pulse glitches. Discontinuities are found to occur when the data acquisition is started soon after one IV curve is acquired. We remedy this by waiting five minutes after each IV curve to allow the array to reach thermal equilibrium. The dataset presented here showed no identifiable discontinuities. We also find fast glitches, which are jumps in the readout that last one or a few samples in length. We interpret these jumps as cosmic ray hits and remove them from the data by selecting the longest glitch free data in each time stream on a detector to detector basis. One or two glitches (per array) in a one minute acquisition interval are typical in the data stream.

For each glitch-free section of data we compute the power spectral density (periodogram) via the Welch method \cite{scipy_2001}. We compute the median power spectral density in a band that goes from 10 to $60 \, \rm{Hz}$ in each acquisition to avoid $1/f$ noise on the low frequency end and the digital filter roll-off in the high frequency end, while having a large bandwidth to decrease the uncertainty in the mean power density estimation. After computing the median power spectral density we use the aliasing fraction defined in equation \ref{eq:af} for (decimated) sampling rates of  $300.5 \, \rm Hz$ for PA4 and $395 \, \rm Hz$ for PA5 and PA6 (slow) and  $4.8\, \rm kHz$ for PA4 and $5.4\,\rm kHz$ for PA5 and PA6 (fast).
We estimate the error in the median power spectral density as $\sigma_{mean} = \sigma / \sqrt{N_{bins}}$ where $\sigma$ is the standard deviation taken across frequency bins and $N_{bins}$ is the number of frequency samples in the band. This error in the estimation of the mean power spectral density across the band is $<4\%$ for a single observation, which propagates to a $< 6\%$ uncertainty in the aliasing fraction estimation for each detector individually. Uncertainties can be reduced further assuming detectors of the same frequency (90, 150, $220 \, \rm{GHz}$) have similar noise properties (identically distributed aliasing fractions) and aggregating their statistics.

To discard non-working detectors, we only consider detectors  where the aliasing fraction is higher than 0.5 and lower than 1.5. Aliasing fractions outside this range have a very low probability of belonging to a working detector given the observed distribution as suggested by Monte Carlo simulations \cite{gallardo_optimizing_2019}.

\section{Results}
\label{sec:Results}
\begin{figure}
    \centering
    \begin{minipage}{0.49\textwidth}
        \includegraphics[width=\textwidth]{AR4_AliasingFractionHistogram_pctRn_50_nrows64.pdf}
    \end{minipage}
    \begin{minipage}{0.49\textwidth}
        \includegraphics[width=\textwidth]{AR5_AliasingFractionHistogram_pctRn_50_nrows55.pdf}
    \end{minipage}

    \begin{minipage}{0.5\textwidth}
        \includegraphics[width=\textwidth]{AR6_AliasingFractionHistogram_pctRn_50_nrows55.pdf}
    \end{minipage}
    \caption{Histograms showing the measured aliasing fraction from three arrays split by detector frequency. Bias point is $50\, \rm{\%R_n}$ and aliasing fraction compares the nominal sampling frequency to the sampling frequency obtained when multiplexing the first 4 rows (see table \ref{tab:Aliasing_mceAquisitionParams}).}
    \label{fig:alising_fraction_histograms}
\end{figure}

\begin{figure}
    \centering
   \includegraphics[width=0.8\textwidth]{AliasingFraction_perArray.pdf}
    \caption{Alias fraction as a function of bias point. Nominal observations target $50\, \rm {\%  R_n}$ (gray vertical line). $150\, \rm{GHz}$ detectors show an increasing aliasing fraction with increasing bias point $R_n$, while $90$ and $220\,\rm{GHz}$ detectors show a flat aliasing fraction from 40 to $70\, \rm{ \%R_n}$.}
    \label{fig:aliasFrac_vs_pctRn}
\end{figure}

After applying the aliasing fraction cut discussed in Section \ref{sec:Measurements} we make histograms of the aliasing fractions for the first four rows of detectors splitting them by detector frequency. Histograms are consistent with the expected distribution from Monte Carlo simulations which use the observed mean power spectral density distribution (discussed in Section \ref{sec:Measurements}) to predict the aliasing fraction distribution. At $50\% \, \rm{R_n}$ we obtain aliasing fractions over unity of 4, 3 and 5$\%$ for the $150\, \rm GHz$ detectors in PA4, PA5 and PA6 respectively. PA4 showed an over unity aliasing fraction of $10 \, \rm \%$ at $220 \, \rm{GHz}$. PA5 and PA6 showed 7 and 10$\, \rm \%$ respectively at 90GHz. Uncertainties are in the 1-2\% level at $1\sigma$. Figure \ref{fig:alising_fraction_histograms} shows aliasing fraction histograms for the three arrays in Advanced ACTPol split by detector frequency at $50\% \,\rm{R_n}$. Labels indicate the detector frequency, the number of detectors that survived the cut and the estimate of the mean of the distribution. Uncertainties are computed as $\sigma/\sqrt{N}$, where $\sigma$ is the aliasing fraction standard deviation across detectors and $N$ is the number of detectors that survived the cut for each detector type.

Figure \ref{fig:aliasFrac_vs_pctRn} shows the mean aliasing fraction estimate as a function of the bias point expressed as a percentage of the normal resistance for the three arrays. Note that the $150\, \rm{GHz}$ detectors show increasing aliasing fractions with increasing point in the transition ($\rm{\%R_n}$). Detectors operating at  $220\, \rm GHz$  and $90 \, \rm GHz$ show a flat response with respect to bias point between 40 and $70 \% \, \rm {R_n}$. Observed trends support the use of a $50\%$ bias point from aliasing fraction estimates, time constants and stability criteria.
\section{Conclusion}
\label{sec:Conclusion}

We present detector noise aliasing measurements to characterize the performance of time-division SQUID multiplexing in the Advanced ACTPol TES arrays. We demonstrate how these measurements can be done in a MCE system by adjusting the readout parameters to read out only a fraction of the array, and thereby increase the sampling frequency. Excess over unity aliasing fractions (in $\left[\frac{dac^2/Hz}{dac^2/Hz}\right]$) are lower than $\sim 10\%$ for the three arrays studied at $50 \, \rm {\%R_n}$. The $150\, \rm {GHz}$ detectors showed aliasing fractions lower than $ 5 \pm 2 \%$.  More precise measurements could be made by increasing the measurement time, which provides more independent realizations of the noise in the system. Photon loading during these measurements was higher than typical observing conditions (we estimate an equivalent precipitable water vapor of $2\, \rm{mm}$ for our covers-on tests), which suggests that the typical aliasing levels during observations are $2$ to $4\%$ higher than presented here.

Time-division SQUID multiplexing is one of the most mature readout approaches for TES arrays and has been adopted as part of the reference design for the next generation CMB-S4 project \cite{2019arXiv190704473A}.  While other projects are implementing different TES readout techniques with unprecedented detector counts that will be considered for CMB-S4, the performance of these alternative approaches will need to be characterized, optimized and compared to benchmarks achieved by the current generation of CMB experiments, like those presented here.

\begin{acknowledgements}
The authors thank Max Fankhanel and Rodrigo Quiroga for helping with the site operations during the measurements presented here. Comments on this draft were received from two anonymous reviewers.
\end{acknowledgements}

\pagebreak

\bibliographystyle{ltd_style}