Likelihood-Free Dynamical Survival Analysis Applied to the COVID-19 Epidemic in Ohio

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 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. 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 [5, 19, 6, 36, 18, 35, 21, 2, 32], a survival analytic statistical method derived from dynamical systems, holds some promise. The present paper is about a likelihood-free means of performing DSA in the context of non-Markovian models of infectious disease epidemiology. More specifically, we develop a DSA method for a non-Markovian epidemic model with vaccination based on the Approximate Bayesian Computation (ABC) framework.

Why non-Markov models? The standard compartmental Markovian models, which employ Continuous Time Markov Chains 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 [17, 6] 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 [6, 34] for a detailed discussion on this point. Also, see [17, 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.

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 [6, 13, 15, 30]. The crux of the 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 DSA method in Section 6.

As an illustration of the ABC-based 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 ABC-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 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 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.

Our stochastic model is adapted from [6]. 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 [6]. Such an approach has been previously used to model age-stratified Birth-Death (BD) processes [33, 9], spatially stratified populations [12], delays in Chemical Reaction Networks [20]. 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 [6, 17] 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 [6, 33, 12, 9, 20]. 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.

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 $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 [6], the limiting system can be described as

	$\displaystyle(\partial_{t}+\partial_{s})\,y_{S}(t,s)$	$\displaystyle{}=-\left(\lambda(t,s)+\int_{0}^{\infty}\beta(u)y_{I}(t,u)\mathrm{d}u\right)y_{S}(t,s),$		(1)
	$\displaystyle(\partial_{t}+\partial_{s})\,y_{V}(t,s)$	$\displaystyle{}=-\left(\alpha(s)+(1-\alpha(s))\int_{0}^{\infty}\beta(u)y_{I}(t,u)\mathrm{d}u\right)y_{V}(t,s),$
	$\displaystyle(\partial_{t}+\partial_{s})\,y_{I}(t,s)$	$\displaystyle{}=-\gamma(s)y_{I}(t,s)$
	$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}y_{R}(t)$	$\displaystyle{}=\int_{0}^{\infty}\left(\alpha(s)y_{V}(t,s)+\gamma(s)y_{I}(t,s)\right)\mathrm{d}s,$		(2)

with initial and boundary conditions

	$\displaystyle y_{S}(t,0)$	$\displaystyle{}=0\mbox{, for all }t\geq 0,$		(3)
	$\displaystyle y_{S}(0,s)$	$\displaystyle{}=f_{S}(s),$
	$\displaystyle y_{V}(t,0)$	$\displaystyle{}=\int_{0}^{\infty}\lambda(t,s)y_{S}(t,s)\mathrm{d}s$		(4)
	$\displaystyle y_{V}(0,s)$	$\displaystyle{}=0\mbox{, for all }s\geq 0,$
	$\displaystyle y_{I}(t,0)$	$\displaystyle{}=\int_{0}^{\infty}y_{S}(t,s)\int_{0}^{\infty}\beta(u)y_{I}(t,u)\mathrm{d}u\mathrm{d}s+\int_{0}^{\infty}(1-\alpha(s))y_{V}(t,s)\int_{0}^{\infty}\beta(u)y_{I}(t,u)\mathrm{d}u\mathrm{d}s,$		(5)
	$\displaystyle y_{I}(0,s)$	$\displaystyle{}=\rho f_{I}(s),,$
	$\displaystyle y_{R}(0)$	$\displaystyle{}=0.$		(6)

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

$\displaystyle\int_{0}^{\infty}f_{S}(s)\mathrm{d}s=1,\quad\int_{0}^{\infty}f_{I}(s)\mathrm{d}s=1.$

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

$\displaystyle\int_{0}^{\infty}\left(y_{S}(t,s)+y_{V}(t,s)+y_{I}(t,s)\right)\mathrm{d}s=1+\rho-y_{R}(t),\quad\text{ for all }t\geq 0.$

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

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 (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]$.

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

	$\displaystyle\alpha(s)$	$\displaystyle{}\coloneqq\frac{\min(s,\alpha_{L})}{\alpha_{L}}\times\alpha_{\mathrm{eff}},\quad s\geq 0,$
	$\displaystyle\beta(s)$	$\displaystyle{}\coloneqq\frac{\beta_{\alpha}}{\beta_{\theta}}\left(\frac{s}{\beta_{\theta}}\right)^{\beta_{\alpha}-1},\quad s\geq 0,$
	$\displaystyle\gamma(s)$	$\displaystyle{}\coloneqq\frac{\gamma_{\alpha}}{\gamma_{\theta}}\left(\frac{s}{\gamma_{\theta}}\right)^{\gamma_{\alpha}-1},\quad s\geq 0.$

We inferred the parameters in (1)-(5) using COVID-19 epidemic data for the US state of Ohio, [28], 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 Ohio Department of Health (ODH) and Centers for Disease Control and Prevention (CDC) data, [28, 10], 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, [4].

Here we demonstrate how the PDE-DSA model can be used in conjunction with with the ABC method to infer model parameters from epidemic data. As discussed in Section 4.2, the 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 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 [11] 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. Figure 4 shows the model predictions together with the ODH reported daily incidence. The best fit obtained across all 5000 ABC samples captures the nontrivial empirical trends. Table 1 lists the model parameters with their best fit values and the 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 [6].

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 [24, 27, 14, 16, 1]. 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 [23]. The 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 [19] on the eve of an outbreak of the global COVID-19 pandemic, the 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) [32, 35, 21, 36, 18]. As it was argued in [19], the 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 [17].

In this work we discuss a relatively simple yet powerful non-Markovian extension of the original DSA method formulation introduced in [19]. 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 DSA model along with its increased computational cost to fit empirical data. Here we have chosen to apply a likelihood-free ABC approach, which allows us to avoid the large numerical overhead usually associated with DSA likelihood-based methods [6]. 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 DSA, namely evaluating its nonlocal and implicit boundary conditions and associated flow formulation of (9). The analytic properties of (1) and (9) 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 Section 5, using 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.