A local sensitivity and regularity analysis for the Vlasov-Poisson-Fokker-Planck system with multi-dimensional uncertainty and the spectral convergence of the stochastic Galerkin method

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 [14,23,38,36,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 [13,43,42,16,1,39]. Other approaches took into account the network structure of connections [10], social heterogeneity aspects [21], and the multiscale nature of the pandemic [6]. See also [37,32] and the survey [2] 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 [3,4,9,7].

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,35,30,18,19]. Specifically, general multi-fidelity approaches for kinetic equations have been developed in [18,19], while bi-fidelity techniques with greedy sample selection in [30,31]. We refer also to the recent survey in [17].

The present paper is devoted to extend the bi-fidelity method for transport equations in the diffusive limit developed in [30] to the case of compartmental systems of multiscale equations designed to model mobility dynamics in an epidemic setting with uncertainty [13,7]. For this purpose, the corresponding hyperbolic system recently introduced in [10,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 [13]. In addition, both models share the same diffusive limit in which it is possible to recover classical models of diffusive type [26]. 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 modeling 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 models with random inputs

For simplicity, we first illustrate the modeling 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 $ {\bf z} = (z_1,\ldots,z_{d})^T \in {\mathbb{R}}^d $ characterizing possible sources of uncertainty due to the independent stochastic parameters $ z_1,\ldots,z_{d} $, which may affect variables of the mathematical model, as well as parameters or initial conditions. Individuals at position $ x\in\Omega\subset {\mathbb{R}} $ at time $ t $ moving with velocity $ v \in [-1,1] $ are denote by $ f_S = f_S(x,v,t,{\bf z}) $, $ f_I = f_I(x,v,t,{\bf z}) $ and $ f_R = f_R(x,v,t,{\bf z}) $, which are the respective kinetic densities of susceptible $ S $ (individuals who may be infected by the disease), infectious $ I $ (individuals who may transmit the disease) and removed $ R $ (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 $ \rho(x,t,{\bf z}) = S(x,t,{\bf z})+I(x,t,{\bf z})+R(x,t,{\bf z}) $, denote the density fractions of the population at position $ x $ and time $ t>0 $ that are susceptible, infected and recovered respectively. In this setting, the kinetic densities satisfy the epidemic transport equations [13]

where $ v_S = {\lambda_S}(x) v $, $ v_I = {\lambda_I}(x) v $, $ v_R = {\lambda_R}(x) v $, $ {\lambda_S},{\lambda_I},{\lambda_R} \geq 0 $ take into account the heterogeneities of geographical areas, and are thus chosen dependent on the spatial location. Similarly, also the relaxation times $ \tau_S = \tau_S(x) $, $ \tau_I = \tau_I(x) $ and $ \tau_R = \tau_R(x) $. The quantity $ \gamma = \gamma(x,{\bf z}) $ is the recovery rate of infected, while the transmission of the infection is governed by an incidence function $ F(\cdot,I) $ modeling the transmission of the disease [25]

with the classic bi-linear case corresponding to $ p = 1 $, $ \kappa = 0 $, even though it has been observed that an incidence rate that increases more than linearly with respect to the number of infected $ I $ can occur under certain circumstances [15,5,29]. The parameter $ \beta = \beta(x,t,{\bf z}) $ 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 $ \kappa = \kappa(x,t,{\bf z}) > 0 $ 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,22]. Notice that both $ \beta $ and $ \kappa $ 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 $ R_0 $, which defines the average number of secondary infections produced when one infected individual is introduced into a host population in which everyone is susceptible [25]. 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 $ R_0 > 1 $. Assuming no inflow/outflow boundary conditions in $ \Omega $, integrating over velocity/space and summing up the second equation in (2) we have

when

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 $ t>0 $. This definition naturally extends locally by integrating over any subset of the computational domain $ \Omega $ 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 $ v $, 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 $ \tau_{S,I,R} $ 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 [34,40,45,26]

We remark that the model's capability to account for different regimes, hyperbolic or parabolic, accordingly to the space dependent values $ \tau_S $, $ \tau_I $, $ \tau_R $, 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 epidemic transport model

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

In the above system, individuals $ S(x,t,{\bf z}) $, $ I(x,t,{\bf z}) $ and $ R(x,t,{\bf z}) $ 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 $ R_0 $ 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 $ \tau_{S,I,R}\to 0 $, 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, $ E $ and $ A $, resulting in a SEIAR model [9,7]. Subjects in the $ E $ 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 $ I $ who will develop severe symptoms and a group of individuals $ A $ 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 $ A $ group are very difficult to detect and isolate, contributing more strongly to the spread of the virus than the more easily detectable $ I $ individuals [23,36,41].

3.1.   The high-fidelity SEIAR transport model

Let us now consider a population at position $ x\in\Omega $ moving with velocity directions $ v \in [-1,1] $, still taking into account possible uncertainties related to the random vector $ {\bf z} $. Defining the kinetic densities of susceptible $ f_S = f_S(x,v,t,{\bf z}) $, exposed $ f_E = f_E(x,v,t,{\bf z}) $, severe symptomatic infected $ f_I = f_I(x,v,t,{\bf z}) $, mildly symptomatic or asymptomatic infected $ f_A = f_A(x,v,t,{\bf z}) $ and removed (healed or deceased) $ f_R = f_R(x,v,t,{\bf z}) $, 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 [7]

The quantities $ \gamma_I = \gamma_I(x,{\bf z}) $ and $ \gamma_A = \gamma_A(x,{\bf z}) $ are the recovery rates of symptomatic and asymptomatic infected (inverse of the infectious periods), respectively, while $ a = a(x,{\bf z}) $ represents the inverse of the latency period and $ \sigma = \sigma(x,{\bf z}) $ is the probability rate of developing severe symptoms [41,23,14]. In this model, the transmission of the infection is governed by two different incidence functions, $ F(\cdot,I) $ and $ F_A(\cdot,A) $, simply to distinguish between the behavior of $ I $ and $ A $ individuals. Analogously to (3),

where a different contact rate, $ \beta_A = \beta_A(x,t,{\bf z}) $, and coefficient $ \kappa_A = \kappa_A(x,t,{\bf z}) $ are taken into account for mildly/no symptomatic people. For the derivation of the reproduction number $ R_0 $ of this SEIAR kinetic model, which results

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

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

and for the fluxes

Moreover, defining also $ D_E = \frac13\lambda_E^2\tau_E $ and $ D_A = \frac13\lambda_A^2\tau_A $ 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 $ \pm \lambda $ [7]:

When defining fluxes $ J_E $ and $ J_A $ 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 $ \tau_{S,E,I,A,R}\to 0 $, while keeping the diffusion coefficients, $ D_i = \lambda_i^2 \tau_i, i\in\{S,E,I,A,R\} $, 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 [35,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 [20,28]. 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 $ \tau_i $. 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 [13] for an SIR kinetic transport model in 2D domains. According to [28,24], for $ v>0 $, we can define the even and odd parities for the high-fidelity SIR kinetic transport model (2) as follows:

An equivalent formulation of (2) can be written as

where

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 $ w_i $ and $ \zeta_i $ are the $ N_G $ standard Gauss-Legendre quadrature weights and points in $ [-1,1] $, and $ N_G = N_v $, number of chosen discrete velocities.

Assuming for simplicity of notation that $ \tau_{S,I,R} = \tau $, we can write the above system in the following compact form

with

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

and evaluating the final numerical solution

To properly compute the implicit terms $ \frac{\boldsymbol{R}^{(j)}}{2} $ at each Runge-Kutta internal step explicitly, we refer to [27]. Matrices $ \tilde A = (\tilde a_{kj}) $, with $ \tilde a_{kj} = 0 $ for $ j\geq k $, and $ A = (a_{kj}) $, with $ a_{kj} = 0 $ for $ j > k $ are $ s \times s $ matrices, with $ s $ number of Runge-Kutta stages, defining respectively the explicit and the implicit part of the scheme, and vectors $ \tilde b = (\tilde b_1, ...,\tilde b_s)^T $ and $ b = (b_1, ...,b_s)^T $ are the quadrature weights. Furthermore, referring to [11,12], if the following relations hold,

the method is said to be globally stiffly accurate (GSA). It is worth to notice that this definition states also that the numerical solution of a GSA IMEX Runge-Kutta scheme coincides exactly with the last internal stage of the scheme. Since the GSA property is fundamental to preserve the correct diffusion limit and to achieve asymptotic-preservation stability in stiff regimes [13,10], in the sequel, the GSA BPR(4, 4, 2) scheme presented in [12] is chosen, characterized by $ s = 4 $ stages for the implicit part, 4 stages for the explicit part and 2nd order of accuracy.

At each internal stage of the IMEX scheme (26), we apply a Total-Variation-Diminishing (TVD) Finite Volume discretization to evaluate the numerical fluxes [10,8]. To achieve second order accuracy also in space, while avoiding the occurrence of spurious oscillations, a classical minmod slope limiter has been adopted.

The same AP-IMEX Finite Volume scheme is adopted also to solve the low-fidelity SIR model (6)-(13), as fully presented in [10]. The reader is invited to refer to [10,9,7,13] for further details on the properties of the chosen numerical scheme applied to epidemic models.

4.2.   Bi-fidelity stochastic collocation

In this Section, we briefly review the bi-fidelity method developed in [35,48]. Bi-fidelity methods make use of low-fidelity models to effectively inform the selection of representative points in the parameter space and then employ this information to construct accurate approximations to high-fidelity solutions. To facilitate future discussion, we denote the expensive high-fidelity solution $ u^H({\bf{z}}) $ and the cheap low-fidelity solution $ u^L({\bf{z}}) $ for any given random parameter $ {\bf{z}}\in I_z \subset {\mathbb{R}}^d $, where $ I_z $ is the domain of the random parameter $ {\bf{z}} $.

The basic idea of the bi-fidelity approximation is to construct an inexpensive surrogate $ u^B({\bf z}) $ of the high-fidelity solution in the following non-intrusive manner:

where $ n $ is the number of the selected parameter points in the parameter space. If $ n $ is small and the coefficient $ c_k({\bf{z}}) $ can be efficiently and accurately approximated, an efficient bi-fidelity approximation can be constructed. Once the surrogate is constructed, the statistics of high-fidelity solutions can be quickly approximated by evaluating the bi-fidelity surrogate via Monte Carlo methods or quadrature rules.

There are two major key questions for the performance of the above bi-fidelity algorithms: (a) how to select the representative points $ {{\bf{z}}}_k $ effectively? (b) how to efficiently construct the bi-fidelity approximation for any new given $ {\bf{z}} $ but avoiding requiring a high-fidelity simulation?

Subset Selection. Existing predefined or structured nodes (e.g., sparse grids, cubature rules, etc.) often grows fast in high dimensions. Therefore, these options cannot easily accommodate the current situation where we would like the size $ n $ of $ \gamma_n $ to be small and also arbitrary. Alternatively, adaptive approaches based on high-fidelity models to explore the parameter space suit our needs. However, this typically requires a large number of high-fidelity samples, which might not be computationally affordable. In contrast, the low-fidelity model is inexpensive to evaluate and it mimics the variations of high-fidelity solutions in the parameter space. This motivates us to employ the inexpensive low-fidelity model to learn and explore the behaviors of the high-fidelity model in the parameter space.

We shall identify important points iteratively by a greedy approach [35,48]. Specifically, we denote a candidate sample set $ \Gamma_N = \{{\bf{z}}_1, {\bf{z}}_2, \cdots , {\bf{z}}_N\} $, which is assumed to be large enough to cover the parameter space $ I_z $. Initially, denote $ \gamma_0 = \{\} $ and assume we have the first $ k $ important points $ \gamma_k = \{{\bf{z}}_{i_1}, {\bf{z}}_{i_2}, \cdots , {\bf{z}}_{i_k}\} $ available at the $ k $-th iteration. Denote the snapshot matrix $ u^L(\gamma_k) = \{u({\bf{z}})| {\bf{z}} \in \gamma_k\} $ and the corresponding spanned approximation space $ U^L(\gamma_k) = {\rm{span}}\{u^L({\bf{z}})|{\bf{z}}\in\gamma_k\} $. The corresponding high-fidelity approximation space can be defined similarly, $ U^H(\gamma_n) = {\rm{span}}\{u^H({\bf{z}})|{\bf{z}}\in\gamma_n\} $. Then we pick the point (from the candidate set $ \Gamma_N $) so that the corresponding low-fidelity solution is farthest away from the existing spanned low-fidelity approximation space $ U^L(\gamma_k) $, to be the next sampling point:

where $ d^L(v,W) $ is the distance between a function $ v\in u^L(\Gamma_N) $ and the space $ W\in u^L(\gamma_k) $. We then repeat this step to select all $ n $ important points. The whole procedure can be efficiently implemented by apply the pivoted Cholesky decomposition on $ u^L(\Gamma_N) $ [35,48].

Bi-fidelity approximation. For any given $ z $, to efficiently compute bi-fidelity approximation in (28), it is desirable to find a cheap yet reasonably accurate approximation of $ c_k({\bf z}) $. The bi-fidelity approach developed in [35,48] learns these coefficients from the inexpensive low-fidelity model $ u^{L}({\bf z}) $. Specifically, for any given $ {\bf z} $, we shall compute the low-fidelity solution $ u^L({\bf{z}}) $ and then construct its best approximation in the spanned low-fidelity approximation space $ U^L(\gamma_n) $ by orthogonal projection:

where the low-fidelity projection coefficients can be computed as follows:

and $ {\bf G}^L $ is the Gramian matrix of $ u^L(\gamma_n) $,

where $ \langle\cdot,\cdot\rangle^L $ is the standard inner product associated with $ U^L(\gamma_n) $.

Under certain conditions [35,48], $ c^L_k({\bf{z}}) $ can be reasonably good approximation for the high-fidelity coefficients in (28). Consequently, we can construct the bi-fidelity approximation of the high-fidelity approximation solution $ u^H({\bf{z}}) $ as follows:

To put things together, we outline the major steps for the bi-fidelity approximation of the high-fidelity sample for a given $ {\bf{z}} $ in Algorithm 1, reported in the following.

Remark 2. As we mentioned above, the bi-fidelity method relies on the assumption that the low-fidelity coefficients are similar to the high-fidelity coefficients in the parameter space under certain conditions, as stated in [35,48]. This may not hold for some problems. In this case, a correction mapping between low-fidelity and high-fidelity coefficient can be constructed by leveraging approximation power of neural network to further improve the accuracy, if the additional high-fidelity data are available. We refer readers to [33] for details on this approach.

To construct the bi-fidelity approximation of the high-fidelity mean, with the bi-fidelity surrogate $ u^B({\bf z}) $, we can employ the Monte Carlo or other quadrature-based methods to compute the statistical moments quickly by evaluating bi-fidelity surrogates. To further reduce the number of bi-fidelity surrogate evaluations, a more efficient extension to this approach is developed in [47]. The general procedure is as follows:

● Compute the low-fidelity sample mean (via the Monte Carlo or quadrature rules):

where the high-fidelity mean $ \mu^H $ is defined similarly.

● Compute its best approximation on the low-fidelity approximation space $ U^L(\gamma_n) $ by orthogonal projection,

where the expansion coefficients $ c^L $ are computed by solving the following linear system:

● Using this coefficient $ c^L_{k} $ as the surrogate of the high-fidelity coefficient of the high-fidelity mean $ \mu^H $, the bi-fidelity approximation of the high-fidelity mean can be constructed as follows:

Note that in this way, only one bi-fidelity surrogate evaluation is required. The standard deviation can be computed similarly. We refer readers to [47] for additional details.

5.   Numerical examples

To examine the performance of the proposed methodology, two benchmark tests are considered: the first concerning the kinetic transport model with the SIR compartmentalization discussed in Section 2, and the second one regarding the extension to the SEIAR modeling presented in Section 3. It is worth to highlight that, in these tests, we collect all the quantities of interests together in a single vector to choose the $ {\bf{z}}_{k} $ points of the bi-fidelity algorithm, for either the 3 compartments $ S,I,R $ in Test 1 or the five compartments $ S $, $ E $, $ I $, $ A $, $ R $ in Test 2. From the experiments, we also find that the results are similar if we choose the $ {\bf{z}}_{k} $ points separately for the different epidemic compartments.

5.1.   Test 1: SIR model with heterogeneous environment

We first consider the initial distributions of the high-fidelity kinetic SIR model (2) as follows:

where $ c = \frac{1}{2} \sum_{i = 1}^{N_G} w_i \, e^{-\frac{\zeta_i^2}{2}} $ is a re-normalization constant, with $ N_G $ number of Gauss-Legendre quadrature points as defined in Section 4.1, and

in the physical domain $ L = [0,20] $. The initial fluxes $ J_S(x,0), J_I(x,0) $ and $ J_R(x,0) $ are null and we consider periodic boundary conditions. The same initial conditions for $ S, I, R $ and $ J_S, J_I, J_R $ are imposed in the low-fidelity SIR model (6)-(13). In the following tests, we consider a 2-dimensional random vector $ {\bf{z}} = (z_1, z_2)^T $, with independent random parameters $ z_1 $ and $ z_2 $ that follow a uniform distribution, $ z_j \sim \mathcal{U}(-1,1) $, $ j = 1,2 $.

We analyze the behavior of the proposed methodology concerning spatially heterogeneous environments, considering a spatially variable contact rate [10,44]

where

and a recovery rate

In the incidence function, we set $ \kappa = 0 $ and $ p = 1 $. Therefore, we simulate an infectious disease characterized by an $ R_0 $ perturbed around the value 1 depending on the random fluctuations. We discretize the spatial domain, both in the high-fidelity and low-fidelity models, with $ N_x = 150 $ cells and consider $ N_v = 8 $ velocities for the high-fidelity model. To compute the reference solutions for the mean and standard deviation, we use a $ 3 $-rd level sparse grid quadrature based on Clenshaw-Curtis rules for the choice of the stochastic collocation nodes, with a total of $ 29 $ points, for both the high-fidelity and low-fidelity models.

Test 1 (a): In this case, a parabolic configuration of speeds and relaxation parameters is considered, setting $ \lambda_i^2 = 10^5 $, $ i \in \{S,I,R\} $, and $ \tau_i = 10^{-5} $ in the low-fidelity model and $ \tau_i = 3\times 10^{-5} $ in the high-fidelity model, to maintain consistency of the two simulations (see Remark 1). The time step size results $ \Delta t = 0.89 \times 10^{-2} $ in both models, nevertheless the low-fidelity code is almost 5 times faster than the high-fidelity one ($ t_{CPU}^{HF} \approx 70.0 \,s $, $ t_{CPU}^{LF} \approx 14.3 \,s $ to simulate until $ t = 5 $ in a machine with processor Intel(R) Core(TM) i7-9750H CPU @2.60 GHz – 6 cores and 12 logical processors).

In the first row of Figure 1, the expectation and standard deviation of the solution of compartment $ I $ for the high-fidelity model, the low-fidelity model and the bi-fidelity approximation at time $ t = 5 $ are shown. As expected, since high-fidelity and low-fidelity models share the same diffusive limit, a perfect agreement of the solutions can be observed. To confirm this, in the second row of the same figure, we also plot the $ L^2 $ errors of the mean and standard deviation between the bi-fidelity and high-fidelity solutions at $ T = 5 $ with respect to the number of selected "important" points $ n $ of the bi-fidelity algorithm. A fast error decay is clearly observed. With only $ n = 8 $ hi-fidelity sample points, bi-fidelity approximation can achieved a relative error of $ \mathcal{O}(10^{-6}) $ for both the mean and standard deviation.

Test 1 (b): In this second case, we consider the hyperbolic regime by letting $ \lambda_i = 1 $, $ i \in \{S,I,R\} $, while $ \tau_i = 1 $ in the low-fidelity model and $ \tau_i = 3 $ in the high-fidelity model, for consistency. The time step size results $ \Delta t = 0.12 $ in both models, but the low-fidelity simulation is again almost 5 times faster than the high-fidelity one ($ t_{CPU}^{HF} \approx 5.6\, s $, $ t_{CPU}^{LF} \approx 1.3\, s $ to simulate until $ t = 5 $ in the same machine previously mentioned).

In Figure 2, the results of Test 1(b) are reported for the infectious compartment $ I $ at time $ t = 5 $. The first row shows expectation and standard deviation of the high-fidelity, the low-fidelity and the bi-fidelity solutions. While low-fidelity model fails to capture details of the hi-fidelity solutions around peaks due to the different velocity setting considered in the two models, bi-fidelity approximation are in almost perfect agreement with the high-fidelity ones, confirming the validity of the proposed methodology even with hyperbolic configurations of the scaling parameters. In addition, we plot the relative $ L^2 $ error convergence of the mean and standard deviation of the bi-fidelity and high-fidelity solutions at $ t = 5 $ with the number of high-fidelity samples, in the second row of Figure 2. Again, we observe a fast error decay for mean and standard deviation with $ \mathcal{O}(10) $ high-fidelity samples.

5.2.   Test 2: SEIAR model with distinguished epidemic hotspots

Next, we analyze the effectiveness of the proposed methodology also with the extended SEIAR compartmentalization examining a more realistic epidemic scenario. Let us now consider a 2-dimensional random vector $ {\bf{z}} $, this time with $ z_j \sim \mathcal{U}(0,1) $, $ j = 1,2 $. We design an initial condition for the low-fidelity SEIAR model (18)-(22) that simulates the presence of 3 cities aligned in the spatial domain $ L = [0,20] $ with a different number of exposed (infected but not yet infectious) individuals, subjected to randomness, distributed following a Gaussian function,

where $ x_1 = 10/3 $, $ x_2 = 10 $, $ x_3 = 50/3 $ are the coordinates of the city centers and $ \alpha_1 = 0.01(1+z_1) $, $ \alpha_2 = 0.001(1+z_1) $, $ \alpha_3 = 0.004(1+z_1) $ are the stochastic amplitudes. In this case, we set the large stochastic amplitude to model large uncertainty present on the initial amount of exposed people when an epidemic starts spreading. Indeed, it is worth to underline that the initial number of exposed individuals is certainly one of the variables most affected by uncertainty during epidemic outbreaks [23]. We consider that there are no infectious people at $ t = 0 $, being

and $ J_S(x,0) = J_E(x,0) = J_I(x,0) = J_A(x,0) = J_R(x,0) = 0 $, with periodic boundary conditions, to allow a connection also between cities 1 and 3. The initial distributions of the high-fidelity kinetic SEIAR model (15) then read as defined in (38), for $ i \in \{S,E,I,A,R\} $.

Concerning epidemic parameters, to simulate the more challenging scenario in which the incidence function presents sinusoidal oscillations in space as well as being greater in the most populated areas, we consider the following distribution of the contact rate related to asymptomatic individuals:

where $ \beta_A^0 $ is affected by the following random fluctuations:

Assuming that highly infectious subjects are mostly detected in the most optimistic scenario, being subsequently quarantined or hospitalized, we set $ \beta_I(x,{\bf{z}}) = 0.03\, \beta_A(x,{\bf{z}}) $. Then, we fix $ \gamma_I = 1/14 $, $ \gamma_A = 1/7 $, $ a = 1/3 $, $ \sigma = 1/12.5 $, considering these clinical parameters deterministic according to values adopted in [9,7,14,23] to simulate the COVID-19 spread. Finally, in the incidence functions, we set $ \kappa_I = \kappa_A = 0 $, hence assuming that initially individuals are not aware of the epidemic outbreak, and $ p = 1 $ to work with the standard bi-linear case. The resulting setting permits to simulate and epidemic characterized by a baseline reproduction number $ R_0 = 3 $. We consider $ N_x = 150 $ cells for both high-fidelity and low-fidelity simulations, and use $ N_v = 8 $ velocities for the high-fidelity solution.

Test 2 (a): In the first case, an intermediate regime between parabolic and hyperbolic is considered, setting $ \lambda_i^2 = 10 $, $ i\; \in\; \{S,E,A,R\} $, $ \lambda_I = 0 $ and $ \tau_i = 0.25 $, $ i\; \in\; \{S,E,I,A,R\} $, in the low-fidelity model and $ \tau_i = 0.75 $ in the high-fidelity model. The characteristic speed of compartment $ I $ is fixed to zero because we assumed that infectious people with severe symptoms are generally detected and isolated [9,7].

The baseline space-time evolution of the compartments is shown in Figure 3, where the fast propagation of the epidemic can be observed, especially starting from the first city on the left, as expected. In Figure 4, the expectation and standard deviation of infected individuals, divided by each compartment, are presented for each methodology adopted, highlighting the validity of the bi-fidelity approach. In Figure 5, $ L^2 $ error decay of the bi-fidelity approximations, with respect to the number of selected "important" points $ n $, is compared with the error that would be obtained simply using the low-fidelity model.

Test 2 (b): In the second case, we consider a fully hyperbolic regime choosing $ \lambda_i^2 = 1 $, $ i\; \in\; \{S,E,A,R\} $, $ \lambda_I = 0 $, and $ \tau_i = 10 $, $ i\; \in\; \{S,E,I,A,R\} $, and $ \tau_i = 30 $ respectively in the low-fidelity and high-fidelity model.

The baseline space-time evolution of the compartments of major interest is presented in Figure 6. In this figure, the spatial heterogeneity of the epidemic spread related to the prescribed contact rate can be clearly appreciated. Mean and standard deviation obtained with the 3 approaches are shown for compartments $ E $, $ I $ and $ A $ in Figure 7. Again, although the low-fidelity and high-fidelity solutions show clear discordant trends, the effectiveness of the bi-fidelity method is plainly confirmed. In Figure 8, the decay of the $ L^2 $ relative error norm of the bi-fidelity approximation with respect to the number of points $ n $ shows that a proper accuracy can be achieved for both mean and standard deviation even in this challenging test case.

6.   Conclusions

In this work we introduced a bi-fidelity method for the quantification of uncertainty in epidemiological transport models based on an asymptotic-preserving space-time discretization. In detail, after presenting the high-fidelity epidemiological model, we considered the corresponding bi-fidelity model. Both models share the same diffusive limit and permits to recover classical epidemic models based on diffusion equations. The numerical scheme used allows to have an efficient quantification of the uncertainty in different regimes, using few simulations of the high-fidelity model and several runs of the low-fidelity model for points selection in the random space. Results for one-dimensional transport problems based on realistic compartmentalization in relation to the recent COVID-19 pandemic, which also include asymptomatic individuals, show the validity of the presented approach. Further research will be directed toward extending the present approach to realistic contexts such as those studied in [7,9].

Acknowledgments

G.B. and L.P. acknowledge the support of MIUR PRIN 2017, Project No. 2017KKJP4X Innovative numerical methods for evolutionary partial differential equations and applications. G.B. holds a Research Fellowship from the Italian National Institute of High Mathematics, INdAM (GNCS). L.L. holds the Direct Grant for Research supported by Chinese University of Hong Kong and Early Career Scheme 2021/22, No. 24301021, from Research Grants Council of Hong Kong. X.Z. was supported by the Simons Foundation (504054).