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 , social heterogeneity aspects , and the multiscale nature of the pandemic . See also [36, 31] and the survey  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, 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 .
The present paper is devoted to extend the bi-fidelity method for transport equations in the diffusive limit developed in  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 . In addition, both models share the same diffusive limit in which it is possible to recover classical models of diffusive type . 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 vectorcharacterizing 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 
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 
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; whereasacts 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 . 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):
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 
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).
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 
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
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 :
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
4 An asymptotic-preserving bi-fidelity numerical method
. 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  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.