Bi-fidelity stochastic collocation methods for epidemic transport models with uncertainties

Uncertainty in data is certainly one of the main problems in epidemiology, as shown by the recent COVID-19 pandemic. The need for efficient methods capable of quantifying uncertainty in the mathematical model is essential in order to produce realistic scenarios of the spread of infection. In this paper, we introduce a bi-fidelity approach to quantify uncertainty in spatially dependent epidemic models. The approach is based on evaluating a high-fidelity model on a small number of samples properly selected from a large number of evaluations of a low-fidelity model. In particular, we will consider the class of multiscale transport models recently introduced in Bertaglia, Boscheri, Dimarco Pareschi, Math. Biosci. Eng. (2021) and Boscheri, Dimarco Pareschi, Math. Mod. Meth. App. Scie. (2021) as the high-fidelity reference and use simple two-velocity discrete models for low-fidelity evaluations. Both models share the same diffusive behavior and are solved with ad-hoc asymptotic-preserving numerical discretizations. A series of numerical experiments confirm the validity of the approach.



page 20

page 23


A bi-fidelity stochastic collocation method for transport equations with diffusive scaling and multi-dimensional random inputs

In this paper, we consider the development of efficient numerical method...

Multi-fidelity methods for uncertainty propagation in kinetic equations

The construction of efficient methods for uncertainty quantification in ...

GAN for time series prediction, data assimilation and uncertainty quantification

We propose a new method in which a generative adversarial network (GAN) ...

Neural Network Training Using ℓ_1-Regularization and Bi-fidelity Data

With the capability of accurately representing a functional relationship...

Transport away your problems: Calibrating stochastic simulations with optimal transport

Stochastic simulators are an indispensable tool in many branches of scie...

Multifidelity Orbit Uncertainty Propagation using Taylor Polynomials

A new multifidelity method is developed for nonlinear orbit uncertainty ...

A Stochastic Operator Approach to Model Inadequacy with Applications to Contaminant Transport

The mathematical models used to represent physical phenomena are general...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Mathematical modeling in epidemiology has certainly experienced an impressive increase in recent times driven by the overwhelming effects of the COVID-19 pandemic. A large part of the research has been directed towards the construction of models capable of describing specific characteristics associated with the pandemic, in particular the presence of asymptomatic individuals, which were largely underestimated, especially in the early stages of the spread of the disease [13, 22, 37, 35, 41]. Another line of research was aimed at the construction of models capable of describing the spatial characteristics of the epidemic, in order to be able to properly assess the impact of containment measures, in particular with regard to the mobility on the territory and restrictions in areas with higher risk of infection [10, 43, 42, 15, 39, 38]. Other approaches took into account the network structure of connections [8], social heterogeneity aspects [17], and the multiscale nature of the pandemic [5]. See also [36, 31] and the survey [1] for some recent developments in epidemiological modelling based on the use of kinetic equations.

Regardless of the characteristics of the model, a common aspect concerns the uncertainty of the data, which, especially in the first phase of the pandemic, largely underestimated the number of infected individuals. More precisely, the difficulty in correctly identifying all infected individuals in the early stages of the pandemic, due to structural limitations in performing large-scale screening and the inability to track the number of contacts, required the introduction of stochastic parameters into the models and the construction of related techniques for quantifying uncertainty [2, 3, 9, 6].

Among the various techniques of uncertainty quantification, the approaches based on stochastic strategies that do not necessarily require a-priori knowledge of the probability distribution of the uncertain parameters, as needed in the case of methods based on generalized polynomial chaos 

[46], are particularly interesting in view of a comparison with experimental data. On the other hand, the low convergence rate of Monte Carlo type sampling techniques poses serious limitations to their practical use.

In this context, multi-fidelity methods have shown to be able to efficiently alleviate such limitations through control variate techniques based on an appropriate use of low-fidelity surrogate models able to accelerate the convergence of stochastic sampling [48, 34, 29, 19, 20]. Specifically, general multi-fidelity approaches for kinetic equations have been developed in [19, 20], while bi-fidelity techniques with greedy sample selection in [29, 30]. We refer also to the recent survey in [16].

The present paper is devoted to extend the bi-fidelity method for transport equations in the diffusive limit developed in [29] to the case of compartmental systems of multiscale equations designed to model mobility dynamics in an epidemic setting with uncertainty [10, 6]. For this purpose, the corresponding hyperbolic system recently introduced in [8, 9] will be used as a reduced low-fidelity model. The two models allow to correctly describe the hyperbolic dynamics of the movement of individuals over long distances together with the small-scale, high-density, diffusive nature typical of urban areas [10]. In addition, both models share the same diffusive limit in which it is possible to recover classical models of diffusive type [25]. From a numerical viewpoint, a space-time asymptotic-preserving discretization, which work uniformly in all regimes, has been adopted in combination with the bi-fidelity approach. This permits to obtain efficient stochastic asymptotic-preserving methods.

The rest of the manuscript is organized as follows. In Section 2 we introduce the epidemic transport model with uncertainty in the simplified SIR compartmental setting, together with its diffusive limit. The corresponding reduced order two-velocity model used as low-fidelity surrogate is also discussed. Next, we extend the previous modelling to more realistic compartmental settings based on the introduction of the exposed and asymptomatic compartments in Section 3. The details of the bi-fidelity method and the asymptotic-preserving IMEX Finite Volume scheme are then given in Section 4. Section 5 contains several numerical experiments that illustrate the performance of the bi-fidelity approach. Some conclusions are reported at the end of the manuscript.

2 Epidemic transport model with random inputs

For simplicity, we first illustrate the modelling in the case of a classic SIR compartmental dynamic and subsequently we will extend our arguments to a more realistic SEIAR model, designed to take into account specific features of the COVID-19 pandemic, in Section 3.

2.1 The high-fidelity epidemic transport model

Let us consider a random vector

characterizing possible sources of uncertainty due to the independent stochastic parameters , which may affect variables of the mathematical model, as well as parameters or initial conditions. Individuals at position at time moving with velocity are denote by , and , which are the respective kinetic densities of susceptible (individuals who may be infected by the disease), infectious ( individuals who may transmit the disease) and removed (individuals healed or died due to the disease). We assume to have a population with subjects having no prior immunity and neglect the vital dynamics represented by births and deaths due to the time scale considered.

The kinetic distribution is then given by

and we recover their total density by integration over the velocity space

As a consequence,


with , denote the density fractions of the population at position and time that are susceptible, infected and recovered respectively. In this setting, the kinetic densities satisfy the epidemic transport equations [10]


where , , , take into account the heterogeneities of geographical areas, and are thus chosen dependent on the spatial location. Similarly, also the relaxation times , and . The quantity is the recovery rate of infected, while the transmission of the infection is governed by an incidence function modeling the transmission of the disease [24]


with the classic bi-linear case corresponding to , , even though it has been observed that an incidence rate that increases more than linearly with respect to the number of infected can occur under certain circumstances [14, 4, 28]. The parameter

characterizes the average number of contacts per person per time, multiplied by the probability of disease transmission in a contact between a susceptible and an infectious subject; whereas

acts as an incidence damping parameter based on the self-protective behavior of the individual that arises from awareness of the epidemic risk as the disease progresses [44, 21]. Notice that both and may vary in time as a consequence of governmental control actions, such as mandatory wearing of masks or full lockdowns, and the increasing awareness of the epidemic risks among the population.

The standard threshold of epidemic models is the well-known reproduction number , which defines the average number of secondary infections produced when one infected individual is introduced into a host population in which everyone is susceptible [24]. This number determines when an infection can invade and persist in a new host population. For many deterministic infectious disease models, an infection begins in a fully susceptible population if and only if . Assuming no inflow/outflow boundary conditions in , integrating over velocity/space and summing up the second equation in (2) we have



The above quantity, therefore, defines the stochastic reproduction number for system (2) describing the space averaged instantaneous variation of the number of infectious individuals at time . This definition naturally extends locally by integrating over any subset of the computational domain if one ignores the boundary flows.

2.1.1 Macroscopic formulation and diffusion limit

Let us introduce the flux functions


Then, integrating the system (2) against , it is straightforward to get the following set of equations for the macroscopic densities of commuters


whereas the flux functions satisfy


By introducing the space dependent diffusion coefficients


which characterize the diffusive transport mechanism of susceptible, infectious and removed, respectively, and keeping the above quantities fixed while letting the relaxation times to zero, we get from the r.h.s. in (2)

and consequently, from (7), we obtain a proportionality relation between the fluxes and the spatial derivatives (Fick’s law):


Thus, substituting (9) into (6) we get the diffusion system [33, 40, 45, 25]


We remark that the model’s capability to account for different regimes, hyperbolic or parabolic, accordingly to the space dependent values , , , makes it suitable for describing the dynamics of populations composed of human beings. Indeed, it is clear that the daily routine is a complex mixing of individuals moving at the scale of a city and individuals moving among different urban centers. In this situation, it seems reasonable to avoid, due to the lack of microscopic information and the high complexity, the description of the details of movements within an urban area and to describe this aspect through a diffusion operator. On the other hand, individuals when moving from one city to another follow well established connections for which a hyperbolic setting is certainly more appropriate.

2.2 The low-fidelity two-velocity epidemic model

The low-fidelity model, is based on considering individuals moving in two opposite directions (indicated by signs “+” and “-”), with velocities for susceptible, for infectious and for removed, we can describe the spatio-temporal dynamics of the population through the following two-velocity epidemic model [8]


In the above system, individuals , and are defined as

The transmission of the infection is governed by the same incidence function as in the high-fidelity model, defined in (3). Also the definition of the basic reproduction number results the same previously introduced in (4).

2.2.1 Macroscopic formulation and diffusion limit

If we now introduce the fluxes, defined by


we obtain a hyperbolic model equivalent to (11), but presenting a macroscopic description of the propagation of the epidemic at finite speeds, for which the densities follow system (6) and equations of fluxes read


Let us now consider the behavior of the low-fidelity model in diffusive regimes. To this aim, we introduce the diffusion coefficients


As for the previous model, the diffusion limit of the system is formally recovered letting the relaxation times , while keeping the diffusion coefficients (14) finite. Under this scaling, from (13) we get equations (9), which inserted into (6) lead again to the parabolic reaction-diffusion system (10).

Remark 1.

An important aspect in the bi-fidelity approach here proposed, is that the high-fidelity model and the low-fidelity one exactly coincide in the diffusive limit. The only difference lays in the definition of the diffusion coefficients (see eqs. (8) and eqs. (14)). As a consequence, in such a regime the bi-fidelity method achieves the maximum accuracy.

3 Extension to more epidemic compartments

To account for more realistic models to analyze the evolution of the ongoing COVID-19 pandemic, we consider extending the simple SIR compartmentalization by taking into account two additional population compartments, and , resulting in a SEIAR model [9, 6]. Subjects in the compartment are the exposed, hence infected but not yet infectious, being in the latent period. Moreover, among the infectious subjects, we distinguish the population between a group of individuals who will develop severe symptoms and a group of individuals who will never develop symptoms or, if they do, these will be very mild. This feature turns out to be essential to correctly analyze the evolution of COVID-19. In fact, it has been shown that individuals belonging to the group are very difficult to detect and isolate, contributing more strongly to the spread of the virus than the more easily detectable individuals [22, 35, 41].

3.1 The high-fidelity SEIAR transport model

Let us now consider a population at position moving with velocity directions , still taking into account possible uncertainties related to the random vector . Defining the kinetic densities of susceptible , exposed , severe symptomatic infected , mildly symptomatic or asymptomatic infected and removed (healed or deceased) , the kinetic distribution of the population results

In addition to (1), we have that

In this setting, the kinetic densities satisfy the following transport equations [6]


The quantities and are the recovery rates of symptomatic and asymptomatic infected (inverse of the infectious periods), respectively, while represents the inverse of the latency period and is the probability rate of developing severe symptoms [41, 22, 13]. In this model, the transmission of the infection is governed by two different incidence functions, and , simply to distinguish between the behavior of and individuals. Analogously to (3),


where a different contact rate, , and coefficient are taken into account for mildly/no symptomatic people. For the derivation of the reproduction number of this SEIAR kinetic model, which results


the reader is invited to refer to [9, 6].

When introducing the same definition (5) of flux for the additional compartments, and , integrating system (2) in , we get the following set of equations for the macroscopic densities


and for the fluxes


Moreover, defining also and and following the same procedure discussed in Section 2.1.1, we recover this SEIAR system in the diffusive regime:


3.2 The low-fidelity SEIAR transport model

The low-fidelity model is obtained, also with the more complex SEIAR compartmentalization, considering the discrete-velocity case with only two opposite velocities  [6]:


When defining fluxes and as in (12), the equivalent hyperbolic model underlying the macroscopic formulation of the spatial propagation of an epidemic is obtained. The evolution of the densities follows (18), while for the fluxes we have


As for the previous cases, the diffusion limit of the system is formally recovered letting the relaxation times , while keeping the diffusion coefficients, , finite. Thus, the same procedure presented for the SIR compartmentalization in Section 2.2.1 lead to the SEIAR parabolic system (20).

4 An asymptotic-preserving bi-fidelity numerical method

In this Section, we present the details of the numerical method used to solve the stochastic problem following the bi-fidelity asymptotic-preserving scheme proposed in [34, 48]

. For the high-fidelity model, the numerical scheme is structured in agreement with a discrete ordinate method in velocity with the even and odd parity formulation 

[18, 27]. Both the high-fidelity and low-fidelity solvers use Finite Volume Method (FVM) in space and achieve asymptotic preservation in time using suitable IMEX Runge-Kutta schemes [11, 12]. This permits to obtain a numerical scheme able to deal with the diffusion limit of the mathematical models without loosing consistency, and for which the time step size of the temporal discretization is not subject to excessive restrictions related to the smallness of the scaling parameters . For simplicity, we illustrate the numerical method in the case of the simpler SIR model. The application of the same numerical scheme results straightforward for the case of the SEIAR compartmentalization.

4.1 Asymptotic-preserving IMEX Finite Volume scheme

In this Section, we present the AP-IMEX Finite Volume scheme adopted to solve the SIR model at each stochastic collocation point selected for the bi-fidelity approximation.

The asymptotic-preserving IMEX method and the corresponding even and odd parities formulation was introduced in [10] for an SIR kinetic transport model in 2D domains. According to [27, 23], for , we can define the even and odd parities for the high-fidelity SIR kinetic transport model (2.1) as follows:

An equivalent formulation of (2.1) can be written as




The above densities can be approximated by a Gauss-Legendre quadrature rule. This leads to a discrete velocity setting, usually referred to as the discrete ordinate method, where we approximate

where and are the standard Gauss-Legendre quadrature weights and points in , and , number of chosen discrete velocities.

Assuming for simplicity of notation that , we can write the above system in the following compact form



Following [12], the Implicit-Explicit (IMEX) Runge-Kutta discretization that we consider for system (25) consists in computing the internal stages