Research article Special Issues

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

  • The Dynamical Survival Analysis (DSA) is a framework for modeling epidemics based on mean field dynamics applied to individual (agent) level history of infection and recovery. Recently, the Dynamical Survival Analysis (DSA) method has been shown to be an effective tool in analyzing complex non-Markovian epidemic processes that are otherwise difficult to handle using standard methods. One of the advantages of Dynamical Survival Analysis (DSA) is its representation of typical epidemic data in a simple although not explicit form that involves solutions of certain differential equations. In this work we describe how a complex non-Markovian Dynamical Survival Analysis (DSA) model may be applied to a specific data set with the help of appropriate numerical and statistical schemes. The ideas are illustrated with a data example of the COVID-19 epidemic in Ohio.

    Citation: Colin Klaus, Matthew Wascher, Wasiur R. KhudaBukhsh, Grzegorz A. Rempała. Likelihood-Free Dynamical Survival Analysis applied to the COVID-19 epidemic in Ohio[J]. Mathematical Biosciences and Engineering, 2023, 20(2): 4103-4127. doi: 10.3934/mbe.2023192

    Related Papers:

    [1] Shiqiang Feng, Dapeng Gao . Existence of traveling wave solutions for a delayed nonlocal dispersal SIR epidemic model with the critical wave speed. Mathematical Biosciences and Engineering, 2021, 18(6): 9357-9380. doi: 10.3934/mbe.2021460
    [2] Pannathon Kreabkhontho, Watchara Teparos, Thitiya Theparod . Potential for eliminating COVID-19 in Thailand through third-dose vaccination: A modeling approach. Mathematical Biosciences and Engineering, 2024, 21(8): 6807-6828. doi: 10.3934/mbe.2024298
    [3] Miniak-Górecka Alicja, Nowakowski Andrzej . Sufficient optimality conditions for a class of epidemic problems with control on the boundary. Mathematical Biosciences and Engineering, 2017, 14(1): 263-275. doi: 10.3934/mbe.2017017
    [4] Hamdy M. Youssef, Najat A. Alghamdi, Magdy A. Ezzat, Alaa A. El-Bary, Ahmed M. Shawky . A new dynamical modeling SEIR with global analysis applied to the real data of spreading COVID-19 in Saudi Arabia. Mathematical Biosciences and Engineering, 2020, 17(6): 7018-7044. doi: 10.3934/mbe.2020362
    [5] Chayu Yang, Jin Wang . Computation of the basic reproduction numbers for reaction-diffusion epidemic models. Mathematical Biosciences and Engineering, 2023, 20(8): 15201-15218. doi: 10.3934/mbe.2023680
    [6] Qihua Huang, Hao Wang . A toxin-mediated size-structured population model: Finite difference approximation and well-posedness. Mathematical Biosciences and Engineering, 2016, 13(4): 697-722. doi: 10.3934/mbe.2016015
    [7] Biplab Dhar, Praveen Kumar Gupta, Mohammad Sajid . Solution of a dynamical memory effect COVID-19 infection system with leaky vaccination efficacy by non-singular kernel fractional derivatives. Mathematical Biosciences and Engineering, 2022, 19(5): 4341-4367. doi: 10.3934/mbe.2022201
    [8] Cheng-Cheng Zhu, Jiang Zhu . Spread trend of COVID-19 epidemic outbreak in China: using exponential attractor method in a spatial heterogeneous SEIQR model. Mathematical Biosciences and Engineering, 2020, 17(4): 3062-3087. doi: 10.3934/mbe.2020174
    [9] Guo Lin, Shuxia Pan, Xiang-Ping Yan . Spreading speeds of epidemic models with nonlocal delays. Mathematical Biosciences and Engineering, 2019, 16(6): 7562-7588. doi: 10.3934/mbe.2019380
    [10] A. Q. Khan, M. Tasneem, M. B. Almatrafi . Discrete-time COVID-19 epidemic model with bifurcation and control. Mathematical Biosciences and Engineering, 2022, 19(2): 1944-1969. doi: 10.3934/mbe.2022092
  • The Dynamical Survival Analysis (DSA) is a framework for modeling epidemics based on mean field dynamics applied to individual (agent) level history of infection and recovery. Recently, the Dynamical Survival Analysis (DSA) method has been shown to be an effective tool in analyzing complex non-Markovian epidemic processes that are otherwise difficult to handle using standard methods. One of the advantages of Dynamical Survival Analysis (DSA) is its representation of typical epidemic data in a simple although not explicit form that involves solutions of certain differential equations. In this work we describe how a complex non-Markovian Dynamical Survival Analysis (DSA) model may be applied to a specific data set with the help of appropriate numerical and statistical schemes. The ideas are illustrated with a data example of the COVID-19 epidemic in Ohio.



    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 β and γ 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 βθ is not given since (2.7) determines βθ=βθ(βα,fI). The narrow posterior of βθ was expected as reasoned in Section 4.2.
    Parameter Unit Description Value/Prior Best Fit ABC Posterior 95% Credible Interval
    Contact interval
    βα - Weibull shape parameter (5,20) 5.24 (5.06,9.55)
    βθ day Weibull scale parameter - 13.5038 (13.5036,13.5076)
    Infectious period
    γα - Weibull shape parameter (2,8) 6.05 (4.98,7.94)
    γθ day Weibull scale parameter (1,18) 13.6 (13.43,15.03)
    Starting infection
    ρ % mass of the initial infected population relative the size of the susceptible population (0.1,5) 4.48 (0.79,4.95)
    Vaccine
    αL day time for vaccine to achieve full efficacy 14 - -
    αeff % vaccine blocking efficacy (70,100) 97.0 (70.7,99.2)

     | 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 λ(t,s), or to the infectious compartment with instantaneous rate Υ(t):=0β(s)yI(t,s)ds. The vaccinated individuals are also subject to infection. An infected individual at age s can be removed with instantaneous rate γ(s). The vaccine efficacy α is assumed age-dependent.

    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.

    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 λ(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 β. Therefore, an infectious individual of age of infection s makes an infectious contact with a susceptible or vaccinated individual at an instantaneous rate β(s). The hazard function of the probability law of the infectious period is denoted by γ. 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 α(s), while with probability (1α(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 β's evaluated at the individual ages of infection.

    In the limit of a large population (n) and under suitable constraints on the hazard functions β,γ,λ, the measure-valued processes, when appropriately scaled by n1, 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

    yS(t,s): The density of susceptible individuals with age s at time t;

    yV(t,s): The density of vaccinated individuals with age (of vaccination) or time since vaccination s at time t;

    yI(t,s): The density of infectious individuals with age (of infection) or time since infection s at time t;

    yR(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 yS,yV,yI,yR 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 λ are defined over a common domain RS. The quantities yV and α are defined over a common domain RV, and the quantities yI, γ, and β are defined over a common domain RI. When these distinctions are not needed, we subsume these s-axes into the common interval [0,).

    Analogous to [3], the limiting system can be described as

    (t+s)yS(t,s)=(λ(t,s)+0β(u)yI(t,u)du)yS(t,s), (2.1)
    (t+s)yV(t,s)=(α(s)+(1α(s))0β(u)yI(t,u)du)yV(t,s),(t+s)yI(t,s)=γ(s)yI(t,s)ddtyR(t)=0(α(s)yV(t,s)+γ(s)yI(t,s))ds, (2.2)

    with initial and boundary conditions

    yS(t,0)=0, for all t0, (2.3)
    yS(0,s)=fS(s),yV(t,0)=0λ(t,s)yS(t,s)ds (2.4)
    yV(0,s)=0, for all s0,yI(t,0)=0yS(t,s)0β(u)yI(t,u)duds+0(1α(s))yV(t,s)0β(u)yI(t,u)duds, (2.5)
    yI(0,s)=ρfI(s),,yR(0)=0. (2.6)

    We assume there are no vaccinated or removed individuals initially. The nonnegative functions fS and fI describe the initial age distributions of the susceptible and the infected individuals, and satisfy

    0fS(s)ds=1,0fI(s)ds=1.

    The parameter ρ is the initial proportion of infected individuals in the population. It is worthwhile to point out that the system is mass conserved, i.e.,

    0(yS(t,s)+yV(t,s)+yI(t,s))ds=1+ρyR(t), for all t0.

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

    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

    fS(0)=0,λ(0,s)=0, for all s0,fI(0)=0β(u)fI(u)du. (2.7)

    The last equality generally reduces the combined degrees of freedom for β and fI by one. We note also that the above constraint on λ could be relaxed to holding over the support of fS only.

    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 {S,V,I,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 {S,V,I,R}×[0,) 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

    μS(t,A)=AyS(t,s)ds1+ρ,μV(t,A)=AyV(t,s)ds1+ρ,μI(t,A)=AyI(t,s)ds1+ρ, (2.8)

    for Borel subsets A of [0,). Here, the quantity μS(t,A) describes the probability of a randomly chosen individual with age in the set A[0,) 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.

    For the remainder of the paper, we will make a simplifying assumption that the initial data fS and fI 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 RS, RV, and RI, 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 RS=[0,T]×[0,LS], RV=[0,T]×[0,LV], and RI=[0,T]×[0,LI].

    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):

    tyV(t,0)=LS0(tλyS+λtyS)ds=LS0(tλyS+λ(sySλySySLI0β(u)yI(t,u)du))ds=[λyS](t,s)|LSs=0+LS0[(tλ+sλ)yS](t,s)dsLS0[λ2yS](t,s)ds(LS0[λyS](t,s)ds)(LI0[βyI](t,s)ds).

    In the second equality, we have substituted for tyS using the yS 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=LS is vanishing and substitute (2.3) 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 tyI(t,0). Below we present the resulting evolution equations for the solution's boundary values. Of course, yS(t,0)0 by (2.3), so we need only consider the remaining two:

    tyV(t,0)=LS0yS(t+s)λdsLS0λ2ySds(LS0λySds)(LI0βyIds),tyI(t,0)=(LI0βyIds)[LS0λySds(LS0ySds)(LI0βyIds)LV0yV(t+s)αds+[(1α(t,0))yV(t,0)]LV0(1α)αyVds(LV0(1α)2yVds)(LI0βyIds)]+(LS0ySds+LV0(1α)yVds)(LI0yI(t+s)βds+β(0)yI(t,0)LI0βγyIds). (3.1)

    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 (t+s)α=sα and similarly for the hazard function β. 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 [s=0] boundary. The explicit boundary data assigned at [t=0] 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.

    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 – yS, yV, and yI – by its own nodal basis expansion

    y(t,s)ni=1yi(t)ϕi(s).

    The {ϕi(s)}i functions are piecewise linear, C0 basis splines determined by a given mesh resolution. Each basis spline is associated with an s-axis node, si, 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,si) is given by yi(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 {yi(t)}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+Δt along the boundary axis where the implicit data is prescribed according to the ODEs of (3.1). This determines the first value y1(t+Δ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 y1. 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 {yi(t+Δt)}ni=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 Δt and two consecutive smaller Euler steps of size Δ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.

    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 y1(t) to y1(t+Δt) by taking an explicit Euler step of size Δt along (3.1),

    3. Propagate from {yi(t)}ni=1 and y1(t+Δt) to {yi(t+Δt)}ni=2 by taking Euler steps of size at most 2Δt along the characteristics of (2.1),

    4. Repeat Steps 2-3 for two consecutive Euler steps with Δt=Δ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+Δt by storing its nodals and continue with Δt=2Δt. Otherwise, repeat steps 2-4 with Δt=Δ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].

    For the vaccine efficacy, α, we begin with a function that linearly increases from 0 to a maximum of αeff over a period of αL days. After αL, it was taken to be constant. For the contact interval and infectious period hazard functions, β and γ, we choose Weibull distributions as the Weibull distribution has been shown to be a flexible choice in modelling epidemics [3]. Therefore, we have

    α(s):=min(s,αL)αL×αeff,s0,β(s):=βαβθ(sβθ)βα1,s0,γ(s):=γαγθ(sγθ)γα1,s0.

    Further, we smoothed α by a moving integral average to ensure its derivative was everywhere classical. We took the density for the initial infected population, fI, 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, fI was constant up to day 13.25 and then linearly decayed to 0 by day 13.75. As with α, we smoothed fI 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, fS, consistent with US census information on age distributions in Ohio, [28].

    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 fS, which is also shown. Additionally, we normalize fS 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, fS. Later we additionally make this density decay to 0 at the origin consistent with (2.7).

    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, λ, 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 C2 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 λ(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 C2 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.

    The shape of the hazard function λ 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.

    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 (βα,γα,γθ,ρ,α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 αL was fixed to be two weeks. The parameter βθ 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 βθ is the length of the support of fI 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 n0yI(t,s)ds, 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.

    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 α-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 β and the infectious period characterized by γ, along with estimate of the initial amount of infection ρ. The vaccine efficacy parameter α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 (12.25,13.0) days while the mean time to recovery was approximately estimated to be between (12.25,15.0) 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.

    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 α-variant of the SARS-CoV-2 virus in Ohio, the infectious individuals seem to have attained the highest viral load no later than (12.25,13.0) 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=0Sγ(s)β(s)ds as in [2], where Sγ(s) is the survival function of the probability distribution characterized by the hazard rate γ. 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 (0.89,1.94), which is consistent with other studies of reproduction numbers associated with the COVID-19.

    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.

    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

    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.

    The authors declare no conflict of interest.

    In this section, we begin by assuming the equation coefficients α, β, γ, and λ are continuous and that α, β, and λ also have directional derivatives along characteristics which are continuous as well. We then also assume that the initial data fS and fI 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.

    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 [0,Li]. Thus, we initially solve a separate problem than (2.1) where instead of integrals over [0,) 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=(yS,yV,yI) and C(R) stand for the set of continuous functions on the domain R, we now rewrite (2.1) in the form

    (t+s)yS=Λ1[y],(t+s)yV=Λ2[y],(t+s)yI=Λ3[y], (A.1)

    for functionals Λi:3j=1C(Rj)C(Ri). Moreover, for positive quantities C1 and L1, the Λi satisfy boundedness and Lipschitz conditions

    Λi[y]C1(λ,α,γ,LI,y), (A.2)
    Λi[y1]Λi[y2]L1(α,β,LI,y1,y2)y1y2. (A.3)

    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

    tyV(t,0)=Γ2[y],tyI(t,0)=Γ3[y], (A.4)

    for functionals Γi:3j=1C(Rj)C([0,T]). Their boundedness and Lipschitz conditions are

    Γi[y]C2(Dχλ,Dχα,Dχβ,γ,LS,LV,LI,y), (A.5)
    Γi[y1]Γi[y2]L2(Dχλ,Dχα,Dχβ,γ,LS,LV,LI,y1,y2)y1y2. (A.6)

    Above we let Dχf stand for the max supremum norm of both f and (t+s)f.

    Remark 1. Without loss of generality, we may assume C1, C2, L1, and L2 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 C1=C1(data,y). Similar holds for C2, L1, and L2.

    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 τ=min(s,t) and define FT:3j=1C(Rj)3j=1C(Rj) by

    FT[y](t,s)={fS(sτ)+τ0Λ1[y](tτ+χ,sτ+χ)dχ(tτ)0Γ2[y](θ)dθ+τ0Λ2[y](tτ+χ,sτ+χ)dχfI(sτ)+(tτ)0Γ3[y](θ)dθ+τ0Λ3[y](tτ+χ,sτ+χ)dχ. (A.7)

    Then a fixed point of FT is a classical solution of (2.1) and (3.1) with the initial conditions matching the prescribed explicit boundary data.

    Lemma 1. Let Fh:3j=1C(ˉRj)3j=1C(ˉRj) where ˉRi has the same s-axis as Ri but in the t-axis only spans [0,h]. Then M,c,h all strictly positive and c<1 which depend only on the data so that

    1. The restriction Fh:ˉBM(0)ˉBM(0) and so is bounded on a Banach space,

    2. For y1,y2ˉBM(0), Fh(y1)Fh(y2)cy1y2, and so Fh is a contraction mapping.

    Proof: For the first part, begin by picking an M larger than twice the supremum norm of fS and fI which we label f for brevity. From the sup bounds (A.2) and (A.5), we see that on ˉBM(0) the vector functionals Λ,Γ are bounded in terms of this constant M. We may now pick h so that hmax(Λ,Γ)<f/3. Now examining the sums in (A.7), we see these can be no more than 53f<M owing to the integral domains having length O(h).

    For the second part, we may now use that y1,y2 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 (fS,0,fI) and an M strictly larger than the supremum norm of the f's, then over a time span h=h(data,M), the coupled flow equations (2.1) and (3.1) have a unique bounded solution in the space 3j=1C(ˉRj). 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 fV 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 [kh,(k+1)h] along the t-axis.

    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 s and 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 ϕV(t)=LS0λ(t,s)yS(t,s)ds. We aim to show ϕV(t)=yV(t,0). We will also notate Δvf(t,s)=f(t+v1,s+v2)f(t,s), which is the finite difference of the function f along the vector v. We may now compute

    Δhe1ϕV(t)=LS0(Δhe1λ(t,s)yS(t+h,s)+λ(t,s)Δhe1yS(t,s))ds=LS0(Δhe1λ(t,s)yS(t+h,s)+λ(t,s)Δhe1+he2yS(t,s)λ(t,s)Δhe2yS(t+h,s))ds.

    We next apply summation by parts to the last term:

    LS0λ(t,s)Δhe2yS(t+h,s)ds=LS0λ(t,s)yS(t+h,s+h)ds+LS0λ(t,s)yS(t+h,s)ds=LShλ(t,sh)yS(t+h,s)ds+LS0λ(t,s)yS(t+h,s)ds=LShΔhe2λ(t,s)yS(t+h,s)ds+h0λ(t,s)yS(t+h,s)ds.

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

    tϕV(t)=LS0(yS(t,s)(t+s)λ(t,s)+λ(t,s)(t+s)yS)ds+λ(t,0)yS(t,0)=LS0((t+s)λ(t,s)yS(t,s)+λ(t,s)(λ(t,s)LI0β(u)yI(t,u)du)yS(t,s))ds.

    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 yS(t,0)=0 by (2.3). Upon inspection, we see this is exactly the flow equation for the boundary data of yV in (3.1). It now follows from the uniqueness of solutions for ODEs that ϕV(t) and yV(t,0) must coincide. These arguments may be repeated for the yI solution to show its intended implicit boundary data also evolves according to the flow of (3.1).

    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 yV, and the right panel is for yI. (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.

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

    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, μ, and standard deviations, σ, 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.


    [1] L. An, V. R. Grimm, A. Sullivan, B. L. Turner II, N. Malleson, A. Heppenstall, et al., Challenges, tasks, and opportunities in modeling agent-based complex systems, Ecolog. Model., 457 (2021), 109685. https://doi.org/10.1016/j.ecolmodel.2021.109685 doi: 10.1016/j.ecolmodel.2021.109685
    [2] W. R. KhudaBukhsh, B. S. Choi, E. Kenah, G. A. Rempała, Survival dynamical systems: Individual-level survival analysis from population-level epidemic models, Interface Focus, 10 (2020). https://doi.org/10.1098/rsfs.2019.0048
    [3] F. Di Lauro, W. R. KhudaBukhsh, I. Z. Kiss, E. Kenah, M. Jensen, G. A. Rempała, Dynamic survival analysis for non-markovian epidemic models. J. Royal Soc. Interf., 19 (2022), 20220124. https://doi.org/10.1098/rsif.2022.0124
    [4] B. S. Choi, S. Busch, D. Kazadi, B. Ilunga, E. Okitolonda, Y. Dai, et al., Modeling outbreak data: Analysis of a 2012 Ebola virus disease epidemic in DRC, Biomath, 8 (2019). https://doi.org/10.11145/j.biomath.2019.10.037
    [5] H. Vossler, P. Akilimali, Y. H. Pan, W. R. KhudaBukhsh, E. Kenah, G. A. Rempała, Analysis of individual-level epidemic data: Study of 2018-2020 Ebola outbreak in Democratic Republic of the Congo, Sci. Rep., 12 (2022). https://doi.org/10.1038/s41598-022-09564-4
    [6] W. R. KhudaBukhsh, S. K. Khalsa, E. Kenah, G. A. Rempala, J. H. Tien, COVID-19 dynamics in an Ohio prison,, medRxiv, 2021. Available from: https://www.medrxiv.org/content/early/2021/01/15/2021.01.14.21249782
    [7] M. Wascher, P. M. Schnell, W. R. KhudaBukhsh, M. Quam, J. H. Tien, G. A. Rempała, Monitoring sars-cov-2 transmission and prevalence in populations under repeated testing, 2021. Available from: https://www.medrxiv.org/content/10.1101/2021.06.22.21259342v1
    [8] I. Somekh, W. R. KhudaBukhsh, E. D. Root, G. A. Rempała, E. Sim{ o}es, E. Somekh, Quantifying the population-level effect of the COVID-19 mass vaccination campaign in Israel: A modeling study, Open Forum Infect. Diseases, 9 (2022). https://doi.org/10.1093/ofid/ofac087
    [9] E. Kenah, Contact intervals, survival analysis of epidemic data, and estimation of R0, Biostatistics, 12 (2011), 548–566. https://doi.org/10.1093/biostatistics/kxq068 doi: 10.1093/biostatistics/kxq068
    [10] N. G. van Kampen, Remarks on Non-Markov Processes, Brazilian J. Phys., 28 (1998).
    [11] R. Ferrière, V. C. Tran, Stochastic and deterministic models for age-structured populations with genetically variable traits, In CANUM 2008, volume 27 of ESAIM Proc., pages 289–310. EDP Sci., Les Ulis, 2009. https://doi.org/10.1051/proc/2009033
    [12] E. Franco, M. Gyllenberg, O. Diekmann, One dimensional reduction of a renewal equation for a measure-valued function of time describing population dynamics, Acta Appl. Math., 175 (2021), 12. https://doi.org/10.1007/s10440-021-00440-3 doi: 10.1007/s10440-021-00440-3
    [13] J. M. Hyman, J. Li, Infection-age structured epidemic models with behavior change or treatment, J. Biol. Dynam., 1 (2007), 109–131. https://doi.org/10.1080/17513750601040383 doi: 10.1080/17513750601040383
    [14] N. Sherborne, J. C. Miller, K. B. Blyuss, I. Z. Kiss, Mean-field models for non-markovian epidemics on networks, J. Math. Biol., 76 (2018), 755–778. https://doi.org/10.1007/s00285-017-1155-0 doi: 10.1007/s00285-017-1155-0
    [15] V. C. Tran, Large population limit and time behaviour of a stochastic particle model describing an age-structured population, ESAIM. Probab. Stat., 12 (2008), 345–386. https://doi.org/10.1051/ps:2007052 doi: 10.1051/ps:2007052
    [16] N. Fournier, S. Méléard, A microscopic probabilistic description of a locally regulated population and macroscopic approximations, Ann. Appl. Probab., 14 (2004), 1880–1919.
    [17] W. R. KhudaBukhsh, H.-W. Kang, E. Kenah, G. Rempała, Incorporating age and delay into models for biophysical systems, Phys. Biol., 18 (2021), 10. https://doi.org/10.1088/1478-3975/abc2ab doi: 10.1088/1478-3975/abc2ab
    [18] W. R. KhudaBukhsh, C. D Bastian, M. Wascher, C. Klaus, S. Y. Sahai, M. H. Weir, et al., Projecting covid-19 cases and subsequent hospital burden in ohio, medRxiv, 2022. Available from: https://www.medrxiv.org/content/10.1101/2022.07.27.22278117v1.full.pdf+html
    [19] C. D. Bastian, G. A. Rempala, Throwing stones and collecting bones: Looking for poisson-like random measures, Math. Methods Appl. Sci., 43 (2020), 4658–4668. https://doi.org/10.1002/mma.6224 doi: 10.1002/mma.6224
    [20] G. F. Webb, Theory of nonlinear age-dependent population dynamics, volume 89 of Monographs and Textbooks in Pure and Applied Mathematics. Marcel Dekker, Inc., New York, 1985. https://doi.org/10.1007/BF00250793
    [21] M. Iannelli, M. Martcheva, F. A. Milner, Gender-structured population modeling, volume 31 of Frontiers in Applied Mathematics, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2005.
    [22] L. C. Evans, Partial differential equations, volume 19 of Graduate Studies in Mathematics, American Mathematical Society, Providence, RI, second edition, 2010.
    [23] E. DiBenedetto, Partial differential equations, Cornerstones. Birkhäuser Boston, Ltd., Boston, MA, second edition, 2010.
    [24] C Klaus, PDE-DSA github repository, 2022. Available from: https://github.com/klauscj68/PDE-Vax
    [25] J. Bezanson, A. Edelman, S. Karpinski, V. B. Shah, Julia: A fresh approach to numerical computing, SIAM Rev., 59 (2017), 65–98.
    [26] Ohio Department of Health, Ohio Department of Health COVID Dashboard, Available from: https://coronavirus.ohio.gov/wps/portal/gov/covid-19/dashboards/overview
    [27] Centers for Disease Control and Prevention (CDC), US Centers for Disease Control and Prevention: COVID-19 vaccinations in the United States, County, Available from: https://data.cdc.gov/Vaccinations/COVID-19-Vaccinations-in-the-United-States-County/8xkx-amqh
    [28] U. S. Census, County Population Totals 2010-2019, Available from: https://www.census.gov/data/datasets/time-series/demo/popest/2010s-counties-total.html
    [29] L. L. Schumaker, Spline functions: Computational Methods, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2015. https://doi.org/10.1137/1.9781611973907
    [30] T. Kypraios, P. Neal, D. Prangle, A tutorial introduction to bayesian inference for stochastic epidemic models using approximate bayesian computation, {Math. Biosci.}, 287 (2017), 42–53. https://doi.org/10.1016/j.mbs.2016.07.001 doi: 10.1016/j.mbs.2016.07.001
    [31] S..A. Sisson, Y. N. Fan, M. A. Beaumont, Handbook of Approximate Bayesian Computation, CRC Press, Boca Raton, FL, 2020.
    [32] Centers for Disease Control and Prevention (CDC), US Centers for Disease Control and Prevention: SARS-CoV-2 Variant Classifications and Definitions, Available from: https://www.cdc.gov/coronavirus/2019-ncov/variants/variant-classifications.html
    [33] K. Koelle, M. A. Martin, R. Antia, B. Lopman, N. E. Dean, The changing epidemiology of sars-cov-2, Science, 375 (2022), 1116–1121. https://doi.org/10.1126/science.abm4915 doi: 10.1126/science.abm4915
    [34] M. O'Driscoll, G. R. Dos Santos, L. Wang, D. A. T. Cummings, A. S. Azman, J. Paireau, et al., Age-specific mortality and immunity patterns of SARS-CoV-2, Nature, 590 (2021), 140–145. https://doi.org/10.1038/s41586-020-2918-0 doi: 10.1038/s41586-020-2918-0
    [35] I. Holmdahl, C. Buckee, Wrong but useful—what covid-19 epidemiologic models can and cannot tell us, New England J. Med., 383 (2020), 303–305.
    [36] N. P. Jewell, J. A. Lewnard, B. L. Jewell, Predictive mathematical models of the COVID-19 Pandemic: Underlying principles and value of projections, JAMA, 323 (2020), 1893–1894. https://doi.org/10.1001/jama.2020.6585 doi: 10.1001/jama.2020.6585
    [37] N. Barda, D. Riesel, A. Akriv, J. Levy, U. Finkel, G. Yona, et al. Developing a COVID-19 mortality risk prediction model when individual-level data are not available, Nat. Commun., 11 (2020), 4439. https://doi.org/10.1038/s41467-020-18297-9 doi: 10.1038/s41467-020-18297-9
    [38] C. Klaus, M. Wascher, W. R. KhudaBukhsh, J. H. Tien, G. A. Rempała, E. Kenah, Assortative mixing among vaccination groups and biased estimation of reproduction numbers, Lancet Infect. Diseases, 22 (2022), 579–581. https://doi.org/10.1016/S1473-3099(22)00155-4 doi: 10.1016/S1473-3099(22)00155-4
    [39] O. A. Ladyženskaja, V. A. Solonnikov, N. N. Ural'ceva, Linear and quasilinear equations of parabolic type, in Translations of Mathematical Monographs, Vol. 23. American Mathematical Society, Providence, R.I., 1968. Translated from the Russian by S. Smith.
  • This article has been cited by:

    1. Wasiur R. KhudaBukhsh, Sat Kartar Khalsa, Eben Kenah, Gregorz A. Rempała, Joseph H. Tien, COVID-19 dynamics in an Ohio prison, 2023, 11, 2296-2565, 10.3389/fpubh.2023.1087698
  • Reader Comments
  • © 2023 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(2748) PDF downloads(72) Cited by(1)

Figures and Tables

Figures(7)  /  Tables(1)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog