Spatial spread of COVID-19 outbreak in Italy using multiscale kinetic transport equations with uncertainty

1.   Introduction

The advent of the COVID-19 pandemic has caused a strong commitment of many researchers acting in different fields with the scope of trying to understand and give explanations to the global crisis we are experiencing. From the mathematical modeling point of view, several progresses have been done concerning the development of epidemic models capable of taking into account the different facets of this terrible disease. In particular, many recent researches have been addressed to the search of control strategies ^[1,2,3,4] to limit the spread and consequently hospitalizations and deaths, possibly reducing, at the same time, as much as possible the negative impact on the economy of the restrictive measures ^[5,6].

Most of the recently proposed model represent improvements at various levels of the seminal works on compartmental epidemiological modeling proposed originally by Kermack and McKendrick ^[7]. These approaches ^{[1,3,8,9,10,11,12,13]} are typically focused on the epidemiological aspects of the virus spread under the hypothesis of global mixing of the population, hence, without taking into account the role of the spatial component in the evolution of an outbreak.

Despite in many situations the above description is sufficient to delineate the global trend of an epidemic, there are cases in which the spatial homogeneity assumption does not hold true. In these situations, one seeks for local interventions in order to reduce the spread without eventually affecting regions where the incidence of the infection does not require special care, as for instance in the recent case of the COVID-19 in Italy ^[14]. Consequently, from the modelling point of view, the inclusion of the spatial dependence represents a key challenge ^[15].

Indeed, with the increasing amount of information on population mobility and the computational resources available today, the design and simulation of epidemic models based on partial differential equations (PDEs) that include the details of spatial dynamics can be considered a realistic goal ^[16]. We recall that most of the existing epidemiological models taking into account spatial heterogeneities are based on reaction-diffusion equations ^{[17,18,19,20,21,22,23,24,25,26]}. Alternative modelling approaches are represented by the interaction of different homogeneous populations ^[4] or agent-based dynamics ^[27].

Most of these models do not allow for a clear distinction about the possible spatial behaviors of the population inside a given compartment. Consequently, although capable of originating realistic spatial patterns in situations where individuals move indiscriminately through the domain, such approaches are likely to be less effective in the case where one is interested in studying the spread of a virus in the human population. Under these circumstances, it is more realistic to consider only commuting individuals moving in major and preferred directions, not considering overall mass migration between distinct urban areas. Indeed, not all individuals move indiscriminately in the region of interest, as most of the population only interacts at the urban scale. Furthermore, in contrast to the use of diffusion models, the propagation rate of the infection along the spatial domain is obviously finite. Recently, trying to overcome some of the above criticisms, such as the paradox of the infinite propagation speed, alternative models based on hyperbolic PDEs have been proposed ^{[28,29,30,31,32]}.

In this work, following the approach introduced in ^[32], we consider a realistic compartmental structure for the description of the epidemic dynamics of commuters and non commuters individuals in presence of uncertain data. The model consists of a system of kinetic transport equations describing a large-scale (extra-urban) commuting population by a continuous density ^{[33,34,35,36]} in a two-dimensional environment. This density can be interpreted as the probability for an individual to be at a given location and move in a given direction at a given instant of time, in analogy to particle flows in rarefied gases ^{[37,38,39,40,41,42]}. The above system is coupled to a second system consisting of a set of diffusion equations that characterize the movement of the non-commuting population on a small (urban) scale. The epidemic spread is ruled by a Susceptible, Exposed, Infected, Asymptomatic and Removed (SEIAR) compartmental structure, in which both commuters and non commuters can interact. Using a suitable scaling process ^[32,43], the model allows us to highlight, firstly, the relationship with well-known existing approaches based on reaction-diffusion equations and, secondly, to pass naturally from a hyperbolic description in peripheral areas to a diffusive regime when reaching an urban area, with regard to the commuting population.

An important aspect concerns overcoming the limitation caused by standard deterministic models that rely on the assumption that the initial conditions, boundary conditions, and all involved epidemic and mobility parameters are known. However, as observed in the case of the COVID-19 epidemic, this assumption, especially in the early stages of the epidemic, is not reliable. For example, the initial conditions in terms of the number of infected and asymptomatic persons are certainly affected by uncertainty because data are limited and population screening cannot be error-free. Epidemic parameters, although normally estimated or calibrated, are also often candidates for being random variables in a realistic approach. Therefore, in order to take into account these limitations underlying deterministic models, in this paper we resort to a stochastic approach based on the introduction of random terms into the initial modeling ^[1,2,30,44].

Once the model is set up, its numerical solution on a computational domain describing a realistic geographic scenario poses several difficulties. In fact, the model consists of two coupled systems of PDEs, each characterized by five unknown functions living in a multidimensional space characterized by space, velocity, and stochastic variables. Additionally, we have to deal with the irregular shape of the spatial region and the multiscale nature of the dynamics (indeed, as previously stated, hyperbolic and parabolic behaviors coexist). Therefore, a particular care is devoted to the development of efficient and accurate numerical schemes. More precisely, a discretization of the system based on Gaussian quadrature points in velocity space ^[45,46] and a finite volume approach on unstructured grids ^[47,48] is considered. The adoption of asymptotic preserving time discretization techniques permits to avoid time step limitations introduced by the parabolic scaling without degradation of accuracy ^[49,50,51]. Finally, a non-intrusive stochastic collocation method, which guarantees spectral accuracy in the space of the uncertain parameters, is considered to deal with the uncertainty quantification process ^[30,52].

The rest of the paper is organized as follows. In Section 2 the mathematical model is introduced. We first introduce the kinetic transport formulation for the epidemic compartments of commuters together with the corresponding diffusive dynamics of the non commuters in a deterministic setting. Subsequently, we link the two hyperbolic/parabolic dynamics through a formal passage to the limit for the system of commuters. A definition of the basic reproduction number of the epidemic for the resulting model is also reported. Next, we illustrate how to generalize the model in presence of uncertainty. The details of the numerical scheme used to approximate the resulting stochastic system are summarized in Appendix A.1. Section 3 is devoted to present an application of the current modelling to the first outbreak of COVID-19 in Italy and its spread in the Lombardy Region. The capability of the model to accurately represent the first wave of the COVID-19 epidemic in Italy is discussed in detail through comparisons with recorded data reported by official sources ^[53]. Conclusions and future perspectives are finally given in Section 4.

2.   A compartmental kinetic transport model

Let $\Omega\subset {\mathbb{R}}^2$ characterizes a two-dimensional geographical area of interest and assume that individuals have been separated into two different groups: a commuting population, typically moving over long distances (extra-urban), and a non-commuting population, moving only in small-scale urban areas. In the first part of this Section, for ease of presentation, the multiscale kinetic transport model is introduced in a deterministic setting. The relation between the current hyperbolic transport model and classical diffusion models is discussed in the second part. In the third part, details regarding the basic reproduction number associated with the model are given. Finally, in the last part, we discuss the generalization of the deterministic model to the case where uncertainty is taken into account.

2.1.   Characterizing commuter and non commuter dynamics

We consider a population of commuters at position $x\in\Omega$ moving with velocity directions $v \in \mathbb{S}^1$ and denote the respective kinetic densities of susceptible (individuals who may be infected by the disease) by $f_S = f_S(x, v, t)$, exposed (individuals in the latent period, which are not yet infectious) by $f_E = f_E(x, v, t)$, severe symptomatic infected by $f_I = f_I(x, v, t)$, mildly symptomatic or asymptomatic infected by $f_A = f_A(x, v, t)$ and removed (individuals healed or died due to the disease) by $f_R = f_R(x, v, t)$. The kinetic distribution of commuters is then given by

and their total density is obtained by integration over the velocity space

As a consequence, one can recover the number of susceptible, exposed and recovered irrespective of their direction of displacement by integration over the velocity space. This gives

which we refer to as the density fractions of non-infectious individuals, whereas

are the density fractions of infectious individuals.

In this setting, the kinetic densities of the commuters satisfy the transport equations

where the total densities of non infected individuals are defined by

and similarly the total densities of infected by

In the above equations, $S_u(x, t)$, $E_u(x, t)$, $I_u(x, t)$, $A_u(x, t)$, $R_u(x, t)$ are the density fractions of non-commuters who, by assumption, move only on an urban scale. These densities satisfy a diffusion dynamics acting only at the same local scale

In the resulting Eqs (2.1) and (2.2) that couples the commuting and non-commuting dynamics, the velocities $v_i = \lambda_i(x) v$ in Eq (2.1), as well as the diffusion coefficients $D^u_i = D^u_i(x)$ in Eq (2.2), with $i\in\{S, E, I, A, R\}$, are designed to take into account the heterogeneity of geographical areas, and are thus chosen dependent on the spatial location. Similarly, also the relaxation times $\tau_i = \tau_i(x)$, $i\in\{S, E, I, A, R\}$ are space dependent. The quantities $\gamma_I = \gamma_I(x)$ and $\gamma_A = \gamma_A(x)$ are the recovery rates of symptomatic and asymptomatic infected (inverse of the infectious periods), respectively, while $a(x)$ represents the inverse of the latency period and $\sigma(x)$ is the probability rate of developing severe symptoms ^[3,4,9].

The transmission of the infection is governed by the incidence functions $F_I(\cdot, I_T)$ and $F_A(\cdot, A_T)$. We assume local interactions to characterize the nonlinear incidence functions ^[54,55]

where the classic bi-linear case corresponds to $p = 1$, $k_I = k_A = 0$. Parameters $\beta_I = \beta_I(x, t)$ and $\beta_A = \beta_A(x, t)$ characterize the contact rates of highly symptomatic and mildly symptomatic/asymptomatic infectious individuals, accounting for both the number of contacts and the probability of transmission. Hence, they may vary based on the effects of government control actions, such as wearing of masks, shutdown of specific activities or lockdowns ^[1,7,12]. On the other hand, parameters $\kappa_I = \kappa_I(x, t)$ and $\kappa_A = \kappa_A(x, t)$ are the incidence damping coefficients based on the self-protective behavior assumed by the individuals due to the awareness of the epidemic risk ^[11,26,30].

Alternative incidence functions are given by

where $\bar\Omega$ is a chosen portion of the domain which permits to take into account the fact that social distancing may depend on the average level of infection of the region rather than only on the local situation. Since the proposed model separates the contribution of the spread of the disease between the two infectious compartments, through the definition of two different incidence functions, $F_I$ and $F_A$, we chose to keep the same distinction also for the factors multiplying the damping coefficients $\kappa_I$ and $\kappa_A$ in each function respectively. Clearly, the choice could have been different, accounting for the sum $I_T + A_T$ in both the incidence functions, or even considering only the contribution of $I_T$ for both. The resulting Eqs (2.1) and (2.2) will be referred to as multiscale kinetic SEIAR (MK-SEIAR) model. Note that, because of the presence of two populations acting at different scales, the model allows a more realistic description of the typical commuting dynamics involving only a fraction of the population and distinguishes it from the epidemic process affecting the entire population.

2.2.   Commuters behavior in urban areas

The hyperbolic transport model for the commuters deserves some remarks. In fact, while it is clear that a hyperbolic description permits to describe correctly the daily extra-urban commuting part, the same individuals when moving inside the urban area are better described by a traditional diffusion model. A remarkable feature of the transport Eq (2.1), is that it permits to recover a classical diffusion behavior under the hypothesis that the relaxation times $\tau_{S, I, R}$ tend to zero while keeping finite the diffusion coefficients

More precisely, let us introduce the flux functions

Then, integrating system Eq (2.1) in $v$, we get the following set of equations for the macroscopic densities of commuters

whereas the flux functions satisfy

Clearly, the above system is not closed because the evolution of the fluxes in Eq (2.7) involves higher order moments of the kinetic densities. The diffusion limit can be formally recovered, by introducing the space dependent diffusion coefficients Eq (2.5) and letting $\tau_{S, I, R} \to 0$. We get from the r.h.s. in Eq (2.1)

and, consequently, from Eq (2.7) we recover Fick's law

since

and similarly for the other densities. Thus, substituting Eqs (2.8) into (2.6) we get the diffusion system for the population of commuters ^[21,23,56]

which is coupled with Eq (2.2) for the non-commuting counterpart. The capability of the model to account for different regimes, hyperbolic or parabolic, accordingly to the space dependent relaxation times $\tau_i$, $i\in\{S, E, I, A, R\}$, makes it suitable for describing the dynamics 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, due to the lack of microscopic information and the high complexity of the dynamics, it is reasonable to avoid describing the details of movements within an urban area and model this through a diffusion operator. On the other hand, commuters when moving from one city to another follow well-established connections for which a description via transport operators is more appropriate.

2.3.   Basic reproduction number

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 ^[7] during its entire period of infectiousness. This number determines when an infection can invade and persist in a new host population. For many deterministic epidemic models, an infection begins in a fully susceptible population if and only if $R_0 > 1$. Its definition in the case of spatially dependent dynamics is not straightforward, particularly when considering its spatial dependence. In the following, assuming no inflow/outflow boundary conditions in $\Omega$, integrating over velocity and space, we derive the following definition for the average reproduction number value on the domain $\Omega$

The derivation of the above expression for $R_0(t)$, computed following the next-generation matrix approach ^[57], is presented in detail in ^[30] using a suitable linearization of the corresponding nonlinear process for the space averaged quantities.

It is worth to underline that from definition Eq (2.10) it can be deduced that it is a combination of the growth of $E_T, I_T$ and $A_T$ that determines the persistence of the epidemic, not solely the growth of $E_T$ in time, neither the growth of the simple sum $E_T+I_T+A_T$. If, additionally, compartments $I$ and $A$ are considered homogeneously mixed in a unique compartment, allowing $\beta_I = \beta_A = \beta$, $\kappa_A = \kappa_I = \kappa$ and $\gamma_I = \gamma_A = \gamma$, we recover a SEIR-type compartmental model and the reproduction number results as in ^[32]:

Let us finally observe that, under the same no inflow/outflow boundary conditions, integrating in $\Omega$ Eqs (2.1) and (2.2) yields respectively the conservation of the total populations of commuters and non-commuters

2.4.   Including data uncertainty

To extend the Eqs (2.1) and (2.2) to the case in which uncertainties are taken into account, let us suppose that the population of commuters depend on an additional random vector $\boldsymbol{z} = (z_1, \ldots, z_{d})^T \in \Omega_{\boldsymbol{z}} \subseteq \mathbb{R}^{d}$, where $z_1, \ldots, z_{d}$ are independent random variables. This vector is used to characterize possible sources of uncertainty in the physical system due to lack of information on the actual number of infected or specific epidemic characteristics of the infectious disease.

Thus, in the system we have the following high-dimensional unknowns

The same considerations apply to the non-commuter population, yielding

Notice that besides the introduction of a new vector of variables, the structure of the Eqs (2.1) and (2.2) does not change, i.e., there is no direct variation of the unknowns with respect to $\boldsymbol{z}$, which, instead, have to be intended as parameters into the equations. We will further assume that also some epidemic parameters are affected by uncertainty. Therefore, for instance, parameters acting inside the incidence function may have an additional dependence of the form

In the next section, we discuss several numerical examples based on the Eq (2.1) and (2.2) with uncertainty. We will consider that the initial number of detected infected $I$, derived from the data at disposal ^[53,58], represents only a lower bound while the true values are not known but affected by uncertainty. Consequently, also initial conditions of exposed $E$, asymptomatic $A$ and susceptible $S$ will contain a stochastic dependence.

3.   Application to COVID-19 spread in Italy

To validate the proposed methodology in a realistic geographical and epidemiological scenario, a numerical test reproducing the epidemic outbreak of COVID-19 in the Lombardy Region of Italy, from February 27, 2020 to March 22, 2020, is designed, taking into account the uncertainty underlying initial conditions of infected individuals. We underline that the discretization of the multiscale system of PDEs Eqs (2.1) and (2.2) presented in this work is not trivial and requires the construction of a specific numerical method able to correctly describe the transition from a convective to a diffusive regime in realistic geometries. Additionally the method should be capable to deal efficiently with the uncertainty characterized by the stochastic nature of the model. Details on these numerical aspects are given in Appendix A.1 together with the tables of the data used for population mobility in Appendix A.2.

3.1.   Computational setup

The computational domain is defined in terms of the boundary that circumscribes the Lombardy Region, which is available in ^[58] as a list of georeferenced points in the ED50/UTM Zone 32N reference coordinate system. To avoid ill-conditioned reconstruction matrices and other related problems that arise when dealing with large numbers in finite arithmetic, all coordinates are re-scaled by a factor of $10^6$. The resulting computational grid is composed of $N_E = 10, 792$ triangular control volumes. The mesh grid is presented in Figure 1a. No-flux boundary conditions are imposed in the whole boundary of the domain, assuming that the population is not moving from/to the adjacent Regions. The domain is then subdivided in the $\mathcal{N}_c = 12$ provinces of Lombardy: Pavia (PV), Lodi (LO), Cremona (CR), Mantua (MN), Milan (MI), Bergamo (BG), Brescia (BS), Varese (VA), Monza-Brianza (MB), Como (CO), Lecco (LC) and Sondrio (SO). The identification of these cities is shown in Figure 1b.

The units of measure chosen for this numerical test can be summarized as follows:

with [L], [P] and [T] being the length, person and time units used in the simulation, respectively. Notice that the normalization of the population is made with respect to the total number of individuals of the Region (taken from ^[59]), which is $M = 10.027.602$, to properly work in a context in which the total population is equal to the unit.

To avoid the mobility of the population in the entire territory and to simulate a more realistic geographical scenario in which individuals travel along the main traffic paths of the Region, different values of propagation speeds are assigned in the domain which reflect, as close as possible, the real characteristics of the territory. Along the main connections of the Region, a mean value of $\lambda = 0.04$ is prescribed for compartments $S, E, A$ and $R$, which ensures a maximum travel distance of 80 km within a day; while, for the same compartments, $\lambda = 0.02$ is ulterior fixed in the urbanized circles. A spatial width of $h = 0.5$ km is assigned to the traveling paths. On the other hand, assuming that highly infectious subjects are mostly detected in the most optimistic scenario, being subsequently quarantined or hospitalized, the speed assigned to compartment $I$ is set null. However, the infected people, even if limited by quarantines and social distancing, can still contribute to the spread the disease via the diffusion process at the urban scale (mimicking for instance the still possible infections happening at the family level). A null value $\lambda = 0$ is set in the rest of the computational domain for all the epidemiological compartments. The resulting distribution of the characteristic speeds is visible from Figure 1c.

The relaxation time is set $\tau_i = 10^4$ for $i\in\{S, E, I, A, R\}$, so that the model recovers a hyperbolic regime in the entire region, apart from the main cities, where a parabolic setting is prescribed to correctly capture the diffusive behavior of the disease spreading which typically occurs in highly urbanized zones. Hence, to obtain a smooth change of the relaxation time between extra-urban and urban scale, in the urban area of radius $r_c$ of each city $c = 1, \ldots, \mathcal{N}_c$ the following relaxation time $\tau_{c, i}$, with Gaussian-like shape centered in the city center $(x_c, y_c)$, is prescribed:

with a diffusive relaxation time chosen to be $\tau_0 = 10^{-4}$. The resulting distribution of the relaxation times is presented in Figure 1d.

3.2.   Uncertain initial data and epidemic parameters

Considering $p_c$ the number of citizens of a generic city (province) denoted with subscript $c$, the initial spatial distribution of the generic population $f(x, y)$ is assigned, for each province and each epidemiological compartment, as a multivariate Gaussian function with the variance being the radius of the urban area $r_c$:

with $(x_c, y_c)$ representing the coordinates of a generic city center. The initial population setting, for each province of the Lombardy Region, is taken from ^[59] and reported in Table A.1, Appendix A.2. Note that, the radius $r_c$ associated to each city, defined in Table A.1 (first column), permits to exactly re-obtain the population $p_c$ when integrating over the computational domain the initial spatially distributed population.

Since at the beginning of the pandemic, tracking of positive individuals in Italy was very scarce, in this numerical test we consider that the initial amount of infected people is the leading quantity affected by uncertainty. To this aim, we introduce a single source of uncertainty $z$ having uniform distribution, $z \sim \mathcal{U}(0, 1)$ so that the initial conditions for compartment $I$, in each control volume, are prescribed as

with $I_{T, 0}$ initial amount of highly infectious corresponding to the values reported by February 27, 2020 in the GitHub repository ^[53] daily updated by the Civil Protection Department of Italy for each city, listed in the last column of Table A.1. This choice effectively associates all infected individuals detected with the $I$ compartment, as a result of the screening policy adopted during February–March 2020 in Italy. In fact, tests to assess the presence of SARS-CoV-2 were performed almost exclusively on patients with consistent symptoms and fever at the beginning of this pandemic. Regarding the uncertainty of these data, we impose $\mu = 1$, assuming that at least half of the actual highly symptomatic infected were detected at the beginning of the pandemic outbreak. For all the cities with zero infected detected by February 27, 2020 (e.g., Mantua, Varese, Como and Lecco), we choose to fix $I_{T, 0} = 1$ in order to assign an effective uncertainty.

Based on the estimations reported in ^[30], the expected initial amount of exposed $E_{T, 0}$ and asymptomatic/mildly symptomatic individuals $A_{T, 0}$ is imposed so that $E_{T, 0} = 10\, I_{T, 0}$ and $A_{T, 0} = 9\, I_{T, 0}$ in each location. Therefore, also initial conditions for compartments $E$, $A$ and $S$ become stochastic, depending on the initial amount of severe infectious at each location:

Finally, removed individuals are initially set null everywhere in the network, $R_T(0, z) = 0.0$.

To properly subdivide the population in commuters and non-commuters, regional mobility data are considered. In particular, the matrix of commuters presented in Table A.2, Appendix A.2, reflects mobility data provided by the Lombardy Region for the regional fluxes of year 2020 ^[60], which is in agreement with the one derived from ISTAT data released in October, 2011, as also confirmed in ^[4]. Therefore, to each control volume, the total percentage of commuters referred to the province where it is located is assigned, and non-commuting individuals are computed as a result of conservation principles. From Table A.2 it can be noticed that some connections are not taken into account simply because the amount of commuters along these routes is negligible if compared to the amount of individuals traveling in the other paths and with respect to the dimension of the populations.

Concerning epidemiological parameters, accordingly with values reported in ^[4,9,30], we set $\gamma_I\; = \; 1/14$, $\gamma_A = 2\gamma_I = 1/7$, $a = 1/3$, considering these parameters as clinical ones and therefore deterministic. We also assume the probability rate of developing severe symptoms $\sigma = 1/12.5$, as in ^[9,30].

From the first day simulated in this test, the population was aware of the risk associated with COVID-19 and recommendations such as washing hands often, not touching their faces, avoiding handshakes and keeping distances had already been disseminated by the government, hence, we initially fix coefficients $k_I = k_A = 50$.

The initial value of the transmission rate of asymptomatic/mildly infectious people is calibrated as the result of a least square problem with respect to the observed cumulative number of infected $I_T(t)$, through a deterministic SEIAR ODE model set up for the whole Lombardy Region, which provide the estimate $\beta_A = 0.58\times 10^{-3}$. As previously mentioned, since we are assuming that highly infectious subjects are mostly detected in the most optimistic scenario, being subsequently quarantined or hospitalized, the transmission rate of $I$ is set $\beta_I = 0.03\, \beta_A$, as in ^[4,9,30]. Finally, in the incidence function we fix $p = 1$.

With this parametric setup, we obtain an initial expected value of the basic reproduction number for the Region, evaluated as from definition Eq (2.10), $R_0(0)\; = \; 3.2$, which is in agreement with estimations reported in ^[4,9,30,61]. Nevertheless, with the proposed methodology it is possible to present the heterogeneity underlying the basic reproduction number at the local scale, as shown in Figure 2 for $\mu = 0$, together with the initial global amount of infected people $E_{T, 0}(x, y)+I_{T, 0}(x, y)+A_{T, 0}(x, y)$ present in the domain.

To model the effects of the lockdown imposed by the government from March 9, 2020 in north of Italy ^[14], from that day the transmission rate $\beta_A$ is reduced by 50%, also increasing the coefficients $k_I = k_A = 80$ as a result of the population becoming increasingly aware of the epidemic risks. In addition, the percentage of commuting individuals is reduced by 60% for each compartment according to mobility data tracked through GPS systems of mobile phones and made temporarily available by Google ^[61,62].

3.3.   Results and discussion

Numerical results of the test are reported in Figures 3–7. In Figure 3, the expected evolution in time of the infected individuals, together with 95% confidence intervals, is shown for exposed $E$, highly symptomatic subjects $I$ and asymptomatic or weakly symptomatic people $A$, for each city of the Lombardy Region. Here it is already appreciable the heterogeneity of the diffusion of the virus. Indeed, from the different y-axis scales adopted for the plot of the provinces, it can be noticed that Milan, Bergamo and Brescia present a consistently higher contagion.

From Figure 4 it can be observed that the lower bound of the confidence interval of the cumulative amount in time of highly symptomatic individuals is in line with data reported by the Civil Protection Department of Italy ^[53]. As expected, due to the uncertainty taken into account, the mean value of the numerical result in each city is higher than the registered one.

The comparison between the expected evolution in time of the cumulative amount of severe infectious with respect to the effective cumulative amount of total infectious people, including asymptomatic and mildly symptomatic individuals, is shown in Figure 5. From this figure it is clear that the number of infections recorded during the first outbreak of COVID-19 in Lombardy represents a clear underestimation of the actual trend of infection suffered by the Region and by Italy as a whole, and how the presence of asymptomatic subjects, not detected, has affected the pandemic evolution.

Numerical results presented in terms of integrated variables for the whole Lombardy are shown in Figure 6, together with the expected temporal evolution of the reproduction number $R_0(t)$, again with 95% confidence bands. It is here highlighted that the drop of $R_0(t)$ on March 9 reflects the imposition of lockdown restrictions, as presented in the previous Section. Moreover, in this plot, results are reported starting from February 28, instead of February 27, to permit to the system to achieve a correct initialization of the commuters (who are totally placed in the origin location at the beginning of the simulation) during the first day simulated.

In Figure 7, final expectation and variance of the cumulative amount of infected people, namely $E_T+A_T+I_T$, are reported in the 2D framework (first row). If comparing Figure 7a with 2a, it can be noticed that, at the end of March, the virus is no longer majorly affecting the province of Lodi and Cremona, but has been spread arriving to hit most of all Brescia, Milan and Bergamo. Finally, in the second row of Figure 7, the expectation of the susceptible population $S_T$ on the initial day simulated (February 27, 2020) is compared with the one obtained at the end of the simulation (March 22, 2020). Here it can be verified that the majority of the population does not leave their home city, as per real behavior, but there is only a small percentage of commuters who move within the domain, along the prescribed routes. Similar results are obtained for compartments $E$ and $A$, whose commuting part, even though small, strongly contributes to the spatial spread of the epidemic.

4.   Conclusions

In this paper we introduced a realistic model for the spatial spread of a virus with a focus on the case of COVID-19. Unlike models currently in the literature, which typically ignore spatial details or alternatively introduce them as simple diffusive dynamics, in our approach we have tried to capture the essential characteristics of the movements of individuals, which are very different if we consider commuting individuals who for work reasons move over long distances, from one city to another, to individuals who instead carry out their activities on an urban scale. The separation of individuals into these two classes, and the use of different spatial dynamics characterized by appropriate systems of transport and diffusive equations, allows in particular to avoid mass migration phenomena typical of models based on a single population and the instantaneous propagation of infectious disease over long distances. From the epidemiological point of view, modeling is developed in a compartmental context described by a SEIAR-type framework capable of describing the effect of exposed and asymptomatic individuals within the spatial spread of the disease. In addition, given the high uncertainty on the actual amount of individuals in the various compartments able to propagate the infection, the model was developed taking into account the presence of stochastic variables that therefore require an appropriate process of quantification. The resulting multiscale system of partial differential equations was then solved on realistic spatial geometries by a numerical method combining finite volume IMEX techniques for the deterministic part, with a non-intrusive collocation approach for the stochastic component. After a careful calibration of the model parameters based on the available data, an in-depth analysis of the results is reported in the case of the initial phase of the spread of COVID-19 in Italy occurred in the Lombardy region. The results show the ability of the model to correctly describe the epidemic dynamics and the importance of a heterogeneous spatial description and of the inclusion of stochastic parameters.

Acknowledgments

This work has been written within the activities of the GNCS group of INdAM (National Institute of High Mathematics). The support of MIUR-PRIN Project 2017, No. 2017KKJP4X "Innovative numerical methods for evolutionary partial differential equations and applications" is acknowledged.

Conflict of interest

The Authors declare that they have no conflict of interest.

A.   Supplementary material

In this appendix, we report the details of the numerical scheme adopted to approximate the MK-SEIAR Eqs (2.1) and (2.2) together with the data tables concerning the details of the population and of the commuting flows.

A.1.   Numerical method

We first give the details of the method in the case in which uncertainty is not present and successively we explain how the Eqs (2.1) and (2.2) is solved in the case of stochasticity. For the commuters, the numerical scheme for the deterministic case is based on a discrete ordinate method in velocity in which the even and odd parity formulation is employed ^[46,51,63]. The details of such approach are given in A.1.1. Then a finite volume method working on two-dimensional unstructured meshes ^[48,64] for the discrete ordinate approximation of the commuters is introduced in A.1.2. The full discretization of the Eq (2.1) is obtained through the use of suitable IMEX Runge-Kutta schemes ^[49,50]. In particular, the above choices permit to obtain a scheme which consistently captures the diffusion limit from the kinetic system when the scaling parameters $\tau_{S, I, R}$ tends toward zero. This part of the method is discussed in A.1.3. Finally, the discretization of the stochastic part for the Eq (2.1) and (2.2) when uncertainty is present is explained in A.1.4.

A.1.1.   Even and odd parities formulation

We rewrite Eq (2.1) by using the so-called even and odd parities formulation. In other words, we denote $v = (\eta, \xi)\in\mathbb{S}^1$ and then we obtain four equations with non-negative $\xi, \eta \geq 0$ for each compartment of the commuters. The change of variables reads, omitting the time and space dependence for simplicity, as ^[46]

while for the scalar fluxes one has

with $i = S, E, I, A, R$. An equivalent formulation with respect to Eq (2.1) can then be obtained thanks to this change of variables and reads as

and

One can observe that due to symmetry, we need to solve these equations for $\xi$, $\eta$ in the positive quadrant only. Thus the number of unknowns in Eq (2.1) and Eqs (A.1) and (A.2) is effectively the same.

A.1.2.   Space discretization on unstructured grids

We consider now a spatial two-dimensional computational domain $\Omega$ which is discretized by a set of non overlapping polygons $P_i, i = 1, \dots N_p$. Each element $P_i$ exhibits an arbitrary number $N_{S_i}$ of edges $e_{j, i}$ where the subscripts indicates that this is the edge shared by elements $P_i$ and $P_j$. The boundary of the cell is consequently given by

$\partial P_i = \bigcup \limits_{j = 1}^{N_{S_i}}{e_{ji}}$.

The governing equations for the commuters rewritten in the odd and even formulation are then discretized on the unstructured mesh by means of a finite volume scheme which is conveniently rewritten in condensed form as

where $\mathbf{Q}$ is the vector of conserved variables

while $\mathbf{F}(\mathbf{Q})$ is the linear flux tensor and $\mathbf{S}(\mathbf{Q})$ represents the stiff source term defined in Eqs (A.1) and (A.2). As usual for finite volume schemes, data are represented by spatial cell averages, which are defined at time $t^n$ as

where $|P_i|$ denotes the surface of element $P_i$.

A first order in time finite volume method is then obtained by integration of the governing Eq (A.3) over the space control volume $|P_i|$

Higher order in space is then achieved by substituting the cell averages by piecewise high order polynomials. We refer to these polynomial reconstructions to as $\mathbf{w}_i(\mathbf{x})$ which are obtained from the given cell averages Eq (A.4). In particular, our choice is to rely on a second order Central WENO (CWENO) reconstruction procedure along the lines of ^[64,65]. We omit the details for brevity. The numerical flux function $\mathbf{F}_{ij} \cdot \mathbf{n}_{ij}$ is given by a simple and robust local Lax-Friedrichs flux yielding

where $\mathbf{w}_{i, j}^+, \mathbf{w}_{i, j}^-$ are the high order boundary extrapolated data evaluated through the CWENO reconstruction procedure. The numerical dissipation is given by $s_{\max}$ which is the maximum eigenvalue of the Jacobian matrix in spatial normal direction,

Let us notice now that in the diffusion limit, i.e., as $(\tau_S, \tau_I, \tau_R)\to 0$, the source term $\mathbf{S}(\mathbf{Q})$ becomes stiff, therefore in order to avoid prohibitive time steps we need to discretize the commuters system implicitly. To this aim, a second order IMEX method which preserves the asymptotic limit given by the diffusion Eq (2.9) is proposed and briefly described hereafter.

A.1.3.   Time integration and numerical diffusion limit

We consider again Eq (2.1) formulated using the parities Eqs (A.1) and (A.2). We also assume $\tau_{S, I, R} = \tau$ and rewrite Esq (A.1) and (A.2) in partitioned form as

in which

and $\mathbf{f}({\mathbf{u}})$, $\mathbf{g}({\mathbf{u}})$, ${\mathbf{E}}({\mathbf{v}})$ are defined similarly. Now, following ^[50], the Implicit-Explict Runge-Kutta (IMEX-RK) approach appid to system (A.8) reads as

where ${\mathbf{u}}^{(k)}, {\mathbf{v}}^{(k)}$ are the so-called internal stages. The numerical solution reads

In the above equations, the 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, we choose the Runge-Kutta scheme in such a way that the following relations hold true

The scheme Eqs (A.10) and (A.11) treats implicitly the stiff terms and explicitly all the rest. Moreover, one can prove that the above scheme is a consistent discretization of the limit system in the diffusive regime. In fact, assuming for simplicity $D_{S, I, R}$ independent from space, the second equation in Eq (A.10) can be rewritten as

therefore, assuming Eq (2.5), the limit $\tau\to 0$ yields

where ${\mathbf{D}} = {\rm diag}\left\{D_S, D_S, D_E, D_E, D_I, D_I, D_A, D_A, D_R, D_R\right\}$ and where we used the fact that from the first equation in Eq (A.10) as $\tau\to 0$ we have ${\mathbf{u}}^{(j)} = {\mathbf{U}}^{(j)}$. Note also that (A.13) implies that $j_{S, E, I, A, R}^{(1)} = j_{S, E, I, A, R}^{(2)}$ in ${\mathbf{v}}^{(j)}$, i.e., we restore perfect symmetry in direction of propagation of the information. Using now the identity ${\mathbf{u}}^{(j)} = {\mathbf{U}}^{(j)}$ into the first equation in Eq (A.10) we get

and using Eqs (A.13) into (A.14) thanks to the definitions of $\mathbf{f}$ and $\mathbf{g}$ gives

Finally, integrating over the velocity field one has

and thus, the internal stages correspond to the stages of the explicit scheme applied to the reaction-diffusion Eq (2.9). To conclude the proof one has to notice that thanks to the choice Eq (A.12), the last stage is equivalent to the numerical solution. Thus, this is enough to guarantee that the scheme is a consistent discretization of the limit equation.

Note that the limit system is consistent with the discretization of the non commuters diffusive Eq (2.2). In this latter case, we adopt the same finite volume setting for the unknowns

This simply reads

with

Then, the same CWENO reconstruction and the same numerical local Lax-Friedrichs flux is employed where however, we account for a dissipation proportional to the diffusive terms. In other words, the numerical viscosity is given by the maximum eigenvalue of the viscous operator $s_{max}^V = \max \left(D^u_S, D^u_E, D^u_I, D^u_A, D^u_R \right)$. Concerning the time discretization, the explicit part of the Runge-Kutta scheme introduced in the previous paragraph is employed.

A.1.4.   Stochastic collocation method

In the case in which we deal with the stochastic Eqs (2.1) and (2.2), we employ a generalized Polynomial Chaos (gPC) expansion technique ^[52,66]. We restrict to the case in which there is only one stochastic variable $z$ in the system. The extension to the case of a vector of random variables is straightforward. The probability density function of the single random input is supposed known: $\rho_z: \Gamma \rightarrow \mathbb{R}^+$. In this case, the approximate solution for the commuters $\mathbf{Q}_M(x, v, t, z)$ and the non commuters $\mathbf{Q}_M^{u}(x, t, z)$ are represented as truncations of the series of the orthonormal polynomials describing the random space, i.e..

where $M$ is the number of terms of the truncated series and $\phi_j(z)$ are orthonormal polynomials, with respect to the measure $\rho_z(z)\, d z$. The expansion coefficients are obtained by

Then, the exact integrals for the expansion coefficients in Eq (A.19) are replaced by a suitable quadrature formula characterized by the set $\{z_n, w_n \}_{n = 1}^{N_p}$, where $z_n$ is the $n$-th collocation point, $w_n$ is the corresponding weight and $N_p$ represents the number of quadrature points. We then have

where ${\mathbf{Q}}(x, v, t, z_n)$ and ${\mathbf{Q}}^{u}(x, t, z_n)$ with $n = 1, \ldots, N_p$ are the solutions of the problem evaluated at the $n$-th collocation point for the commuters and non commuters. Thanks to the computation of the above coefficients than

it is possible to compute all quantities of interest concerning the random variable. For example, the expectations

are approximated as

and

In the same way, all other quantities of interest such as the variance of $\mathbf{Q}$ and $\mathbf{Q}^{u}$ can be computed.

A.2.   Population and mobility data

In this appendix we report tables containing data on the initial distribution of populations in the various urban areas considered in the Lombardy region, see Table A.1, and data on the corresponding flows of commuters between cities, see Table A.2.