Likelihood-Free Dynamical Survival Analysis applied to the COVID-19 epidemic in Ohio

1.   Introduction

As of October, 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 Ordinary Differential Equation (ODE)/Partial Differential Equation (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. We refer the interested readers to this recent paper ^[1] that highlights some of the challenges and opportunities for complex agent-based models. While there is no obvious way to completely bypass such trade-offs between these often diametrically opposed modeling approaches, the Dynamical Survival Analysis (DSA) method ^[2,3], a survival analytic statistical method derived from dynamical systems, holds some promise. The present paper is about a likelihood-free means of performing Dynamical Survival Analysis (DSA) in the context of non-Markovian models of infectious disease epidemiology. More specifically, we: 1) develop a likelihood-free Dynamical Survival Analysis (DSA) method based on the Approximate Bayesian Computation (ABC) framework for a non-Markovian epidemic model accounting for vaccination of susceptible individuals; and 2) propose a computationally efficient numerical scheme for solving the mean-field semi-linear system of PDEs based on a flow formulation. Although we illustrate the usefulness of our results through the analysis of COVID-19 epidemic data from Ohio during the initial vaccination campaign in late 2020 and early 2021, we note that the methodology can be applied to other epidemics. In fact, variations of the Dynamical Survival Analysis (DSA) methodology have been applied to analyze the foot-and-mounth disease outbreak in the United Kingdom ^[3], the Ebola epidemic in 2012 ^[4], or the recent multiple waves of the Ebola epidemic in the Democratic Republic of Congo (DRC) ^[5]. Moreover, the method is quite flexible in that it can work with various types of epidemic data, such as mass-testing data ^[6], repeated testing data of a closed (campus) population ^[7], repeated testing and vaccination data on a large population ^[8].

Why non-Markov models? The standard compartmental Markovian models, which employ Continuous Time Markov Chains (CTMCs) to keep track of counts of individuals in different compartments (e.g., individuals with different immunological statuses), assume the infectious period and the contact interval ^[3,9] are exponentially distributed and are thus characterized by a constant hazard function. Because covariates such as age of infection, time since vaccination can be thought of as proxies for biological features such as viral load, amount of antibody etc, the simplistic assumption of a constant hazard for the contact interval and infectious period distributions almost always misrepresents the underlying biology of the disease and hence, is untenable. See ^[3,10] for a detailed discussion on this point. Also, see [9,Table 1 and Figure 1] 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. The implications of relaxing the constant hazard assumption appear to be important not only for the present problem of modeling epidemics but also more generally, for developing more realistic biological models that are relatively simple but at the same time are to much larger extent capable of representing the biological heterogeneity (see, for instance, the discussions in ^[3,11]).

Table 1.  Model parameters used to fit the ODH daily case data from Nov 15, 2020 - Jan 15, 2021. We additionally report their assigned Approximate Bayesian Computation (ABC) priors, best fit values, and Approximate Bayesian Computation (ABC) posterior credible intervals from the numerical simulations presented in Section 5. The Approximate Bayesian Computation (ABC) posteriors were obtained by conditioning 5000 Approximate Bayesian Computation (ABC) samples to the best 10% of observed RMSE between the 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 (2.7) determines $\beta _\theta = \beta _\theta(\beta _\alpha , f_I)$. The narrow posterior of $\beta _\theta$ was expected as reasoned in Section 4.2.

Parameter	Unit	Description	Value/Prior	Best Fit	ABC Posterior 95% Credible Interval
Contact interval
$\beta _\alpha$	-	Weibull shape parameter	$\left(5, 20\right)^*$	5.24	$\left(5.06, 9.55\right)$
$\beta _\theta$	day	Weibull scale parameter	-	13.5038	$\left(13.5036, 13.5076\right)$
Infectious period
$\gamma _\alpha$	-	Weibull shape parameter	$\left(2, 8\right)^\ast$	6.05	$\left(4.98, 7.94\right)$
$\gamma _\theta$	day	Weibull scale parameter	$\left(1, 18\right)^*$	13.6	$\left(13.43, 15.03\right)$
Starting infection
$\rho$	%	mass of the initial infected population relative the size of the susceptible population	$\left(0.1, 5\right)$	4.48	$\left(0.79, 4.95\right)$
Vaccine
$\alpha _L$	day	time for vaccine to achieve full efficacy	14	-	-
$\alpha _{\mathrm{eff}}$	%	vaccine blocking efficacy	$\left(70, 100\right)$	97.0	$\left(70.7, 99.2\right)$

---

 | Show Table

DownLoad: CSV

Figure 1.  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) := \int_0^\infty\beta (s)y_{I}(t, s){\rm{d}}{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.

DownLoad: Full-Size Img PowerPoint

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 CTMC counterparts, as we discuss in Section 2, the age-stratified population densities can be described by a comparatively simple system of PDEs in the limit of a large population ^[3,12,13,14]. In addition to being more realistic models, the non-Markovian models have been shown (e.g., in ^[3] and elsewhere) to have better predictive ability. The crux of the Dynamical Survival Analysis (DSA) method is to interpret this mean-field limiting system of PDEs as describing probability distributions of transfer times (the time required to move from one compartment to another). We shall make this notion clear in Section 2.2 and discuss the statistical benefits of the Dynamical Survival Analysis (DSA) method in Section 6.

As an illustration of the Approximate Bayesian Computation (ABC)-based Dynamical Survival Analysis (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 Section 5. In addition to the introduction of Approximate Bayesian Computation (ABC)-Dynamical Survival Analysis (DSA) methodology, our second contribution in this paper is a solution method for the mean-field limiting system of PDEs 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 Section 2. The mean-field limit in the form of a system of PDEs and the Dynamical Survival Analysis (DSA) methodology are also described in Section 2. In Section 3, we describe the main technical details of how we solve the limiting mean-field Partial Differential Equation (PDE) system. We describe the statistical approach to parameter inference in Section 4, followed by numerical results in Section 5. Finally, we conclude with a brief discussion in Section 6. For the sake of completeness, additional mathematical details, numerical results, and figures are provided in the Appendix.

2.   Non-markovian mass-action model

Our stochastic model is adapted from ^[3]. As shown in Figure 1, the model has four compartments: Susceptible, Vaccinated, Infectious, and Removed. Individuals are in exactly one of the four compartments. Upon vaccination, susceptible individuals move to the 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 I compartment. Finally, infected individuals either recover or are removed. In either case, they move to the terminal compartment 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.

The age distribution of individuals in different compartments is described in terms of finite, point measure-valued processes whose atoms are individual ages. See ^[3]. Such an approach has been previously used to model age-stratified Birth-Death (BD) processes ^[11,15], spatially stratified populations ^[16], delays in Chemical Reaction Networks (CRNs) ^[17]. 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 ^[3,9] 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 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 ^{[3,11,15,16,17]}. The densities (with respect to the Lebesgue measure) of those limiting measure-valued functions can be then described in terms of a system of PDEs. Let us define

● $y_{S}(t, s)$: The density of susceptible individuals with age $s$ at time $t$;

● $y_{V}(t, s)$: The density of vaccinated individuals with age (of vaccination) or time since vaccination $s$ at time $t$;

● $y_{I}(t, s)$: The density of infectious individuals with age (of infection) or time since infection $s$ at time $t$;

● $y_{R}(t)$: The proportion of removed individuals.

From a practical perspective, it is important to note that the timescales for chronological age compared to infection age are greatly separated. The age of individuals, in our application, practically differentiates individuals according to their decade of life while a COVID-19 infection lasts for a few weeks. Therefore, we only keep track of the most important covariates, such as the age of infection for the infected individuals. The functions $y_{S}, y_{V}, y_{I}, y_{R}$ are taken over rectangular domains that share a $t$-axis but, in general, may have different length $s$-axes. In particular, the quantities $y_{S}$ and $\lambda$ are defined over a common domain $R_S$. The quantities $y_{V}$ and $\alpha$ are defined over a common domain $R_V$, and the quantities $y_{I}$, $\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 ^[3], the limiting system can be described as

with initial and boundary conditions

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

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, i.e.,

A consequence of the above conservation law is that the equation (2.2) is in fact redundant.

2.1.   Continuity constraints on boundary conditions

When we seek solutions of (2.1)-(2.5) 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

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.

2.2.   Dynamical Survival Analysis (DSA) perspective of the mean-field PDEs

Before describing how we solve the limiting mean-field Partial Differential Equation (PDE) system in the next section, let us briefly discuss how Dynamical Survival Analysis (DSA) interprets the limiting PDEs as probabilistic quantities.

The Dynamical Survival Analysis (DSA) method ^{[2,3,5,6,7,8,18,19]} 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 ODEs 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 [3,Section 3.4] for the standard Susceptible-Infected-Recovered (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, Dynamical Survival Analysis (DSA) interprets the mean-field limiting PDEs 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 PDEs (or ODEs in the Markovian case) describing population-level densities (or counts). It is precisely in this sense that Dynamical Survival Analysis (DSA) turns an ecological model into an agent-based model! Moreover, this agent-based description gives us the following probability measures

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 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 Dynamical Survival Analysis (DSA)-likelihood. See [3,Section 3.3] for an explicit example of the Dynamical Survival Analysis (DSA)-likelihood. The Dynamical Survival Analysis (DSA)-likelihood can be used for likelihood-based approaches to parameter inference, such as Maximum Likelihood Estimate (MLE), or Bayesian methods employing Markov Chain Monte Carlo (MCMC) techniques. In this paper, we focus on an alternative likelihood-free approach based on the Approximate Bayesian Computation (ABC) method, which we describe in Section 5. In the next section, we describe how we solve the mean-field limiting Partial Differential Equation (PDE) system.

3.   Solving the limiting Partial Differential Equation (PDE) system

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 $y_{S}$, $y_{V}$, and $y_{I}$, 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 Section 3.1 and the forcing terms on the right-hand side of (2.1) 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]$.

3.1.   Solution method

The governing equations, (2.1)-(2.5), constitute a semilinear Partial Differential Equation (PDE) system of nonlocal conservation laws. Similar equations of McKendrick-von Foerster type have been studied extensively, such as in ^[20] by semigroups and in ^[21] by integral equations. The latter also presents numerical constructions. Here we give an alternative solution method that accommodates our particular set of nonlocal boundary conditions.

Solutions in the domain interior are naturally suited for construction by the method of characteristics. For more details on this construction, we refer readers to standard references, such as [22,Ch 3] and [23,Ch 7]. In order to handle the nonlocal and implicit boundary conditions (2.4)-(2.5), 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 (2.4):

In the second equality, we have substituted for $\partial_ty_{S}$ using the $y_{S}$ equation of (2.1), 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 (2.3) into the boundary term at $s = 0$. For the particular case of $y_{V}$, we see both these terms are vanishing. A similar argument works for $\partial_ty_{I}(t, 0)$. Below we present the resulting evolution equations for the solution's boundary values. Of course, $y_{S}(t, 0)\equiv 0$ by (2.3), so we need only consider the remaining two:

Note that in (3.1), it is understood that all integrands are evaluated at $(t, s)$ and integrated in $s$. In our particular application, we also have $\left(\partial_t+\partial_s\right)\alpha = \partial_s\alpha$ and similarly for the hazard function $\beta$. We have chosen to present (3.1) in generality to highlight the underlying structure of (2.1)-(2.5). Note also, the validity of the equations (3.1) 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 (2.1)-(2.5) now consists of coupling (2.1) in the interior domain together with the alternative flow formulation (3.1) of the implicit boundary conditions at the $\left[s = 0\right]$ boundary. The explicit boundary data assigned at $\left[t = 0\right]$ in (2.3)-(2.5) 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 D.1. 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 ^[22,23]. We present a numerical demonstration of this method approximating the implicit boundary conditions, (2.4)-(2.5), in Appendix A.2.

3.2.   Numerical implementation

For numerically approximating (2.1) and (3.1), 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 – $y_{S}$, $y_{V}$, and $y_{I}$ – by its own nodal basis expansion

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 (2.1) and (3.1). 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 ODEs of (3.1). 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 ODEs that correspond to (2.1) 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, $y_{S}$, $y_{V}$, and $y_{I}$ 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 $y_{S}$, $y_{V}$, and $y_{I}$ to a common scale.

Below we also provide a pseudocode of the algorithm described in this section.

Pseudocode 1. Numerical solver for (2.1) and (3.1)

1. Initialize the solution at $t = 0$ by storing values of the initial data at the nodes,

2. Propagate from $y_1(t)$ to $y_1(t+\Delta t)$ by taking an explicit Euler step of size $\Delta t$ along (3.1),

3. 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 (2.1),

4. Repeat Steps 2-3 for two consecutive Euler steps with $\Delta t = \Delta t/2$,

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

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 ^[24]. Additional information on Julia may be found at ^[25].

4.   Model parameters

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 ^[3]. Therefore, we have

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 (2.1)-(2.5) using COVID-19 epidemic data for the US state of Ohio, ^[26], 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 ODH and Centers for Disease Control and Prevention (CDC) data, ^[26,27], 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, ^[28].

4.1.   Empirical estimates

4.1.1.   The distribution for the initial susceptible population, $f_S$

Figure 2 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.

Figure 2.  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 (2.7).

DownLoad: Full-Size Img PowerPoint

4.1.2.   The hazard rate for time to vaccination $\lambda$

Figure 3 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 ODH and the Centers for Disease Control and Prevention (CDC), ^[26,27]. We observe that vaccine rollout did not begin until Dec 15th. It also shows the derived empirical estimates 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 ^[29]. 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.

Figure 3.  (Top) Cumulative vaccine doses administered in Ohio. (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. (Bottom right) The nonparametric hazard rates, computed as the derivatives of the spline curves in the adjacent panel, are shown.

DownLoad: Full-Size Img PowerPoint

The shape of the hazard function $\lambda$ corresponding to vaccination depends largely on the government decisions to prioritize vaccination of certain vulnerable population ahead of the rest rather than the underlying biology of the disease. Recall that, by assumption, there are no vaccinated individuals initially, which is in agreement with the chosen time window for our numerical example. However, modifications can be made in a straightforward fashion should a pre-vaccinated population needs to be incorporated in the model. The parameter estimation method according to the Dynamical Survival Analysis (DSA) method described below will still work because the Dynamical Survival Analysis (DSA) method works with the notion of an effective population size, rather than the actual population size of a region.

4.2.   Estimation by Approximate Bayesian Computation

The remaining model parameters were not explicitly given by data, and we used instead an Approximate Bayesian Computation (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.

Approximate Bayesian Computation (ABC) provides a method to approximate posterior samples without computing the likelihood. The simplest way to perform Approximate Bayesian Computation (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 Approximate Bayesian Computation (ABC) approximates the usual acceptance procedure based on the likelihood ratio that is often used in MCMC. For a more thorough discussion of Approximate Bayesian Computation (ABC), we direct the reader to ^[30,31].

We proposed the vector of parameters $(\beta _{\alpha }, \gamma _{\alpha }, \gamma _{\theta}, \rho, \alpha _{ \rm{eff}})$ by drawing each quantity independently from a uniform prior with the bounds given in Table 1. 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 (2.7). 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 (2.7) implies $\beta _\theta$ is the length of the support of $f_I$ and independent of the other parameters.

We then used the Partial Differential Equation (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 (RMSE). Note that the prevalence at time $t$ is given by $n \int_0^\infty y_{I}(t, s) {\rm{d}}{s}$, where $n$ is the total population size. We generated 5000 Approximate Bayesian Computation (ABC) sample trajectories in this fashion and retained the $10\%$ of sample trajectories with the lowest RMSE. The results are shown in Figure 4.

Figure 4.  The PDE model fit to the ODH daily incidence data as calibrated by ABC. A seven day moving average of the ODH data is shown as a dashed orange line. The solid blue line corresponds to the model's best fit observed across 5000 ABC samples, while the blue band represents the best 10% of all ABC sampled trajectories as measured by the RMSE between the model prediction and empirical incidence.

DownLoad: Full-Size Img PowerPoint

5.   Numerical results and their biological interpretation

Here we demonstrate how the Partial Differential Equation (PDE)-Dynamical Survival Analysis (DSA) model can be used in conjunction with with the Approximate Bayesian Computation (ABC) method to infer model parameters from epidemic data. As discussed in Section 4.2, the Approximate Bayesian Computation (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 Dynamical Survival Analysis (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 ^[32] as well as the start of vaccine roll-out in Ohio. Therefore, it allowed us to study non-standard dynamics without having to adjust for multiple strains of the virus or the long-term effect of vaccination.

The primary parameters of interest were the contact interval characterized by $\beta$ and the infectious period characterized by $\gamma$, along with estimate 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. Figure 4 shows the model predictions together with the ODH reported daily incidence. The best fit obtained across all 5000 Approximate Bayesian Computation (ABC) samples captures the nontrivial empirical trends. Table 1 lists the model parameters with their best fit values and the Approximate Bayesian Computation (ABC) posterior credible intervals. The posterior distributions correspond to the parameter values that produced trajectories within the shaded bands of Figure 4. In Figure 5, 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 ^[3].

Figure 5.  (Left) Posterior distribution of the mean contact interval compared with the posterior distribution for the mean infectious period. (Right) Posterior distribution for $R$.

DownLoad: Full-Size Img PowerPoint

Note that the contact interval distribution and the infectious periods can be thought of as proxies for the viral load after infection. Therefore, one can gain some qualitative insights into the biology of the disease from the shape of the estimated contact intervals and the infectious periods. For instance, for the $\alpha$-variant of the SARS-CoV-2 virus in Ohio, the infectious individuals seem to have attained the highest viral load no later than $\left(12.25, 13.0\right)$ days after infection. Also, they appear to carry an amount of viral load large enough to infect others (with a significantly high probability) only for a short duration of a few days (when the density of the contact interval distribution is high) after which the viral load remains low (below a threshold) and decays slowly as suggested by the long tail of the infectious period. This description, we believe, is more truthful to the underlying biology of the disease than what we could obtain from the traditional CTMC-based Markov models in which the hazard functions corresponding to the contact intervals and the infectious periods would be assumed constant (regardless of the time since infection, and thereby, viral load).

In the second panel of Figure 5, we give the posterior distribution for the reproduction number $R$. This quantity is computed using the formula $R = \int^\infty_0 S_{\gamma }(s)\beta (s){\rm{d}}{s}$ as in ^[2], 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 reproduction number. The model estimated a mean of $R = 1.28$ with a 95% credible interval of $\left(0.89, 1.94\right)$, which is consistent with other studies of reproduction numbers associated with the COVID-19.

6.   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 ^{[33,34,35,36,37]}. 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 ^[38]. The Dynamical Survival Analysis (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 ^[2] on the eve of an outbreak of the global COVID-19 pandemic, the Dynamical Survival Analysis (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) ^[5,6,7,8,18]. As it was argued in ^[2], the Dynamical Survival Analysis (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 ^[9]. In fact, the Dynamical Survival Analysis (DSA) method can estimate an effective population size based on a discount estimator [3,Section 3.5]. Therefore, not only do we not need the size of the total susceptible population, but also can we ignore the size of the removed population.

In this work we discuss a relatively simple yet powerful non-Markovian extension of the original Dynamical Survival Analysis (DSA) method formulation introduced in ^[2] accounting for vaccination. The idea is to transform the non-Markovian system into a measure-valued Markovian system for which a significant body of literature has been developed over the years. In particular, our work is inspired by works such as ^[11,15,16]. 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 Dynamical Survival Analysis (DSA) model along with its increased computational cost to fit empirical data. Here we have chosen to apply a likelihood-free Approximate Bayesian Computation (ABC) approach, which allows us to avoid the large numerical overhead usually associated with Dynamical Survival Analysis (DSA) likelihood-based methods ^[3]. 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 Dynamical Survival Analysis (DSA), namely evaluating its nonlocal and implicit boundary conditions and associated flow formulation of (3.1).

As seen from the numerical examples in Section 5, using Approximate Bayesian Computation (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.

A.   Acronyms

ABC Approximate Bayesian Computation

BD Birth-Death

CDC Centers for Disease Control and Prevention

CRN Chemical Reaction Network

CTMC Continuous Time Markov Chain

DSA Dynamical Survival Analysis

MCMC Markov Chain Monte Carlo

MLE Maximum Likelihood Estimate

ODE Ordinary Differential Equation

ODH Ohio Department of Health

PDE Partial Differential Equation

RMSE Root Mean Square Error

SIR Susceptible-Infected-Recovered

Acknowledgments

WRK acknowledges financial support from the London Mathematical Society (LMS) under a Scheme 4 grant (Ref No: 42118) and an International Collaboration Fund awarded by the Faculty of Science, University of Nottingham (UoN). GAR acknowledges support from the National Science Foundation under grants DMS-2027001 and DMS-1853587.

Conflict of interest

The authors declare no conflict of interest.

A.   Appendix

A.1.   Analysis of the coupled flow system (2.1) and (3.1)

In this section, we begin by assuming the equation coefficients $\alpha$, $\beta$, $\gamma$, and $\lambda$ are continuous and that $\alpha$, $\beta$, and $\lambda$ also have directional derivatives along characteristics which are continuous as well. We then also assume that the initial data $f_S$ and $f_I$ are continuous, compactly supported, and satisfying the compatibility conditions (2.7). As observed in Remark 2 below, we will also come to require that coefficients are bounded.

A.1.1.   Structure conditions

First we truncate the domains of the integrals in (2.1) to span their respective $R$'s which we now index from 1 to 3. These $R$'s share a common t-axis but have different s-axes respectively spanning $\left[0, L_i\right]$. Thus, we initially solve a separate problem than (2.1) where instead of integrals over $\left[0, \infty\right)$ we have finite domains. However, because of the compactness of our initial data, this problem will be equivalent to the original formulation a posteriori. Letting $y = \left(y_{S}, y_{V}, y_{I}\right)$ and $\mathcal{C}(R)$ stand for the set of continuous functions on the domain $R$, we now rewrite (2.1) in the form

for functionals $\Lambda_i:\prod^3_{j = 1}\mathcal{C}\left(R_j\right)\rightarrow \mathcal{C}\left(R_i\right)$. Moreover, for positive quantities $C_1$ and $\mathcal{L}_1$, the $\Lambda_i$ satisfy boundedness and Lipschitz conditions

The boundary flow kernels in (3.1) are more elaborate but they satisfy similar structure conditions. The key point is that the flow equations are algebraic combinations of functionals which themselves satisfy boundedness and Lipschitz conditions. We may rewrite (3.1) as

for functionals $\Gamma_i:\prod^3_{j = 1}\mathcal{C}\left(R_j\right)\rightarrow \mathcal{C}\left(\left[0, T\right]\right)$. Their boundedness and Lipschitz conditions are

Above we let $\left\|D_\chi f\right\|_\infty$ stand for the max supremum norm of both $f$ and $\left(\partial_t+\partial_s\right)f$.

Remark 1. Without loss of generality, we may assume $C_1$, $C_2$, $\mathcal{L}_1$, and $\mathcal{L}_2$ are increasing functions of their arguments.

Remark 2. The argument dependence of the Lipschitz constants in (A.2)-(A.3) and (A.5)-(A.6) on norms of the equation coefficients quantifies boundedness assumptions needed on equation coefficients for the Lipschitz assumptions to hold.

Remark 3. Hereafter, we will abbreviate a constant's dependence on equation coefficients or geometry as data. So $C_1 = C_1({\rm{data}}, \left\|y\right\|_\infty)$. Similar holds for $C_2$, $\mathcal{L}_1$, and $\mathcal{L}_2$.

A.1.2.   A fixed point characterization of (A.1) and (A.4)

Mimicking the standard approach in ODEs, we now recast the solution of (A.1) and (A.4) as a fixed point for an integral equation. To that end, let $\tau^* = \min(s, t)$ and define $\mathcal{F}_T:\prod^3_{j = 1}\mathcal{C}\left(R_j\right)\rightarrow \prod^3_{j = 1}\mathcal{C}\left(R_j\right)$ by

Then a fixed point of $\mathcal{F}_T$ is a classical solution of (2.1) and (3.1) with the initial conditions matching the prescribed explicit boundary data.

Lemma 1. Let $\mathcal{F}_h:\prod^3_{j = 1}\mathcal{C}\left(\bar{R}_j\right)\rightarrow \prod^3_{j = 1}\mathcal{C}\left(\bar{R}_j\right)$ where $\bar{R}_i$ has the same s-axis as $R_i$ but in the t-axis only spans $\left[0, h\right]$. Then $\exists M, c, h$ all strictly positive and $c < 1$ which depend only on the data so that

1. The restriction $\mathcal{F}_h:\bar{B}_M(0)\rightarrow \bar{B}_M(0)$ and so is bounded on a Banach space,

2. For $y_1, y_2\in\bar{B}_M(0)$, $\left\|\mathcal{F}_h(y_1)-\mathcal{F}_h(y_2)\right\|_\infty\leq c\left\|y_1-y_2\right\|_\infty$, and so $\mathcal{F}_h$ is a contraction mapping.

Proof: For the first part, begin by picking an $M$ larger than twice the supremum norm of $f_S$ and $f_I$ which we label $\left\|f\right\|_\infty$ for brevity. From the sup bounds (A.2) and (A.5), we see that on $\bar{B}_M(0)$ the vector functionals $\Lambda, \Gamma$ are bounded in terms of this constant M. We may now pick $h$ so that $h\max(\left\|\Lambda\right\|_\infty, \left\|\Gamma\right\|_\infty) < \left\|f\right\|_\infty/3$. Now examining the sums in (A.7), we see these can be no more than $\frac{5}{3}\left\|f\right\|_\infty < M$ owing to the integral domains having length $O(h)$.

For the second part, we may now use that $y_1, y_2$ are supremum bound by $M$ to reduce the dependence of the Lipschitz constants in (A.3) and (A.6) to just the data. But now we conclude in the usual manner by seeing that the integral domains in (A.7) have length $O(h)$ and so, by further reduction in $h$, we make the product of $h$ with those Lipschitz constants suitably underneath a threshold $c < 1$.

Corollary 1. For initial data $\left(f_S, 0, f_I\right)$ and an $M$ strictly larger than the supremum norm of the $f$'s, then over a time span $h = h({\rm{data}}, M)$, the coupled flow equations (2.1) and (3.1) have a unique bounded solution in the space $\prod^3_{j = 1}\mathcal{C}\left(\bar{R}_j\right)$. This solution may be further taken to be bounded by M over the time span $h$.

Proof: This follows from the last lemma and the uniqueness of fixed points for contraction mappings in Banach space.

Remark 4. When instead the initial data has a nonzero density $f_V$ for the vaccinated compartment, the same arguments and conclusions hold just with $M$ now also larger than that component's supremum norm.

Remark 5. Since the size of $h$ is in a monotone decreasing relationship with M but otherwise depends on data which are fixed, we see that so long as an a priori global bound can be derived for solutions of (2.1) and (3.1), then they will persist globally and be unique. This follows by the usual iteration of extending the solution by patches $\left[kh, (k+1)h\right]$ along the t-axis.

A.1.3.   On the equivalence of (2.1) and (3.1) with the original formulation (2.1) and (2.3)-(2.5)

It remains to show that fixed points of (A.7) coincide at the t-axis with their intended implicit boundary values (2.3)-(2.5). The technical concern is that while (A.7) is sufficient to imply a fixed point solution is differentiable along characteristics, it does not imply a priori that a solution is differentiable along either the t-direction or s-direction individually. Recall in motivating (3.1), we treated the $\partial_s$ and $\partial_t$ derivatives separately. Nevertheless, we show here fixed point solutions do take the intended implicit boundary data. For further examples on the uses of difference quotients and summation by parts in the theory of PDEs, we refer to standard references such as ^[39].

For the fixed point solution $y$ of (A.7), Let us notate $\phi_V(t) = \int^{L_S}_0 \lambda (t, s)y_{S}(t, s)ds$. We aim to show $\phi_V(t) = y_V(t, 0)$. We will also notate $\Delta_{v}f(t, s) = f(t+v_1, s+v_2)-f(t, s)$, which is the finite difference of the function $f$ along the vector $v$. We may now compute

We next apply summation by parts to the last term:

If we now combine these equations, divide by $h$, and then pass to the limit as $h\rightarrow 0$, we obtain from the regularity of all terms involved that

Notice in the last equality that we have used directional derivatives of the solution along characteristics are classical and given by (2.1). Also, we used that $y_{S}(t, 0) = 0$ by (2.3). Upon inspection, we see this is exactly the flow equation for the boundary data of $y_{V}$ in (3.1). It now follows from the uniqueness of solutions for ODEs that $\phi_V(t)$ and $y_{V}(t, 0)$ must coincide. These arguments may be repeated for the $y_{I}$ solution to show its intended implicit boundary data also evolves according to the flow of (3.1).

A.2.   Demonstrating numerical convergence of the implicit boundary conditions

In Figure 6 we demonstrate the Partial Differential Equation (PDE) numerical solution better satisfying the implicit boundary conditions as the mesh resolution and ode integrator tolerances become increasingly fine. For the figure legend, we note that nnd is the number of mesh nodes along the s-axis while atol and rtol are respectively the adaptive Euler schemes absolute and relative tolerances for refining the step. Recall that, since the magnitude of a probability density varies inversely with the length of its domain, the absolute errors shown are for the solution density multiplied by the size of its domain. In this sense, the absolute error reported is for the solution density if it were first rescaled to a density on the unit interval.

Figure 6.  Convergence demonstrations of the PDE Dynamical Survival Analysis (DSA) solver at the implicit boundary axis as the mesh resolution and ode solver tolerances become increasingly fine when using the best fit parameters obtained by Approximate Bayesian Computation (ABC) in Figure 4 and Table 1 of the main body. (Top) The absolute error, scaled by the domain length, between the solution's value evaluated at the boundary where implicit data is prescribed and the implicit integral quantity which it is supposed to equal. These are respectively the left and right hand sides of (2.4)-(2.5). The left panel is for $y_V$, and the right panel is for $y_I$. (Bottom) Similar as the above panels but for the relative error as measured by the magnitude error divided by the solution value at that boundary point.

DownLoad: Full-Size Img PowerPoint

We observe a slight early increase in the relative error of $y_V$'s implicit boundary data as resolutions become increasingly fine. We expect this is due to the true magnitude of $y_V$ becoming increasingly small, which thus leads to division by a smaller number.

A.3.   Parameter posteriors

Figure 7.  Posterior distributions for ABC estimated model parameters. The best 10% of all 5000 ABC samples were kept. The parameter descriptions are given in Table 1 of the main body; however, we give the means, $\mu$, and standard deviations, $\sigma$, for the contact interval and infectious period rather than their shape and scale parameters. This amounts to an equivalent parameterization of the two-parameter Weibull distributions.

DownLoad: Full-Size Img PowerPoint