key: cord-0773947-52wghbi4
authors: Bertaglia, Giulia; Pareschi, Lorenzo
title: Hyperbolic compartmental models for epidemic spread on networks with uncertain data: application to the emergence of Covid-19 in Italy
date: 2021-05-29
journal: Mathematical Models & Methods in Applied Sciences
DOI: 10.1142/s0218202521500548
sha: d177e096fbab96b62d71ccffd00954219cc04278
doc_id: 773947
cord_uid: 52wghbi4

The importance of spatial networks in the spread of an epidemic is an essential aspect in modeling the dynamics of an infectious disease. Additionally, any realistic data-driven model must take into account the large uncertainty in the values reported by official sources, such as the amount of infectious individuals. In this paper we address the above aspects through a hyperbolic compartmental model on networks, in which nodes identify locations of interest, such as cities or regions, and arcs represent the ensemble of main mobility paths. The model describes the spatial movement and interactions of a population partitioned, from an epidemiological point of view, on the basis of an extended compartmental structure and divided into commuters, moving on a suburban scale, and non-commuters, acting on an urban scale. Through a diffusive rescaling, the model allows us to recover classical diffusion equations related to commuting dynamics. The numerical solution of the resulting multiscale hyperbolic system with uncertainty is then tackled using a stochastic collocation approach in combination with a finite-volume IMEX method. The ability of the model to correctly describe the spatial heterogeneity underlying the spread of an epidemic in a realistic city network is confirmed with a study of the outbreak of COVID-19 in Italy and its spread in the Lombardy Region.

The advent of the COVID-19 pandemic, which still afflicts the entire world, has prompted many researchers in the mathematical field and beyond to propose epidemic models increasingly suitable for the study of the evolution of this specific coronavirus, to support the definition of the best strategies for the control of its spread. Most of the proposed models are rooted in the compartmental epidemiological modeling proposed by Kermack and McKendrick [21, 32] , based on systems of ordinary differential equations (ODE), which describe only the temporal evolution of the spread of the epidemic, neglecting the spatial component in favor of an assumption of homogeneity of population and territory [20, 27, 30, 31, 36, 40, 43, 44, 50, 52] . Recently, models have also been proposed in which we have moved away from the idea of a homogeneously mixed population to allow a better description in terms of social contact patterns, vulnerability and case fatality ratio, but still neglecting the spatial development of the epidemic spread [8, 11] .

Generally, the concept of the average behavior of a population is sufficient to have a first reliable description of the development of an epidemic. However, the inclusion of the spatial component in epidemiological systems is crucial especially when there is a need to consider spatially heterogeneous interventions, as was and still is the case for the control of the spread of COVID-19 [45, 49] . Some recent works have attempted to fill these gaps by proposing epidemiological models based on partial differential equations (PDE), in which the spatial structure is taken into account. In [17, 22, 53, 54] a two-dimensional (2D) space dependence is considered, obtaining models able to properly capture the complexity of the spatial transmission of the virus. On the other hand, meta-population network approaches have also been proposed [23, 29] , which result in systems of ODE that provide no information about the spatial characteristics of the epidemic spread along the pathways connecting the different populations [49] .

Another important aspect concerns the wide use of deterministic models, which, although more computationally efficient, are based on the assumption that initial conditions, boundary conditions and all the parameters involved are known. However, in practical applications, this assumption is rarely true. In the context of epidemic modeling, for example, initial conditions are certainly affected by uncertainty because data are limited and screening policy is always a matter of compromises. Also epidemic parameters, although normally estimated or calibrated, are candidates for being random variables. Therefore, when attempting to solve the problem of interest numerically, one must take into account these limitations, recurring to more consistent stochastic models [8, 23, 29, 41, 44] .

In this work, a hyperbolic transport model on networks with uncertain data capable to describe the spatial spread of an epidemic phenomenon defined by an extended compartmental dynamics is proposed. The model takes inspiration by the multiscale transport model recently introduced in [17] by extending it to include more realistic compartmentalizations. In addition, the spatial structure of the system lays on the network modeling originally adopted in [15] , in which nodes of the network identify locations of interests, depending on the spatial scale considered (villages, provinces, regions, nations), and arcs constitute the mobility paths connecting them (hence roads, railways, airline connections, etc.). Individuals moving on the network are subdivided into commuting and non-commuting, the latter acting only at the urban, nodal scale they belong to, not contributing in the transmission of the virus in the extra-urban domain. This particular feature prevents unrealistic mass migration effects where the entire population travels through the network [17] . Furthermore, thanks to its hyperbolic structure the model avoids the unphysical feature of infinite propagation speeds typical of reaction-diffusion systems, still recovering the parabolic behaviour in the diffusive limit [10, 15] .

The resulting model is solved numerically through an IMEX finite volume stochastic scheme that permits to preserve uniform accuracy with respect to the scaling parameters [15, 35] . More precisely, the numerical scheme combines a stochastic collocation method, a non-intrusive method which guarantees spectral convergence in the stochastic space and ease of implementation, avoiding the loss of important structural properties of the original system of governing equations [41, 48, 57] , with an Implicit-Explicit (IMEX) Runge-Kutta finite-volume scheme that maintains the consistency in the asymptotic limit (i.e. the asymptotic-preserving (AP) property) [16] .

The rest of the paper is organized as follows. In Section 2 the mathematical model is introduced, first in its kinetic transport formulation and then deriving its macroscopic representation. Subsequently, the diffusion limit is formally computed and the details of the model extension to networks are presented. Moreover, a definition of the basic reproduction number of the epidemic is proposed, with the details of its derivation discussed in Appendix A. In Section 3 we present an application to the first outbreak of COVID-19 in Italy and its spread in the Lombardy Region. First, to validate the proposed numerical approach a convergence study is performed. Next, the capacity of the proposed model to capture the spatial heterogeneous characteristics underlying COVID-19 epidemics is assessed through several simulations based on data reported by official sources. The details of the stochastic collocation approach are summarized in Appendix B, whereas the IMEX finite-volume scheme is described in Appendix C. Finally, conclusions and future perspectives are drawn in Section 4. The epidemiological starting point of our model is given by a compartmental structure with a SEIAR partitioning [44, 51, 52] , in which the population is assumed to be divided in susceptible individuals S, exposed individuals in the latent period, which are not yet infectious E, infectious individuals manifesting severe symptoms I, infectious individuals with no/mild symptoms A and removed (deceased or healed and immune) individuals R. We assume to have a population with subjects having no prior immunity. The vital dynamics represented by births and deaths is neglected because of the time scale considered. This model differs from the classic SEIR [32] for the presence of a subgroup of people who will never develop symptoms or very mild ones, which is an essential feature to take into account if one aims at analyzing the evolution of the COVID-19 pandemic. Indeed, individuals belonging to this group are very hard to be detected, deeply impacting the efficiency of the monitoring and making isolation, containment and tracing of individual cases very challenging [20, 29, 30, 44] . They tend to behave like simple susceptible, but in fact highly contributing in augmenting the spread of the virus.

To account for the spatial movement of the population, the model has its roots within discrete velocity models in kinetic theory [15, 39] , and follows the approach in [17] , where individuals of each compartment are subdivided in three classes, S ±,0 , E ±,0 , I ±,0 , A ±,0 , R ±,0 , traveling in a one-dimensional space domain Ω ⊆ R with characteristic speeds +λ i , −λ i and 0 respectively, with i ∈ {S, E, I, A, R}. The total compartmental densities are defined as the sum of all the components of the subgroups

The discrete-velocity system of the SEIAR epidemic transport model for commuters, associated to relaxation times τ i , i ∈ {S, E, I, A, R}, then reads

This system is coupled with the following SEIAR model describing the evolution of a stationary population of non-commuters

In the above system (2.2)-(2.3) all the epidemic densities S ±,0 , E ±,0 , I ±,0 , A ±,0 , R ±,0 depend on (x, t, z), where (x, t) are the physical variables of space x ∈ Ω ⊆ R and time t > 0, while z = (z 1 , . . . , z d ) T ∈ R d is a random vector characterizing the possible sources of uncertainty due to independent parameters z 1 , . . . , z d . The same applies for the contact rate function f , defined with respect to the infectious compartments, I and A respectively, as

where β I (x, t, z) and β A (x, t, z) are the transmission rates, accounting for both number of contacts and probability of transmission, hence may vary based on the effects of S ±,0 E ±,0 A ±,0 Figure 1 : Flow chart of the multi-population SEIAR dynamic based on five compartments: susceptible (S), exposed (E), severe symptomatic infectious (I), mildly symptomatic/asymptomatic infectious (A), and removed -healed or deceased-population (R), each one subdivided in three classes of individuals traveling in the domain with characteristic speeds +λ i , −λ i and 0, with i ∈ {S, E, I, A, R}. The transition rates between the compartments, β j , a, γ j are the inverse of the contact period 1/β j , the latent period 1/a and the infectious period 1/γ j , respectively, with distinct values for groups j ∈ {I, A}, except for the latent period. The rate of probability for the exposed to enter in the symptomatic and asymptomatic subgroups of the infectious population are σ and (1 − σ), respectively. Coefficients k j affect the contact rates in response to social distancing and other control actions.

government control actions, such as mandatory wearing of masks, shutdown of specific work/school activities, or full lockdowns [8, 30, 32, 54] ; k I (x, t, z) and k A (x, t, z) act as incidence damping coefficients based on the self-protective behavior of the individual that arises from awareness of the risk associated with the epidemic [15, 27] . Note that, the classic bilinear case corresponds to p = 1 and k I = k A = 0, even though it has been observed that an incidence rate that increases more than linearly with respect to the number of infectious can occur under certain circumstances [10, 21, 37] . Parameters γ I (x, t, z) and γ A (x, t, z) are the recovery rates of highly symptomatic infected and of asymptomatic or mildly infected (inverse of the infectious periods), respectively, while a(x, t, z) represents the inverse of the latency period and σ(x, t, z) is the probability rate of developing severe symptoms [20, 29, 52] . The flow chart of the SEIAR model considered in this work is also illustrated in Fig. 1 , where transition rates between compartments are clearly displayed.

Remark 2.1.

• The model allows to describe more realistically the typical dynamic of commuters, which regards only a small fraction of individuals, and to distinguish this from the epidemic process which, instead, involves the whole population, including noncommuters. The presence of a group of non-commuting population, indeed, permits to avoid that the whole population in a compartment moves indiscriminately in the full space originating an unrealistic mass migration effect.

• The presence of uncertainty in the data included from the beginning in the modeling process could allow the compartmentalization of asymptomatic individuals to be eliminated by implicitly including them in the uncertainty about the number of infected individuals. This approach was proposed in [8] . In our case, however, in order to highlight the link with similar models used in the literature, we chose to keep the asymptomatic compartment, which consequently will be the one affected by the highest level of uncertainty.

Introducing now the macroscopic variables S c , E c , I c , A c , R c for the commuters, with

and defining the fluxes

a hyperbolic model underlying the macroscopic formulation of the spatial propagation of an epidemic, equivalent to the mesoscopic one [9] , presented in system (2.2), is obtained:

Note that here the above system is coupled with the equations for the non commuter population (2.3) through identities (2.1).

Formally, system (2.4) is a 1D system of stochastic balance laws which can be rewritten in compact form as:

in which

It is easy to verify that system (2.5) is symmetric hyperbolic in the sense of Friedrichs-Lax [28] , with real finite characteristic velocities (eigenvalues) λ i , i ∈ {S, E, I, A, R}, and a complete set of linearly independent eigenvectors. All the eigenvectors are associated with genuinely non-linear fields, defining shocks and rarefactions, and the Riemann invariants of the system correspond to the kinetic transport variables

(2.6)

The standard threshold of epidemic models is the well-known basic reproduction number R 0 , also called the basic reproduction ratio or the basic reproductive rate, which defines the expected number of secondary cases produced, in a completely susceptible population, by a typical infected individual during its entire period of infectiousness [24, 32] . For many deterministic infectious disease models, this number determines whether an infection can invade and persist in a new host population (R 0 > 1) or cannot (R 0 < 1). Its definition in the case of spatially dependent dynamics, as already noted in [54] , is not straightforward particularly when considering its spatial dependence. In the following we will consider the following definition for the average value of the reproduction number on the domain Ω for t > 0:

The derivation of the above expression for R 0 (t), computed following the next-generation matrix approach [24] , is presented in detail in Appendix A.

In the two contributions of eq. (2.7) it is possible to identify:

• the transmission rates for compartment I and A: f I (S T , I T ) and f A (S T , A T ), respectively;

• the mean time of staying in compartment I and A: γ −1 I and γ −1 A , respectively; • the fraction of individuals passing from compartment E to I and A: aσ/a = σ and a(1 − σ)/a = 1 − σ, respectively.

It is worth to underline that from definition (2.7) 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 . Remark 2.2. If no spatial dependence is assigned to variables and parameters, hence the ODE version of system (2.4) is considered, and no social distancing effects are taken into account (i.e. k I = k A = 0) the reproduction number results in accordance with [29, 51, 52] :

From a formal viewpoint, it can be shown that the proposed model recovers the parabolic behavior expected from standard space-dependent epidemic models in the diffusion limit [10, 15] . Introducing the diffusion coefficients

that characterize the diffusive transport mechanism of S, E, I, A, R respectively, and letting τ i → 0, i ∈ {S, E, I, A, R}, while keeping the diffusion coefficients finite, from the last five equations of system (2.4) we recover Fick's law getting

Finally, inserting these results in the rest of the equations of system (2.4) yields the following parabolic reaction-diffusion system for the commuters

Therefore, the relaxation times can modify the nature of the behavior of the solution [10, 15, 16] , which can result either hyperbolic or parabolic (when considering small relaxation times and large speeds). This feature of the model makes it particularly suitable for the description of the dynamics of human populations, which are characterized by movement at different spatial scales [17] . It is therefore natural to assume τ i = τ i (x), i ∈ {S, E, I, A, R}, since in geographic areas densely populated we can assume a diffusive dynamics while in other areas or along the main arteries of communication a hyperbolic description will be more appropriate avoiding propagation of information at infinite speed.

In the model description adopted in this Section, the relaxation times are assumed to be space dependent but independent of the state variables. More generally, it is possible to consider the relaxation times

leading, asτ i → 0 and taking λ 2

A classical example is represented by the choice

which corresponds to a generalization of Carleman model. In the limitτ i → 0, for α = 0 we recover again the linear diffusion model whereas assuming α ∈ (−1, 0) we obtain a fast diffusion process considered by other authors in epidemiology [13] . From a mathematical viewpoint we refer to [39] for rigorous results concerning these kind of diffusion limits for generalized Carlemann models. Although interesting, here we do not explore further this direction.

The transport model here proposed can be extended to network approaches in the sense of those presented in [18, 19, 46, 47] . In the sequel we summarize the details of the network structure we adopted to characterize arcs and nodes with different sizes. Following [15] , it is possible to structure a 1D network considering that the nodes of the network identify locations of interest such as municipalities, provinces or, in a wider scale, regions or nations, while the arcs, enclosing the 1D spatial dynamics, represent the paths linking each location to the others. In this way, the epidemic state of each node evolves in time influenced by the mobility of the commuting individuals, moving from the other locations included in the network, always considering a part of the population composed by non-commuting individuals which remain at the origin node.

In order to prescribe the proper coupling between nodes and arcs, ensuring the conservation of total density (population) in the network and of fluxes at the interface, it is necessary to impose appropriate transmission conditions at each arc-node interface.

A network or a connected graph G = (N , A) is composed of a finite set of N nodes (or vertices) N and a finite set of A bidirectional arcs (or edges) A, such that an arc connects a pair of nodes [47] . Let us parametrize the A arcs of the network as intervals a i = [0, L i ], i = 1, . . . , A. Arcs are bidirectional, as the network is non-oriented, but an artificial orientation needs to be fixed in order to define a sign for the velocities and therefore the fluxes. For an incoming arc, L i is the abscissa of the node, whereas for an outgoing arc the same abscissa is 0.

To define transmission conditions at a generic node n ∈ N having a i ∈ A, i = 1, . . . , N a,L incoming arcs and a j ∈ A, j = 1, . . . , N a,0 outgoing arcs, we need to consider two kind of interfaces at the node: the interfaces neighboring incoming arcs (L, i) and the interface neighboring outgoing arcs (0, j). If variables are discontinuous across these interfaces, at time t + ∆t, for each of them (1 + N a,L ) new states originate at interfaces (L, i) and (1 + N a,0 ) new states originate at interfaces (0, j) [47] . To compute them, we need to solve (2 + N a,L + N a,0 ) Riemann problems, using the Riemann Invariants (kinetic variables) of the system, defined in Eqs. (2.6) , and the principle of conservation of fluxes at interfaces [18, 19] .

For each one of the compartments of commuting individuals of the SEIAR model here discussed, let us indicate, for ease of notation, with u the number of individuals of the compartment, with v the corresponding analytical flux, with λ its characteristic velocity, and with u ± the Riemann Invariants. Introducing transmission coefficients, α i,j ∈ [0, 1], which represent the probability that an individual at a generic arc-node interface decides to move across that interface, from the j-th location to the i-th location, transmission conditions at the interfaces with an incoming arc (L, i) results for the arcs side

and for the node side (with the subscript n indicating the variable -or the location, when concerning transmission coefficients-of the node)

(2.10)

On the other hand, for the transmission conditions at the interfaces with an outgoing arc (0, j), we have for the arcs side

(2.11) and for the node side

Notice that the condition differs when considering an incoming flow or an outgoing flow, due to the artificial orientation that has been set. Indeed, for each incoming arc, we need to use u + L,i from the arc and u − n from the node; while for each outgoing arc we consider u − 0,j from the arc and u + n from the node [19] . Furthermore, to guarantee the conservation of fluxes at the interface, ensuring that the global mass (population) of the system is conserved, the following must hold [18, 19] 

which is fulfilled imposing at interfaces (L, i)

and at interfaces (0, j)

Nodes located at the inlet (outlet) end of the domain are without any incoming (outgoing) arcs. At these nodes, in order to ensure that there are no individuals entering or leaving the network (thus ensuring the preservation of the total population), we simply enforce the standard no-flux boundary condition [47] , which consists of imposing at inlet nodes v * L,n = 0, u * L,n = u n − v n λ n , (2.16) and at outlet nodes v * 0,n = 0, u * 0,n = u n + v n λ n .

(2.17)

The multiscale transport SEIAR model (2.2)-(2.3) is solved using a second-order IMEX Runge-Kutta Finite Volume Collocation method (see Appendices B and C for details). In particular, we will show that the numerical scheme is capable to reach spectral accuracy in the stochastic space, if the solution is sufficiently smooth in that space, and to preserve this accuracy in the diffusive (stiff) limit (i.e. stochastic AP property) [14, 35] . An advantage related to the choice of the stochastic collocation method lies in its non-intrusive nature [57] . This feature ensures ease of implementation, since the method requires only the evaluation of the solution of the corresponding deterministic problem, followed by a postprocessing step. Thus, no major manipulation efforts of the deterministic computational code are required and the loss of important structural properties of the original problem is avoided [48] . . Each entry is given as number of people and percentage of individuals with respect to the total amount of commuters of the origin Province. The last column shows the amount of total commuters of the origin Province and corresponding percentage with respect to the total population of the origin Province. This matrix is extracted from the origin-destination matrix provided by the Lombardy Region for the regional fluxes of year 2020 [3] . 

To analyze the effectiveness of the proposed approach in a realistic geographical and epidemic scenario, we design a numerical setting reproducing the evolution of the first outbreak of COVID-19 in the Lombardy Region of Italy, from February 27, 2020 to March 27, 2020, with respect to uncertainties underlying the initial conditions and chosen epidemic parameters.

A five-node network is considered, whose nodes represent the 5 main provinces interested by the epidemic outbreak in the first months of 2020: Lodi (n 1 ), Milan (n 2 ), Bergamo (n 3 ), Brescia (n 4 ) and Cremona (n 5 ). The arcs a j connecting each node to the others identify the main set of routes and railways viable by commuters each day. A schematic representation of this network is shown in Fig. 2 . Routes that connect cities that are not direct neighbors and that require crossing other provinces are to be considered as the sum of the two route sections. Thus, for instance, it can be seen from Fig. 2 that the Milan-Brescia connection is actually identified by the sum of the Milan-Bergamo (a 2 ) and Bergamo-Brescia (a 3 ) arcs, indeed indicated as a 2+3 . Since the space unit of measurement adopted is 1 L = 10 2 km, the lengths of the single sections of arc result: L a 1 = 0.20, L a 2 = 0.30, L a 3 = 0.35, L a 4 = L a 5 = 0.40 . These arcs are discretized with a grid size ∆x a = 0.01. The length of the spatial cell associated to each node is proportional to the dimension of the urbanized area of the corresponding province, hence: ∆x n 1 = 0.025, ∆x n 2 = 0.420, ∆x n 3 = 0.100, ∆x n 4 = 0.085, ∆x n 5 = 0.035. Transmission coefficients α i,j , as well as percentage of commuters belonging to each province, are imposed recurring to official national assessment of mobility flow. In particular, the matrix of commuters presented in Table 1 reflects mobility data provided by the Lombardy Region for the regional fluxes of year 2020 [3] , which is in agreement with the one derived from ISTAT data released in October, 2011 [4] , as also confirmed in [29] . Notice that connections Lodi-Bergamo and Lodi-Brescia are not taken into account (as observable also from Fig. 2 ) because the amount of commuters along these routes is very low compared to the amount of individuals traveling the rest of the routes.

The characteristic speed associated to each arc is fixed to permit a full round trip in each origin-destination section within a day. Since the time unit of measurement adopted is 1 T = 1 day, for all the compartments we impose: λ a 1 = 0.4, λ a 2 = 0.6, λ a 3 = 0.7, λ a 4 = λ a 5 = 0.8, respectively for each arc. The characteristic speed of compartment I is fixed to zero in all the nodes of the network, λ n = 0, while for the rest of the compartments, namely S, E, A, R, λ n = 1 . In the arcs, the relaxation time is assigned so that the model recovers a hyperbolic regime, hence τ r,a = 10 3 . On the other hand, a parabolic setting is prescribed in the cities for commuters in order to correctly capture the diffusive behavior of the disease spread which typically occurs in highly urbanized zones, with τ r,n = 10 −3 . The reader is invited to refer to [15] for further details on the sensitivity of the model to relaxation and speed parameters.

We simulate the spread of COVID-19 starting from February 27, 2020 until March 27, 2020, included. At the beginning of the pandemic, tracking of positive individuals cannot be considered reliable. Thus, the initial amount of infected people is the first information affected by uncertainty. To this aim, we introduce a single source of uncertainty z having uniform distribution, z ∼ U(0, 1) and the initial conditions for compartment I, at each node, are prescribed as

with i 0 density of infectious people on February 27, 2020, as given by data recorded by the Civil Protection Department of Italy [2] , reported in Table 2 for each node, together with the total inhabitants of each province given by ISTAT data of 2019 [1] . We chose to associate all infected persons detected to the I compartment. This choice is justified by the ongoing screening policy. Indeed, during the first wave of this pandemic in Italy, tests to assess the presence of SARS-CoV-2 virus were performed almost only on patients with consistent symptoms and fever. We select µ 1 = 1, assuming no more than half of the actual highly symptomatic infected were detected at the very beginning of the outbreak of the pandemic.

According to values used in [20, 29] , we fix γ I = 1/14, γ A = 2γ I = 1/7, a = 1/3, considering these clinical parameters deterministic. We also assume the probability rate of developing severe symptoms σ = 1/12.5, as in [20, 36] . Since the first day of the simulation, the public was aware of the epidemic outbreak and recommendations such as washing hands often, not touching the face, avoiding handshakes and keeping distance had already been disseminated. Therefore, we initially set coefficients k I = k A = 30. The initial value of the transmission rate of asymptomatic/mildly infectious people is calibrated as the result of a least square problem, namely the L2 norm of the difference between the observed cumulative number of infected I(t) and the numerical evolution of [29] , a reduced contact rate is considered for the node of Lodi, with β A = 0.50. In the above fitting, the expected initial amount of exposed and asymptomatic/mildly symptomatic individuals is also estimated imposing an initial condition of one single exposed individual, obtaining e 0 ≈ 10.23 i 0 and a 0 ≈ 9.15 i 0 . In particular this permits to have a rough estimate of the presumed day of outbreak of the epidemic which, in our case, would appear to have started on January 14, 2020. With the above setup, we obtain an initial expected value of the basic reproduction number in the whole network (as from definition given in Section 2.2.1) R 0 = 3.6, which is in agreement with estimations reported in [20, 29, 55] .

On the other hand, the transmission rate of compartment I, β I , is considered a random parameter, given that, at each node, the initial amount of severe symptomatic subjects I is affected by uncertainty, as previously discussed. Therefore we impose:

Assuming that highly infectious subjects are mostly detected in the most optimistic scenario, being subsequently quarantined or hospitalized, the minimum value of the transmission rate of I is set β I,0 = 0.03 β A , as in [20, 29] . Furthermore, µ 2 = 0.06 −1 , indicating that the more the number of infected in I increases relative to the observed value i 0 , the more the transmission rate of this compartment tends to that of compartment A, proportionally to the error in I.

As a consequence, also initial conditions for compartments E, A and S are stochastic, depending on the initial amount of severe infectious at each location:

E(x, 0, z) = 10.23 I(x, 0, z) ,

A(x, 0, z) = 9.15 I(x, 0, z) ,

Finally, removed individuals are considered initially null everywhere in the network, R(x, 0, z) = 0.0, with all arcs assumed empty at the beginning of the simulation.

To model the escalation of lockdown restrictions, starting from March 9, 2020, initial day of the northern Italy lockdown, the transmission rate β A is reduced by 15% (with a consequent update also of β I,0 and therefore of β I ) and coefficients k I = k A = 60, due to the public being increasingly aware of the epidemic risks. Furthermore, the percentage of commuting individuals is reduced by 60% in the whole network according to mobility data tracked through mobile phones and made available by Google [6, 55] . After the additional restrictions in place as of March 22, 2020, β A is ulteriorly reduced by 10% (again with a consequent rearrangement also of β I,0 and β I ) and coefficients k are increased to k I = k A = 70.

Lodi (n 1 ) Bergamo (n 3 ) Figure 3 : Exponential decay of L ∞ norm of highly infectious I at node n 1 (left) and of asymptomatic/weakly symptomatic individuals at node n 3 (right) in the Lombardy network test, as the number of collocation points M p increases. Table 3 : Error estimates and empirical order of accuracy of expectation E and variance V of susceptible S and highly symptomatic individuals I in the Lombardy network test, at node n 1 (Lodi). Results obtained applying the sAP-IMEX Runge-Kutta FV Collocation method. M p indicates the number of points used for the stochastic Collocation method. 

As a first numerical test we analyze the numerical convergence of the stochastic Collocation method, discussed in details in Appendix B. We consider the above setting, but with a larger level of uncertainty, with µ 1 = 10, to emphasize the convergence rates. The numerical results evaluated with an increased number of collocation points M p are compared to a reference solution obtained using M p = 50, in terms of expectation and variance of the state variables. The expected exponential convergence is shown in Tables  3-4 for chosen state variables, respectively for nodes n 1 and n 3 (taken as representative nodes), where L 2 error norms and the related order of accuracy are presented. The result is highlighted by Fig. 3 , in which the spectral decay of the L ∞ norm is observed in terms of both expected value and variance as the number of collocation points increases.

In the following, we consider several scenarios regarding the spread of the COVID-19 epidemic in the Lombardy region based on the data and parameters determined in the previous section. In particular, we will consider the baseline scenario, corresponding to the actual spread of the epidemic observed from the data, and two hypothetical scenarios corresponding to the absence of long-range mobility and the absence of restrictions.

Numerical results of the Lombardy network test, in its baseline configuration (presented in details at the beginning of Section 3.1), obtained using the sAP IMEX Runge-Kutta FV Collocation method with M p = 6 points, are reported in Figs. 4-7. In Fig. 4 a first qualitative trend of the spatial and temporal spread of the epidemic is presented. Here it can be observed the heterogeneity of the diffusion of the virus, which has firstly mostly affected the city of Bergamo and only in a second moment prevailing in Milan. In Fig. 5 the expected evolution in time of the infected individuals, together with 95% confidence intervals, is shown for each node and for the whole Lombardy network, including exposed E, highly symptomatic subjects I and asymptomatic or weakly symptomatic people A. Each plot is also associated with the expected temporal evolution of the reproduction number R 0 (t), computed as described in Section 2.2.1.

One can see the capacity of the model to reproduce a very heterogeneous epidemic trend in the network analyzed, which is also reflected in the different ranges and patterns shown for the R 0 of each province. It can also be observed the agreement between the evolution of the reproduction number and the epidemic spread. In particular, it is confirmed the decline of the daily number of infected as R 0 reaches values below 1, as shown in the plots for Lodi, Bergamo and Cremona. On the other hand, the persistence of the virus in the complete network, and especially in Milan, is noticed until March 27, 2020 (last day of the simulation), where the reproduction number remains R 0 > 1, confirming the consistency of the definition proposed for R 0 (t).

As visible from Fig. 6 , the lower bound of the confidence band of the cumulative amount in time of highly symptomatic individuals is comparable with the observed data of the Civil Protection Department of Italy [2] . As expected, the mean value of the numerical result reports an higher amount of I, especially in Milan, the province most affected by the virus, due to the uncertainty of available data, which surely underestimate the real amount of infectious people, as discussed in Section 3.1.

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 Fig. 7 . Here, it can be noticed how much of the spread of COVID-19 has actually been lost from the data of the first outbreak in Lombardy and the impact that the presence of asymptomatic or undetected subjects has had on the epidemic evolution. By the end of the simulation, indeed, E[A + R A ]/E[I + A + R] = 0.92 in the network, indicating with R A the amount of recovered coming from the compartment A. This result is in line with WHO guidelines given during the first wave of COVID-19, stating that approximately 80% of the infected population was asymptomatic [5] . It is here remarked, indeed, that the compartment A in the proposed model includes not only the asymptomatic people but also the mildly symptomatic, which would approximately be 12% of the total cases. Remark 3.1. In the results here presented, transmission coefficients α i,j (x), which define the behavior of commuters in the network, are imposed as deterministic and constant in time, meaning that the amount of individuals exiting from (and entering in) each city is the same in each instant of the simulation. Clearly, this assumption leads to a simplification of the description of the phenomenon of commuting, generally characterized by a peak in the early morning hours and a drastic drop in the night hours. What we represent is indeed the mean commuting trend during the day, which tends to assume a stationary solution in time. However, a more realistic behavior, characterized for example by periodic sinusoidal functions, leads to a slightly oscillatory trend in the reported curves of each compartment and also in the reproduction number R 0 . By selecting α i,j (x) values that define an average trend over the day, allows us to avoid this misleading representation, while maintaining the consistency of the results.

To assess the effects of the mobility of commuters, a second scenario is investigated in which the whole population is not allowed to move out of the residence city. Numerical results of this test are reported in Fig. 8 for two representative provinces of the network: Lodi and Milan. It can be observed the different evolution of the spread of the epidemic comparing trends presented in Fig. 8 with the corresponding ones in Fig. 5 (for the first column of plots) and Fig. 7 (for the second column of plots). As a consequence of the absence of commuting people in the network, the spread of COVID-19 in Milan results consistently damped. The province of Milan, in fact, is the one with the highest number of daily incomes of commuting workers and students, followed by Bergamo and Brescia, as shown in Table 1 . On the other hand, Lodi is the city with the highest amount of daily outcomes of commuters (30% of its total population). Thus, preventing the population from leaving the province shows a slight increase in local infections with respect to the baseline scenario. This result is not intended to suggest that a constraint on the mobility of people would be disadvantageous in fighting the spread of the epidemic. In fact, one must consider that the contagions shown in these Figures are due to a still present interaction of people at the provincial level, which has not been reduced in any way with respect to the previous scenario. Similar conclusions are drawn also in [26] . These results are primarily intended to demonstrate how much the evolution of a pandemic can vary in response to changes in people's mobility.

Finally, results in the event that no restrictions of any kind were applied by the government are presented in Fig. 9 for the province of Milan, Bergamo and the complete Lombardy network. Comparing these results with the corresponding in Fig. 5 (for the first column of plots) and Fig. 7 (for the second column of plots), it can be seen immediately the fundamental role that government restrictions have played in containing the spread of the disease (notice the different y-axis scale in the Figures). In fact, the cumulative number of infected people across the region in this scenario is almost 1.8 times higher than in the baseline scenario at the end of the simulation. This value indicates that the restrictions imposed during the first wave of COVID-19 in Italy, and the consequent increasing risk awareness in the population, have contributed to attenuate the spread of the virus by 43%, which is in agreement with the result presented in [29] . Also the basic reproduction number R 0 shows a totally different evolution in time, being far from reaching the value 1 at the end of the simulation. Indeed, on March 27, 2020 the virus still persists in its propagation, with mean value R 0 = 1.50 when considering the whole network and, in particular, R 0 = 1.90 in the province of Milan.

In this paper, a stochastic transport model for the spread of an epidemic phenomenon described by a multi-population SEIAR compartmental dynamics on networks is presented. The starting point for the description of spatial motion and the interactions of the individuals has its roots in the kinetic theory of discrete velocity models [9, 12, 39] . The model is structured on spatial networks, where nodes identify the locations of interest (cities in this case) and arcs represent the set of major mobility paths (roads and railways). Individuals are divided into commuters, who move on a suburban scale, and non-commuters, who act on an urban scale. In this way, we avoid unrealistic mass migration effects in which the entire population moves indiscriminately through the network. Thanks to the hyperbolicity of the resulting system, the unphysical feature of instantaneous diffusive effects, which is typical of parabolic reaction-diffusion models proposed in literature, is removed. Nevertheless, we show that for small relaxation times and large characteristic speeds, in the diffusive regime, the proposed model recovers its parabolic nature.

Since the derivation and the consequent definition of the basic reproduction number is well-established in the literature for ODE epidemic models, we resort on these results to introduce a derivation of R 0 for spatial epidemic models in the case of no-inflow/outflow boundary conditions, validating the effectiveness of the resulting definition as an indicator of the viral growth rate. The model has been tested for the analysis of the emergence of COVID-19 in Italy, by simulating the propagation and evolution of the virus in the months of February and March 2020 in Lombardy. Indeed, it is in these early stages of the epidemic that uncertainty in data and transport dynamics played a key role. In order to study the effects of the uncertainties of the initial conditions and of the parameters involved in the model on the solution of the problem, a second-order stochastic asymptotic-preservative IMEX Runge-Kutta Finite Volume Collocation method has been used. It is shown that the proposed numerical scheme achieves spectral accuracy in the stochastic space by achieving an exponential convergence rate even in the case where uncertainty is present not only in the initial data but also in the nonlinear interaction terms.

Numerical results of the baseline scenario were compared with observed data made available by the Italian Civil Protection Department, demonstrating that the proposed model is suitable to adapt to real settings and applications and showing its ability on capturing the heterogeneity underlying epidemics. Also alternative scenarios have been evaluated, considering both a total lockdown of extra-urban mobility and a total absence of restrictions at governmental level.

The proposed methodology has the potential to be applied to epidemiological models structured in more compartments, such as those proposed in [20, 29, 30] . In this context, it was decided to keep the compartmentalization as simple as possible given the lack of more specific observed data at the provincial level. As a matter of fact, for the single provinces during the first epidemic wave, the Italian Civil Protection Department has made available only the cumulative trend of the subjects tested positive to COVID-19 [2] . There is no distinction between individuals simply quarantined, hospitalized or even in intensive care. Similarly, no structured data are publicly available at the level of provinces for individuals who have died.

Future developments foresee an extension of the model to include the age structure of the population, as this is an essential feature to correctly describe mobility flows (which mainly affect the younger part of the population) and consequently a more correct spread of viruses such as COVID-19 [7, 8] . Another interesting aspect, especially for analysis related to recent developments of the COVID-19 pandemic, is certainly the extension of the model to take into account the effects of the viral load of individuals, by analyzing the effects of the so called super-spreaders [38] , and the immunization of susceptible individuals through vaccination campaigns [31] . differential equations and applications", code 2017KKJP4X.

Recurring to the next-generation matrix (NGM) procedure, it is possible to compute the dominant eigenvalue of a positive linear operator, called next-generation operator, which permits to define the reproduction number when there are several compartments contributing to the spread of the infection [20, 24, 29, 30, 52, 54] . The generalization to space dependent models is non straightforward due to the nonlinear nature of the interactions. In the following we consider the situation where no-inflow/outflow boundary conditions are assumed.

Given the spatial model defined in (2.4), let us consider the following subsystem of balance equations of the infectious compartments including both commuters and noncommuters

Considering, a deterministic framework, in which the epidemic parameters only depend on the variable x, thanks to our assumptions, integration over Ω leads to

If we define the total densities

and the averaged epidemic operators,

, system (A.1) can be rewritten in the following ODE forṁ

Following the NGM approach [24, 29] , the Jacobian matrix of the above obtained ODE system results

which can be decomposed into the following transmission matrix T and transition matrix Σ, so that

The basic reproduction number R 0 is the spectral radius of the next-generation operator, K L = −TΣ −1 , and results composed by the sum of two components, deriving from the two compartments of the model contributing to the spread of the epidemic (I and A)

and therefore, recalling the previous definitions

Note that, the dependence of additional stochastic parameters will not modify the above reasoning leading to the same definition (A.5) including dependence from the uncertainty variable z = (z 1 , . . . , z d ) T ∈ R d .

Following a stochastic collocation approach [14, 56] , the solution of the stochastic problem (2.5) can be computed employing a generalized Polynomial Chaos (gPC) expansion [33, 34, 41] . Let us consider, for simplicity, a single source of uncertainty z ∈ R and that the probability density function (PDF) of this random input, ρ z : Γ → R + , is known. The approximated solution of the problem,

is the state vector of non-commuters, referring to system (2.3), can be expressed as a finite series of orthonormal polynomials in terms of the stochastic parameter where M is the number of terms of the truncated series and φ j (z) are orthonormal polynomials, with respect to the measure ρ z (z) dz. The expansion coefficientsQ j (x, t) is obtained asQ

Following the stochastic collocation method [34, 57] , the integrals for the expansion coefficients in Eq. (B.2) are replaced by a suitable quadrature U Mp characterized by the set {z n , w n } Mp n=1 , where z n is the n-th collocation point, w n is the corresponding weight and M p represents the number of quadrature points. For instance, for a stochastic parameter with a uniform distribution the quadrature is defined by the Gauss-Legendre weights and nodes; while for a random variable associated to a Gaussian PDF, we will rely on a Gauss-Hermite quadrature, which readŝ Similarly, other statistical quantities of interest can be computed [56] , like the variance

When more than one stochastic parameter is considered in the problem, z ∈ R d , the joint PDF ρ z (z) of the random vector composed by the random parameters (assuming independence of the variables) is given by

where ρ z,k is the PDF of the k-th variable. The extension of the collocation method then follows in a similar way [57] .

To compute the solution at each collocation point, an IMEX Runge-Kutta Finite Volume method for hyperbolic systems with multiscale relaxation is adopted [15, 16] . The IMEX discretization for the commuters' dynamics (2.5), coupled with the equations for the non commuter population (2. Here ∆t = t n+1 − t n is the time step size that follows the less restrictive between the standard hyperbolic Courant-Friedrichs-Levy condition, ∆t ≤ CFL ∆x max{λ i } , and the parabolic stability restriction, ∆t ≤ ν ∆x 2 max{D i } , given by the diffusive components of the system, with CFL and ν suitable stability constants [16] and ∆x = x i+ 1 2 − x i− 1 2 space grid size. In this work, we fix CFL = 0.9 and ν = max{D i } 2 .

The numerical scheme is characterized by two s×s matrices,à = (ã kj ), withã kj = 0 for j ≥ k, and A = (a kj ), with a kj = 0 for j > k, and by the weights vectorsb = (b 1 , ...,b s ) T , b = (b 1 , ..., b s ) T (with s identifying the number of the Runge-Kutta stages), which can be represented in the following Butcher notation [42] The distribution of matricesà and A permits to treat implicitly the stiff terms (hence those depending on the scaling parameters, τ i and λ i ) and explicitly all the rest. Moreover, if a proper set of matricesÃ, A and vectorsb, b are chosen, the AP property is satisfied, which means that the scheme maintains a consistent discretization of the asymptotic system in the diffusive (parabolic) regime (see [15, 16] ). For example, the second order GSA BPR(4,4,2) scheme proposed in [16] satisfies the AP property and is defined by the following double Butcher tableau (explicit on the left and implicit on the right) At each internal stage of the IMEX Runge-Kutta scheme (C.1), we apply a TVD Finite Volume discretization [14, 25] . To achieve second order accuracy also in space, while avoiding the occurrence of spurious oscillations, a classical minmod slope limiter has been adopted.

Finally, let us remark that, given the non-intrusive nature of the stochastic collocation method which only requires the evaluation of the solutions of the corresponding deterministic problem at each collocation point, the AP property of the deterministic IMEX scheme is preserved, leading to a stochastic asymptotic-preserving scheme [35] .

Dipartimento della Protezione Civile, Italia. COVID-19 epidemiological data in Italy

Q&A: Similarities and differences COVID-19 and influenza

Google COVID-19 Community Mobility Reports: Anonymization Process Description (version 1.1)

Modelling lockdown measures in epidemic outbreaks using selective socio-economic containment with uncertainty

Control with uncertain data of socially structured compartmental epidemic models

A unified multiscale vision of behavioral crowds

Spread of infectious diseases in a hyperbolic reaction-diffusion susceptible-infected-removed model

A multiscale model of virus pandemic: Heterogeneous interactive entities in a globally connected world

From the Boltzmann Equation to Discretized Kinetic Models

Simulation of an epidemic model with nonlinear crossdiffusion

Uncertainty quantification of viscoelastic parameters in arterial hemodynamics with the a-FSI blood flow model

Hyperbolic models for the spread of epidemics on networks: kinetic description and numerical methods

A unified IMEX Runge-Kutta approach for hyperbolic systems with multiscale relaxation

Modeling and simulating the spatial spread of an epidemic through multiscale kinetic transport equations

Numerical approximations of a traffic flow model on networks

A hyperbolic model of chemotaxis on a network: a numerical study

Effects of information-induced behavioural changes during the COVID-19 lockdowns: The case of Italy: COVID-19 lockdowns and behavioral change

A generalization of the Kermack-McKendrick deterministic epidemic model

An age and space structured SIR model describing the Covid-19 pandemic

A network model of Italy shows that intermittent regional strategies can alleviate the COVID-19 epidemic

On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations

On universal Osher-type schemes for general nonlinear hyperbolic conservation laws

Mobility restrictions for the control of epidemics: When do they work?

A feedback SIR (fSIR) model highlights advantages and limitations of infection-based social distancing

Systems of Conservation Equations with a Convex Extension

Spread and dynamics of the COVID-19 epidemic in Italy: Effects of emergency containment measures

Modelling the COVID-19 epidemic and implementation of populationwide interventions in Italy

Vaccination and SARS-CoV-2 variants: how much containment is still needed? A quantitative assessment

Efficient stochastic asymptotic-preserving implicitexplicit methods for transport equations with diffusive scalings and random inputs

Uncertainty Quantification for Hyperbolic and Kinetic Equations

Asymptotic-preserving methods for hyperbolic and transport equations with random inputs and diffusive scalings

Beyond just "flattening the curve": Optimal control of epidemics with purely non-pharmaceutical interventions

Non-linear incidence and stability of infectious disease models

Superspreading drives the COVID pandemic -and could help to tame it

Diffusive limit for finite velocity Boltzmann kinetic models

Monitoring Italian COVID-19 spread by a forced SEIRD model

An introduction to uncertainty quantification for kinetic equations and related problems

Implicit-explicit Runge-Kutta schemes and applications to hyperbolic systems with relaxation

SUIHTER: A new mathematical model for COVID-19. Application to the analysis of the second epidemic outbreak in Italy

Visualizing the invisible: The effect of asymptomatic transmission on the outbreak dynamics of COVID-19

Eight challenges for network epidemic models

Modeling blood flow in networks of viscoelastic vessels with the 1D augmented fluid-structure interaction system

Uncertainty quantification for systems of conservation laws

Five challenges for spatial epidemic models

A renewal equation model to assess roles and limitations of contact tracing for disease outbreak control

An updated estimation of the risk of transmission of the novel coronavirus (2019-nCov)

Estimation of the Transmission Risk of the 2019-nCoV and Its Implication for Public Health Interventions

Simulating the spread of COVID-19 via a spatially-resolved susceptible-exposed-infected-recovered-deceased (SEIRD) model with heterogeneous diffusion

Diffusion-reaction compartmental models formulated in a continuum mechanics framework: application to COVID-19, mathematical analysis, and numerical study

Using mobility to estimate the transmission intensity of COVID-19

Numerical Methods for Stochastic Computations -A Spectral Method Approach

High-Order Collocation Methods for Differential Equations with Random Inputs

This work was partially supported by MIUR (Ministero dell'Istruzione, dell'Università e della Ricerca) PRIN 2017, project "Innovative numerical methods for evolutionary partial