Hyperbolic models for the spread of epidemics on networks: kinetic description and numerical methods
HHyperbolic models for the spread of epidemics on networks:kinetic description and numerical methods
Giulia Bertaglia ∗ Lorenzo Pareschi † July 9, 2020
Abstract
We consider the development of hyperbolic transport models for the propagationin space of an epidemic phenomenon described by a classical compartmental dynam-ics. The model is based on a kinetic description at discrete velocities of the spatialmovement and interactions of a population of susceptible, infected and recovered in-dividuals. Thanks to this, the unphysical feature of instantaneous diffusive effects,which is typical of parabolic models, is removed. In particular, we formally show howsuch reaction-diffusion models are recovered in an appropriate diffusive limit. The ki-netic transport model is therefore considered within a spatial network, characterizingdifferent places such as villages, cities, countries, etc. The transmission conditionsin the nodes are analyzed and defined. Finally, the model is solved numerically onthe network through a finite-volume IMEX method able to maintain the consistencywith the diffusive limit without restrictions due to the scaling parameters. Severalnumerical tests for simple epidemic network structures are reported and confirm theability of the model to correctly describe the spread of an epidemic.
Keywords:
Kinetic equations, hyperbolic systems, spatial epidemic models, SIRmodel, network models, IMEX Runge-Kutta schemes, diffusive limit
Contents ∗ Department of Mathematics and Computer Science, University of Ferrara, Via Machiavelli 30, 44121Ferrara, Italy ([email protected]) † Department of Mathematics and Computer Science, University of Ferrara, Via Machiavelli 30, 44121Ferrara, Italy ([email protected]) a r X i v : . [ phy s i c s . s o c - ph ] J u l Numerical method 11
The recent COVID-19 pandemic has led to a strong interest from researchers aroundthe world in building and studying new epidemiological models capable of describing theprogress of the epidemic (see for example [1,7,12,14,16,17,22] and the references therein).Mathematical models can help on the prediction of the evolution of a pandemic and arean indispensable tool to support government decision-making on the control of infectiousdiseases. Together with computer simulations, they permit to build and test hypothesis,quantify the uncertainty related to random input parameter values, and estimate keyparameters from collected data. Usually, SIR-type compartmental models, inspired bythe seminal work of Kermack and McKendrick, are adopted [11, 19]. In the classicalsetting the population is assumed to be divided in three classes: susceptible (S), whocan contract the disease, infectious (I), who have contracted the disease and can transmitit, and removed or recovered (R), consisting of those with permanent infection-acquiredimmunity (eventually, even deceased). More complex models include additional classes ofindividuals, like the SEIR model, which considers also a subgroup of exposed subjects (E),individuating people in the latent period, who are infected but not yet infectious; or likethe MSEIR model, which includes also the class M of infants with passive immunity [19].The choice of which compartments to include in a model depends on the characteristicsof the infectious disease analyzed and the aim of the model. Models including severalcompartments have been especially designed to deal with the COVID-19 pandemic (seefor example [16, 17]).All these models represents the spread of the epidemic only concerning the temporalevolution of the disease among the population, but not taking into account spatial effects.In many cases, the concept of the average behavior of a large population is sufficient toprovide useful guidance on the development of the epidemic. However, the importance ofthe spatial component of many transmission systems is being increasingly recognized, espe-cially when there is a need to consider spatially heterogeneous interventions, as happenedfor COVID-19 [31]. Thus, to describe the spread of epidemics in a more detailed way, anddesign effective confinement strategies, several spatially extended models, based on sys-tems of parabolic reaction-diffusion equations, have been proposed [26, 32, 35]. However,2he parabolic character of these models leads the disease to propagate instantaneously inspace, with infinite speeds. This unphysical feature has been avoided through hyperbolicsystems based on relaxation approximations [5, 6] or nonlocal interactions [12].Furthermore, networks (or graphs) are extremely flexible tools for representing complexsystems of interacting components. Each component is represented by a node and eacharc (or edge) between nodes describes some sort of interaction between them. Becauseof their flexibility, networks have been used to model infection spread in different forms.Nodes can describe single individuals, groups of individuals (e.g. households, farms, cities)or locations to which individuals are connected [29, 31]. Links can represent infectiousattempts or transmission events or simply social relationships through which the infectioncan spread, movements of animals between farm, flight routes, streets [4, 24, 29]. We referto [16, 22] for recent applications of these ideas to the specific case of COVID-19.In this work, without claiming to provide an answer to the complex problem of spatialpropagation of an epidemic phenomenon, we try to provide a contribution in this directionby building models based on hyperbolic partial derivative equations that allow a more real-istic description of the dynamics of the individuals involved in the epidemic phenomenon.The model, inspired by two-velocity kinetic equations [23], is based on a kinetic descriptionof the diffusive spread of the epidemic governed by a SIR-type dynamics. In particular,it is shown that under a suitable scaling limit the model permits to recover the classicalreaction-diffusion SIR-type system [32, 35]. Although, to simplify our presentation, wehave focused on a simple SIR-type dynamic, the approach can be naturally extended tomore realistic compartmental descriptions including more populations.Subsequently, the hyperbolic system is considered on a spatial network where thenodes represent groups of individuals, typically villages, cities, regions or even countriesthat evolve through a standard SIR model. In addition, suitable transmission conditionsare derived at each arc-node interface. Note that here, unlike other network modelsbased on hyperbolic balance laws and kinetic equations, such as in chemotaxis and trafficflows [9, 10, 13, 30], the nodes themselves are evolving. The resulting system of equationsis then solved on each arc using a suitable IMEX finite-volume approach that permits toachieve uniform second order accuracy even in stiff regimes where the diffusive behaviordominates [8, 18, 27].The rest of the manuscript is organized as follows. In Section 2 we introduce themathematical model. First, we discuss the one-dimensional model in a bounded domain,its main properties, and formally derive the diffusion limit. Then, in Section 2.3, the modelis considered in a spatial network with appropriate transmission conditions. Section 3 isdedicated to the description of the adopted IMEX finite-volume scheme, which allowsto preserve the limiting parabolic behavior of the system. Several numerical examples,including networks with different arcs and nodes are presented in Section 4, to illustratethe ability of the model to correctly describe the spread of the epidemic. Some conclusionsare reported in the last Section. 3
Mathematical model
One of the standard compartmental models widely used in epidemiology is the so-called SIR model, in which the population is assumed to be divided in three classes: thesusceptible (S), the infectious (I), and the removed or recovered (R). Let us assume tohave a population with similar individuals, without prior immunity, which are uniformlymixed and that do not present behavioral changes during the epidemic and neglecting thevital dynamics represented by births and deaths because of the time scale considered. Inthe simplest one dimensional (1D) case, if we consider individuals moving in two oppositedirections (indicated by signs “+” and “-”), with velocities ± λ S for susceptible, ± λ I forinfectious and ± λ R for recovered, we can describe the spatio-temporal dynamics of thepopulation through the following SIR-type discrete velocity model ∂S + ∂t + λ S ∂S + ∂x = − f ( S + , I ) − τ S (cid:0) S + − S − (cid:1) ∂S − ∂t − λ S ∂S − ∂x = − f ( S − , I ) + 12 τ S (cid:0) S + − S − (cid:1) ∂I + ∂t + λ I ∂I + ∂x = f ( S + , I ) − γI + − τ I (cid:0) I + − I − (cid:1) ∂I − ∂t − λ I ∂I − ∂x = f ( S − , I ) − γI − + 12 τ I (cid:0) I + − I − (cid:1) ∂R + ∂t + λ R ∂R + ∂x = γI + − τ R (cid:0) R + − R − (cid:1) ∂R − ∂t − λ R ∂R − ∂x = γI − + 12 τ R (cid:0) R + − R − (cid:1) . (1)In the above system, individuals S ( x, t ), I ( x, t ) and R ( x, t ), representing, respectively,the number of susceptible, infectious and recovered, at location x ∈ Ω ⊆ R and time t > N = S + I + R ), are defined as S = S + + S − , I = I + + I − , R = R + + R − . The parameter γ represents the recovery rate, while the transmission of the infection isgoverned by an incidence function f ( S ± , I ) modeling the transmission of the disease. Weassume local interactions characterize the general incidence function f ( S ± , I ) = β S ± I p kI , (2)where the classic bilinear case corresponds to p = 1, k = 0, even though it has beenobserved that an incidence rate that increases more than linearly with respect to thenumber of infectious I can occur under certain circumstances [5, 11, 21]. The parameter β characterizes the contact rate [19], whereas the case k > γ and β can depend on the position x , as well as the positive relaxationtimes τ S , τ I , and τ R . It is important to notice that in this model the class of I includes alsothe asymptomatic population, which still acts like infected and can transmit the disease,even if never developing symptoms. Moreover, it is assumed that an individual, afterrecovering, always becomes immune.The model must then be supplemented with appropriate initial and boundary data.A more detailed discussion about boundary conditions is postponed to Section 2.3 wheresystem (1) is considered within a spatial network. A structural property of (1) is thatthe first four equations are independent of the last two. Once I ± are known, the explicitforms of R ± can be determined directly by solving the last two equations.The standard threshold of epidemic models is the well-known reproduction number R , which defines the average number of secondary infections produced when one infectedindividual is introduced into a host population in which everyone is susceptible [19]. Thisnumber determines when an infection can invade and persist in a new host population.For many deterministic infectious disease models, an infection begin in a fully susceptiblepopulation if and only if R >
1. In the kinetic SIR-type model (5) the reproductionnumber is characterized locally by the ratio β ( x ) /γ ( x ). More precisely, assuming noinflow/outflow boundary conditions in Ω, summing up the third and fourth equations in(1) and integrating over space we have ∂∂t (cid:90) Ω I ( x, t ) dx = (cid:90) Ω f ( S, I ) dx − (cid:90) Ω γ ( x ) I ( x, t ) dx ≥ R ( t ) = (cid:82) Ω f ( S, I ) dx (cid:82) Ω γ ( x ) I ( x, t ) dx ≥ . (3)The above quantity therefore, defines the basic reproduction number for system (1) de-scribing the space averaged instantaneous variation of the number of infective individualsat time t > ∂∂t (cid:90) Ω ( S ( x, t ) + I ( x, t ) + R ( x, t )) dx = 0 . If we now introduce the fluxes, defined by J S = λ S (cid:0) S + − S − (cid:1) , J I = λ I (cid:0) I + − I − (cid:1) , J R = λ R (cid:0) R + − R − (cid:1) , (4)we obtain an hyperbolic model equivalent to (1), but presenting in the following formula-5ion a macroscopic description of the propagation of the epidemic at finite speeds ∂S∂t + ∂J S ∂x = − f ( S, I ) ∂I∂t + ∂J I ∂x = f ( S, I ) − γI∂R∂t + ∂J R ∂x = γI∂J S ∂t + λ S ∂S∂x = − f ( J S , I ) − J S τ S ∂J I ∂t + λ I ∂I∂x = λ I λ S f ( J S , I ) − γJ I − J I τ I ∂J R ∂t + λ R ∂R∂x = λ R λ I γJ I − J R τ R . (5)Formally, the above system (5) is a system of balance laws which can be rewritten in acompact form as ∂ t U + ∂ x V = F ( U ) ∂ t V + Λ ∂ x U = G ( U , V ) + H ( V ) , (6)in which U = SIR , V = J S J I J R , Λ = λ S λ I
00 0 λ R , F ( U ) = − f ( S, I ) f ( S, I ) − γIγI , G ( U , V ) = − f ( J S , I ) λ I λ S f ( J S , I ) − γJ Iλ R λ I γJ R , H ( V ) = − J S /τ S J I /τ I J R /τ R . It is immediate to see that system (6) is symmetric hyperbolic in the sense of Friedrichs-Lax [15] with real finite characteristic velocities (eigenvalues) λ S , λ I , λ R , and a completeset of linearly independent eigenvectors. All the eigenvectors are associated with genuinelynon-linear fields, defining shocks and rarefactions, and the Riemann invariants of thesystem correspond to the kinetic variables S ± = 12 (cid:18) S ± J S λ S (cid:19) , I ± = 12 (cid:18) I ± J I λ I (cid:19) , R ± = 12 (cid:18) R ± J R λ R (cid:19) . (7)Furthermore, the symmetric hyperbolicity of the system guarantees the well-posedness ofthe Cauchy problem (i.e. existence, uniqueness and continuous dependence of the solutionson the initial data) for suitable smooth initial data (see [2, 25]). Remark . The hyperbolic model (5) differs from the models proposed in [5, 6] where asimple relaxation approximation was used to avoid the paradox of infinite propagationspeed. In particular, in [6] a different asymptotic behavior was considered, correspondingto the case τ S , τ I , τ R → ∞ , in which the dynamics can be reduced to a pair of coupledsecond order wave equations for S and I eliminating the fluxes J S and J I . It is easy toverify that this simplification is not possible in (5) due to the presence of the incidencefunction and recovery rate on the r.h.s. of the flux equations.6 .2 Diffusion limit Let us now consider the behavior of the model in diffusive regimes. To this aim, let usintroduce the diffusion coefficients D S = λ S τ S , D I = λ I τ I , D R = λ R τ R , (8)which characterize the diffusive transport mechanism of susceptible, infectious and recov-ered, respectively. The diffusion limit of the system is formally recovered letting the relax-ation times τ S , τ I , τ R →
0, and simultaneously the characteristic speeds λ S , λ I , λ R → ∞ ,while keeping the diffusion coefficients (8) finite. Under this scaling, from the last threeequations in (5) we get a proportionality relation between the fluxes and the spatial deriva-tives J S = − D S ∂S∂x , J I = − D I ∂I∂x , J R = − D R ∂R∂x , (9)which inserted into the first three equations, lead to the parabolic reaction-diffusion system[26] ∂S∂t = − f ( S, I ) + ∂∂x (cid:18) D S ∂S∂x (cid:19) ∂I∂t = f ( S, I ) − γI + ∂∂x (cid:18) D I ∂I∂x (cid:19) ∂R∂t = γI + ∂∂x (cid:18) D R ∂R∂x (cid:19) . (10)We refer to [23] concerning rigorous results on the diffusion limit of discrete velocity ki-netic models of type (1). Here we limit ourselves to note that these results cannot bestraightforwardly extended to our case, since the epidemic reaction terms destroy the typ-ical monotone behavior of solutions for such discrete kinetic systems. Although interestingthis aspect goes beyond the scopes of the present manuscript and will be the subject offuture investigations. A network or a connected graph G = ( N , A ) is composed of two finite sets: a set of N nodes (or vertices) N and a set of A arcs (or edges) A , such that an arc connects a pair ofnodes [30]. Arcs are bidirectional, as the graph is non-oriented, but an artificial orientationneeds to be fixed in order to define a sign for the velocities and therefore the fluxes. At anode n ∈ N , two different types of edges are present: incoming and outgoing [9, 10]. a a L L n n n Figure 1: Schematic representation of a simple network composed of 3 nodes ( n , n , n )and 2 arcs ( a , a ). 7 tt+ Δ t t+ Δ t Q n u + L,i u + n u − Q L,i Q Q *L,i Q *0,j Q *L,n Q *0,n interface L interface 0 u − n Figure 2: Layout of the Riemann problems solved to calculate the states at time t + ∆ t ,imposing the transmission conditions at the left ( L ) and right (0) boundary of a genericnode n ∈ N at which, respectively, a i ∈ A , i = 1 , . . . , N a,L and a j ∈ A , j = 1 , . . . , N a, arcs converge.For instance, in the simple network shown in Fig. 1, at node n there is an incomingarc ( a ) and an outgoing one ( a ). The A arcs of the network are parametrized as intervals a i = [0 , L i ] , i = 1 , . . . , A . For an incoming arc, L i is the abscissa of the node, whereas foran outgoing arc the same abscissa is 0 (see Fig. 1).In the model here proposed, the nodes of the network identify a particular location ofinterest (a municipality, a province, a region, a nation...), while each arc represents thewhole set of paths connecting each location to the others. Hence, each node is consideredwith its own initial population and a localized social dynamics, which can be influencedby those individuals moving from and to the different locations considered in the network,traveling through the arcs. This produces a network model in which on each arc system(5) is solved, while on each node the standard SIR compartmental model is evaluated [19].To ensure the correct coupling between nodes and arcs, specific transmission conditionsmust be prescribed at nodes. 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, outgoing arcs, we need to con-sider (1 + N a,L + N a, ) states at time t , Q L,i , Q n and Q ,j , separated by the interface ofincoming arcs (L) and the interface of outgoing arcs (0), as shown in Fig. 2. If variablesare discontinuous across these interfaces, (1 + N a,L ) new states originate at interface L, Q ∗ L,i and Q ∗ L,n , and (1 + N a, ) new states originate at interface 0, Q ∗ ,n and Q ∗ ,j , at time t + ∆ t [30]. To compute them, we sought for the solution of (2 + N a L + N a, ) Riemannproblems, recurring to the Riemann Invariants (or kinetic variables) of the system, definedin Eqs. (7), and to the principle of conservation of fluxes at interfaces [9, 10].For each one of the three compartments of individuals of the SIR-type model discussedin § u the number of individuals of thecompartment, with v the corresponding analytical flux, with λ its characteristic velocity,and with u ± the Riemann Invariants. To impose the transmission conditions at interface8, we need to solve the following system (for each compartment of the model) u ∗ L,i + v ∗ L,i λ i = u + L,i (11a) u ∗ L,i − v ∗ L,i λ i = N a,L (cid:88) k =1 α i,k u + L,k + α i,n u − n (11b) u ∗ L,n − v ∗ L,n λ n = u − n (11c) u ∗ L,n + v ∗ L,n λ n = N a,L (cid:88) k =1 α n,k u + L,k + α n,n u − n ; (11d)while to impose the transmission conditions at interface 0, we need to solve the followingsystem u ∗ ,n + v ∗ ,n λ n = u + n (12a) u ∗ ,n − v ∗ ,n λ n = N a, (cid:88) k =1 α n,k u − ,k + α n,n u + n (12b) u ∗ ,j − v ∗ ,j λ j = u − ,j (12c) u ∗ ,j + v ∗ ,j λ j = N a, (cid:88) k =1 α j,k u − ,k + α j,n u + n . (12d)Constants α i,j ∈ [0 ,
1] are the transmission coefficients and represent the probability thatan individual at a generic interface decides to move across that interface from location i to location j , also including the turnabout on the same location. It is worth to noticethat the condition differs when considering an incoming or an outgoing flux, due to theartificial orientation that has been fixed.Furthermore, to guarantee the conservation of fluxes at the interface, therefore ensuringthat the global mass of the system is conserved along the time, the following must hold[9, 10]: v ∗ L,n = N a,L (cid:88) i =1 v ∗ L,i , v ∗ ,n = N a, (cid:88) j =1 v ∗ ,j . (13)To respect these conditions, it is enough to impose at interface L λ i = N a,L (cid:88) k =1 α k,i λ k + α n,i λ n , λ n = N a,L (cid:88) k =1 α k,n λ k + α n,n λ n , (14)and at interface 0: λ n = N a, (cid:88) k =1 α k,n λ k + α n,n λ n , λ j = N a, (cid:88) k =1 α k,j λ k + α n,j λ n . (15)9t is straightforward to analytically solve system (11), obtaining u ∗ L,i = 12 u + L,i + N a,L (cid:88) k =1 α i,k u + L,k + α i,n u − n (16a) v ∗ L,i = λ i u + L,i − N a,L (cid:88) k =1 α i,k u + L,k − α i,n u − n (16b) u ∗ L,n = 12 u − n + N a,L (cid:88) k =1 α n,k u + L,k + α n,n u − n (16c) v ∗ L,n = − λ n u − n − N a,L (cid:88) k =1 α n,k u + L,k − α n,n u − n . (16d)The same applies for system (12), which results u ∗ ,n = 12 u + n + N a, (cid:88) k =1 α n,k u − ,k + α n,n u + n (17a) v ∗ ,n = λ n u + n − N a, (cid:88) k =1 α n,k u − ,k − α n,n u + n (17b) u ∗ ,j = 12 u − ,j + N a, (cid:88) k =1 α j,k u − ,k + α j,n u + n (17c) v ∗ ,j = − λ j u − ,j − N a, (cid:88) k =1 α j,k u − ,k − α j,n u + n . (17d) Nodes located at the extremes of the domain are without any incoming arc, when con-cerning inlet boundaries, or without any outgoing arc, when concerning outlet boundaries.At these nodes, in order to ensure that there are no individuals entering to or exiting fromthe network (hence considering a population which remains constant in time), we justrecover the standard zero-flux boundary conditions [30], which consists on imposing atinlet nodes v ∗ L,n = 0 , u ∗ L,n = u n − v n λ n ; (18)while at outlet nodes: v ∗ ,n = 0 , u ∗ ,n = u n + v n λ n . (19)10 Numerical method
Numerical methods for hyperbolic balance laws and kinetic equations in the diffusivelimit has a long history (see [8, 18, 20, 27] and the references therein). Here, to solvesystem (5), the Implicit-Explicit (IMEX) Runge-Kutta methods recently proposed in [8]for hyperbolic systems with multiscale relaxation terms are considered.
IMEX Runge-Kutta schemes can be easily represented by a double tableau (expliciton the left, implicit on the right) in the usual Butcher notation [8, 28]˜ c ˜ A ˜ b T c Ab T Matrices ˜ A = (˜ a kj ), with ˜ a kj = 0 for j ≥ k , and A = ( a kj ) are s × s matrices, with s numberof Runge-Kutta stages. It is always preferable in terms of computational efficiency to dealwith diagonally implicit Runge-Kutta (DIRK) schemes, which ensure that the explicitpart of the scheme term is always evaluated explicitly [3], hence a kj = 0 for j > k . Thecoefficient vectors ˜ c and c are given by˜ c k = k − (cid:88) j =1 ˜ a kj , c k = k (cid:88) j =1 a kj . and vectors ˜ b = (˜ b , ..., ˜ b s ) T and b = ( b , ..., b s ) T are the quadrature weights that permitto combine the internal Runge-Kutta stages. Furthermore, referring to [8], if the followingrelations hold, the method is said to be globally stiffly accurate (GSA) Definition 1.
An IMEX-RK method is said to be globally stiffly accurate (GSA) if notonly the corresponding diagonally implicit Runge-Kutta (DIRK) method is stiffly accurate,namely a kj = b j , j = 1 , . . . , s, but also the explicit method satisfies˜ a kj = ˜ b j , j = 1 , . . . , s − . It is worth to notice that this definition states also that the numerical solution of aGSA IMEX-RK scheme coincides exactly with the last internal stage of the scheme.In system (5), the scaling depends on the relaxation times. Indeed, these relaxationterms modify the nature of the behavior of the solution, which can result either hyperbolicor parabolic. Standard IMEX Runge-Kutta methods for hyperbolic systems with relax-ation terms loose their efficiency and a different approach must be adopted to guaranteeasymptotic preservation (AP) in stiff regimes (i.e. the consistency of the scheme withthe equilibrium system is guaranteed and the order of accuracy is maintained in the stifflimit). 11ollowing [8], the IMEX Runge-Kutta approach applied to system (6) consists in com-puting the internal stages U ( k ) = U n − ∆ t k (cid:88) j =1 a kj ∂ x V ( j ) + ∆ t k − (cid:88) j =1 ˜ a kj F (cid:16) U ( j ) (cid:17) V ( k ) = V n − ∆ t k − (cid:88) j =1 ˜ a kj (cid:16) Λ ∂ x U ( j ) − G (cid:16) U ( j ) , V ( j ) (cid:17)(cid:17) + ∆ t k (cid:88) j =1 a kj H (cid:16) V ( j ) (cid:17) , (20)followed by the numerical solution U n +1 = U n − ∆ t s (cid:88) k =1 b k ∂ x V ( k ) + ∆ t s (cid:88) k =1 ˜ b k F (cid:16) U ( k ) (cid:17) V n +1 = V n − ∆ t s (cid:88) k =1 ˜ b k (cid:16) Λ ∂ x U ( k ) − G (cid:16) U ( k ) , V ( k ) (cid:17)(cid:17) + ∆ t s (cid:88) k =1 b k H (cid:16) V ( k ) (cid:17) , (21)with a time step size ∆ t = t n +1 − t n that follows the less restrictive between the standardhyperbolic CFL condition, ∆ t ≤ ∆ xmax { λ S ,λ I ,λ R } , and the parabolic stability restriction,∆ t ∼ ∆ x , given by the diffusive components of the system [8]. The scheme (20)-(21) permits to treat implicitly the stiff terms and explicitly all therest, maintaining a consistent discretization of the limit system in the stiff regime, rep-resented by system (10), i.e. the AP property [8]. To verify the AP property, let usdenote by u ( j ) h , v ( j ) h , f ( j ) h , g ( j ) h , u nh , and v nh , h = 1 , ,
3, the components of U ( j ) , V ( j ) , F ( U ( j ) ), G ( U ( j ) , V ( j ) ), U n , and V n , respectively. Then, the IMEX scheme (20)-(21) canbe rewritten highlighting the scale parameters as follows u h = u nh e − ∆ tA∂ x v h + ∆ t ˜ A f h v h = v nh e − ∆ t ˜ A (cid:18) D h τ h ∂ x u h − g h (cid:19) − ∆ tτ h A v h , (22)where u h = ( u (1) h , . . . , u ( s ) h ) T , v h = ( v (1) h , . . . , v ( s ) h ) T , f h = ( f (1) h , . . . , f ( s ) h ) T , g h = ( g (1) h , . . . , g ( s ) h ) T , e = (1 , . . . , T ∈ R s , and we use notations τ h and D h to denote the relaxation times anddiffusion constant of each compartment. The final solution therefore reads u n +1 h = u nh − ∆ tb T ∂ x v h + ∆ t ˜ b T f h v n +1 h = v nh − ∆ t ˜ b T (cid:18) D h τ h ∂ x u h − g h (cid:19) − ∆ tτ h b T v h . (23)From the second equation in (22) we obtain v h = (cid:16) τ h ∆ t I d + A (cid:17) − (cid:16) τ h ∆ t v nh e − ˜ A ( D h ∂ x u h − τ h g h ) (cid:17) , (24)12hich substituted into the first equation of (22) gives u h = u nh e − ∆ tA∂ x (cid:16) τ h ∆ t I d + A (cid:17) − (cid:16) τ h ∆ t v nh e + ˜ Aτ h g h (cid:17) + ∆ tA (cid:16) τ h ∆ t I d + A (cid:17) − ∂ x ( D h ∂ x u h ) + ∆ t ˜ A f h . (25)Now, as τ → u h = u nh e + ∆ t ˜ A∂ x ( D h ∂ x u h ) + ˜ A f h , (26)and thus, the internal stages correspond to the stages of the explicit scheme applied tothe reaction-diffusion system (10).However, this is not enough to guarantee that the scheme satisfies the AP property,since we need to verify the same consistency property on the final numerical solution.To this aim, let us rewrite (24) as v h = τ h ∆ t A − v nh e + A − (cid:16) I d − τ h ∆ t A − (cid:17) ˜ AD h ∂ x u h + τ h A − ˜ A g h + O (cid:0) τ h (cid:1) , which substituted into the second equation of (23) leads to v n +1 h = (cid:0) − b T A − e (cid:1) v nh − (cid:16) b T A − ˜ A − ˜ b T (cid:17) ∆ tτ h D h ∂ x u h − b T A − ˜ AD h ∂ x u h − ∆ t (cid:16) b T A − ˜ A − ˜ b T (cid:17) g h + O ( τ h ) . In order to pass to the limit τ h →
0, we must require that b T A − ˜ A − ˜ b T = 0, which issatisfied if the IMEX scheme is GSA. Indeed, in GSA methods, b T = e Ts A and ˜ b T = e Ts ˜ A ,therefore b T A − ˜ A − ˜ b T = e Ts ˜ A − ˜ b T = 0. Thus, in the limit τ h → v n +1 h = − b T A − ˜ AD h ∂ x u h , and the resulting numerical solution for u n +1 h satisfies u n +1 h = u nh + ∆ t ˜ b T ∂ x ( D h ∂ x u h ) + ˜ b T f h , (28)which, combined with (26), proves the AP-property of the IMEX scheme with respect tothe reaction-diffusion system (10).In the present work, the GSA BPR(4,4,2) scheme presented in [8] is chosen, character-ized by s = 4 stages for the implicit part, 4 stages for the explicit part and 2nd order ofaccuracy, which can be defined by the following tableau (explicit on the left and impliciton the right)0 0 0 0 0 01/4 1/4 0 0 0 01/4 13/4 -3 0 0 03/4 1/4 0 1/2 0 01 0 1/3 1/6 1/2 00 1/3 1/6 1/2 0 0 0 0 0 0 01/4 0 1/4 0 0 01/4 0 0 1/4 0 03/4 0 1/24 11/24 1/4 01 0 11/24 1/6 1/8 1/40 11/24 1/6 1/8 1/4 (29)13o obtain a fully discrete scheme, we consider a finite volume method for the spatialdiscretization, and uniform grid with mesh spacing ∆ x = x i + − x i − . For each internalstep of the IMEX scheme, numerical fluxes are evaluated following the upwind Godunovflux [33]. Boundary-extrapolated values on the two sides of the interface within cell i arecomputed by piecewise linear reconstruction, through the second-order central differencescheme. We point out that the stages of the IMEX scheme are implemented in the form(25), where the second order derivative is discretized directly using a second order accuratecentral scheme. This avoids the creation of large non compact stencils in the discretizationof the resulting second order terms and the possible nonlinear effects due to the presenceof limiters in the discretization of the hyperbolic terms. In a similar way one computesalso the final numerical solution. Remark . In a similar way, it is possible to construct an IMEX-scheme which avoids theparabolic stability condition ∆ t = O (∆ x ) leading to an IMEX solver for the reaction-diffusion system [8]. This is obtained by modifying the equations for V in the IMEX-scheme (20)-(21) as follows V ( k ) = V n + ∆ t k − (cid:88) j =1 ˜ a kj G (cid:16) U ( j ) , V ( j ) (cid:17) − ∆ t k (cid:88) j =1 a kj (cid:16) Λ ∂ x U ( j ) − H (cid:16) V ( j ) (cid:17)(cid:17) , V n +1 = V n + ∆ t s (cid:88) k =1 ˜ b k G (cid:16) U ( k ) , V ( k ) (cid:17) − ∆ t s (cid:88) k =1 b k (cid:16) Λ ∂ x U ( k ) − H (cid:16) V ( k ) (cid:17)(cid:17) , (30)which, as τ → u h = u nh e + ∆ tA∂ x ( D h ∂ x u h ) + ˜ A f h u n +1 h = u nh + ∆ tb T ∂ x ( D h ∂ x u h ) + ˜ b T f h . (31)The derivation, under the GSA assumption, follows the same lines of the explicit case andwe omit the details for brevity.Finally, it is worth to notice that in the IMEX-schemes here proposed, each Runge-Kutta step can be expressed in explicit form, avoiding the adoption of iterative procedures(e.g. Newton-Raphson method or similar ones) for the evaluation of the implicit part,leading to a consistent reduction of the computational cost. In this Section, we present some numerical results to support the validity of the pro-posed model, both in the simple 1D configuration and in the network characterization.The accuracy of the scheme is verified for different values of the relaxation times, includingthe stiff regime of the purely diffusive system. Furthermore, we analyze the behavior of themodel concerning spatially heterogeneous environments, taking into account a spatiallyvariable contact rate, with respect to two different scenarios: R < R >
1. Twotest cases are then designed and simulated to observe the spread of infectious diseaseswith respect to the mobility of individuals on networks. In all test cases we assume p = 1in (2). 14able 1: L error norms and empirical order of accuracy of the variable S resulting fromaccuracy analysis performed choosing different values of relaxation times (and character-istic velocities). N x indicates the number of cells in the computational domain. τ = 1 . τ = 10 − τ = 10 − N x L O (cid:0) L (cid:1) L O (cid:0) L (cid:1) L O (cid:0) L (cid:1)
15 7.8183e-02 6.3566e-02 5.7368e-0245 1.1983e-02 1.7072 6.8991e-03 2.0214 6.1980e-03 2.0255135 1.4804e-03 1.9035 7.5256e-04 2.0168 6.7867e-04 2.0133405 1.5235e-04 2.0698 7.4838e-05 2.1010 6.7836e-05 2.0963Table 2: L error norms and empirical order of accuracy of the variable I resulting from ac-curacy analysis performed choosing different values of relaxation times (and characteristicvelocities). N x indicates the number of cells in the computational domain. τ = 1 . τ = 10 − τ = 10 − N x L O (cid:0) L (cid:1) L O (cid:0) L (cid:1) L O (cid:0) L (cid:1)
15 7.9779e-02 4.6081e-02 3.9811e-0245 1.1861e-02 1.7349 4.9662e-03 2.0278 4.0831e-03 2.0729135 1.4852e-03 1.8912 5.4088e-04 2.0182 4.4396e-04 2.0197405 1.5351e-04 2.0658 5.4262e-05 2.0930 4.4342e-05 2.0970Table 3: L error norms and empirical order of accuracy of the variable J S resulting fromaccuracy analysis performed choosing different values of relaxation times (and character-istic velocities). N x indicates the number of cells in the computational domain. τ = 1 . τ = 10 − τ = 10 − N x L O (cid:0) L (cid:1) L O (cid:0) L (cid:1) L O (cid:0) L (cid:1)
15 6.7762e-02 7.8869e-02 8.5108e-0245 1.0832e-02 1.6689 9.1743e-03 1.9583 8.9856e-03 2.0465135 1.2654e-03 1.9544 1.0350e-03 1.9861 8.4277e-04 2.1542405 1.2807e-04 2.0849 1.1343e-04 2.0126 9.0422e-05 2.0318Table 4: L error norms and empirical order of accuracy of the variable J I resulting fromaccuracy analysis performed choosing different values of relaxation times (and character-istic velocities). N x indicates the number of cells in the computational domain. τ = 1 . τ = 10 − τ = 10 − N x L O (cid:0) L (cid:1) L O (cid:0) L (cid:1) L O (cid:0) L (cid:1)
15 1.1164e-01 6.2077e-02 6.4906e-0245 1.7386e-02 1.6926 8.3184e-03 1.8295 6.4301e-03 2.1044135 2.1539e-03 1.9009 9.6011e-04 1.9654 5.9580e-04 2.1653405 2.2147e-04 2.0706 1.0974e-04 1.9742 6.3867e-05 2.032715 .1 Accuracy analysis
To verify the order of accuracy of the scheme, even in stiff regimes, an accuracy analysisis conducted considering periodic boundary conditions and the following initial conditions S ( x,
0) = 0 . (cid:18) πxL (cid:19) , I ( x,
0) = 1 − S ( x, , R ( x,
0) = 0 . , with x ∈ [ −
1; 1], hence a domain length L = 2, and no initial fluxes. In these tests, thecontact rate is β = 10 . γ = 4 .
0, with k = 0 as the classicalbilinear case. In each simulation, relaxation times and characteristic velocities are fixedequal for all the compartments of the population λ S = λ I = λ R = λ , taking into accountthree different cases, corresponding to a hyperbolic system ( τ = 1 . λ = 1 . τ = 10 − , λ = 10 ) and a purely diffusive system ( τ = 10 − , λ = 10 ).The stability condition is satisfied imposing ∆ t = ∆ x max { CFL /λ ; ∆ x/ } . The finaltime of the simulations is t end = 0 .
1. The refinement of the computational grid is madewith a factor 3, to work with embedded grids, and the time step ∆ t decreases with ∆ x accordingly to the stability condition. Values obtained with N x = 1215 cells are taken asreference solutions.Results are presented for variables S , I , J S and J I in Tables 1, 2, 3 and 4, respectively,in terms of L error norms and related order of accuracy, evaluated as O (cid:0) L (cid:1) = log (cid:18) || E ∆ x |||| E ∆ x/ || (cid:19) , with || E ∆ x || the relative error computed with grid size ∆ x . We observe that the second-order of accuracy of the method is satisfied uniformly in all regimes. It is worth tounderline that the expected accuracy is maintained in the limit diffusive regime for all thevariables, even for the fluxes, thanks to the GSA property of the scheme and because thefollowing additional conditions are satisfied by the IMEX-scheme (29) (see [8]) b T A − ˜ Ae = 1 , b T A − ˜ Ac = 1 , b T A − ˜ AAc = 1 / , b T A − ˜ A ˜ Ac = 1 / . (32) Following [34], we analyze the behavior of the model concerning spatially heterogeneousenvironments, taking into account a spatially variable contact rate β ( x ) = ˆ β (cid:18) .
05 sin 13 πx (cid:19) . Initial conditions are imposed as follows S ( x,
0) = 1 − I ( x, , I ( x,
0) = 0 . e − ( x − , R ( x,
0) = 0 . , with fluxes J S ( x,
0) = J I ( x,
0) = J R ( x,
0) = 0 . k = 0. Sim-ulations are performed in two different scenarios. In the first one an overall value in the16 a) (b) (c)(d) (e) -3 (f) Figure 3: Numerical results of the spatially heterogeneous test case with hyperbolic con-figuration of relaxation times and characteristic velocities ( τ = 1 . , λ = 1 .
0) and repro-duction number R <
1. Time and spatial evolution of S presented in (a), (b) and (c);evolution of I shown in (d), (e) and (f).domain of the initial reproduction number R = 0 . < β = 8 . γ = 4 .
0, corresponding to an infection which is not able to start spreading. In thesecond scenario, the initial reproduction number is R = 1 . >
1, given by the choiceˆ β = 11 . γ = 4 .
0, which identifies an infection that can persist in a new host popu-lation. Moreover, for each scenario, two different sets of relaxation times are considered,to concern both the hyperbolic and the parabolic limit of the system of equations. In thehyperbolic configuration, the relaxation times of all the compartments of individuals are τ = 1 .
0, with the square of the characteristic velocities λ = 1 .
0; while in the parabolicconfiguration τ = 10 − and λ = 10 . Simulations are run with N x = 150 cells in a domainhaving length L = 20, fixing CFL = 0 .
9, with final time t end = 10.Numerical results of each one of the two tests performed for the scenario representingan infectious disease characterized by R < τ = 1 . , λ =1 . τ = 10 − , λ = 10 , different dynamics of theinfectious spread are noticed, which represents the tendency of the system towards moreand more diffuse behavior as relaxation times reach values close to zero and characteristicspeeds close to infinity.In Figs. 5 and 6, numerical results for the scenario with R slightly > a) (b) (c)(d) (e) -3 (f) Figure 4: Numerical results of the spatially heterogeneous test case with parabolic config-uration of relaxation times and characteristic velocities ( τ = 10 − , λ = 10 ) and repro-duction number R <
1. Time and spatial evolution of S presented in (a), (b) and (c);evolution of I shown in (d), (e) and (f).Here, a temporary persistence of the infectious can be noticed, with oscillations thatreflect the sinusoidal form of the spatially variable contact rate. In this case, differencesof the dynamics of the epidemics in the two configurations of the relaxation times aremore accentuated. In particular, observing the evolution of susceptible individuals, it canbe seen that in the purely diffusive case the amount of susceptible tends to a much lowerequilibrium value than in the hyperbolic case, with almost all the individuals of the systeminfected by the disease.Table 5: Transmission coefficients used for the network test with 3 nodes connected by 2arcs in series (Fig. 7), given for each node-arc interface I . I ( n , a ) α n ,a = 1 . α n ,n = 1 . α a ,a = 0 . α a ,n = 0 . I ( a , n ) α a ,n = 0 . α a ,a = 0 . α n ,n = 0 . α n ,a = 0 . I ( n , a ) α n ,a = 0 . α n ,n = 0 . α a ,a = 0 . α a ,n = 0 . I ( a , n ) α a ,n = 0 . α a ,a = 0 . α n ,n = 1 . α n ,a = 1 . a) (b) (c)(d) (e) -3 (f) Figure 5: Numerical results of the spatially heterogeneous test case with hyperbolic con-figuration of relaxation times and characteristic velocities ( τ = 1 . , λ = 1 .
0) and repro-duction number R >
1. Time and spatial evolution of S presented in (a), (b) and (c);evolution of I shown in (d), (e) and (f). To assess the effects of the mobility of individuals on networks with respect to thespread of an infectious disease characterized by β = 10 . γ = 4 .
0, three tests areperformed concerning different simple networks.Table 6: Transmission coefficients used for the network test with 4 nodes connected by 3arcs (Fig. 9), given for each node-arc interface I . I ( n , a ) α n ,a = 1 . α n ,n = 0 . α a ,a = 0 . α a ,n = 0 . I ( a , n ) α a ,n = 0 . α a ,a = 0 . α n ,n = 0 . α n ,a = 0 . I ( n , a , a ) α n ,a = 0 . α n ,a = 0 . α n ,n = 0 . α a ,a = 0 . α a ,a = 0 . α a ,n = 0 . α a ,a = 0 . α a ,a = 0 . α a ,n = 0 . I ( a , n ) α a ,n = 0 . α a ,a = 0 . α n ,n = 1 . α n ,a = 1 . I ( a , n ) α a ,n = 0 . α a ,a = 0 . α n ,n = 1 . α n ,a = 1 . a) (b) (c)(d) (e) -3 (f) Figure 6: Numerical results of the spatially heterogeneous test case with parabolic config-uration of relaxation times and characteristic velocities ( τ = 10 − , λ = 10 ) and repro-duction number R >
1. Time and spatial evolution of S presented in (a), (b) and (c);evolution of I shown in (d), (e) and (f). In the first two examples, we consider a network composed by 3 nodes connected by 2bidirectional arcs in series, as shown in Fig. 7, having length L = L = 5 and discretizedwith a grid size ∆ x = 0 .
1. No social distancing or control effects are taken into account,fixing k = 0. With this network we want to simulate the spread of an epidemic that startsfrom node n (which represents a big city) and reaches nodes n (a city as big as the onerepresented by n ) and n (which identifies a rather small city). For this reason, initialconditions are null for each variable in the two arcs, while at nodes we fix S ( n ,
0) = 0 . , I ( n ,
0) = 0 . , R ( n ,
0) = 0 . ,S ( n ,
0) = 0 . , I ( n ,
0) = 0 . , R ( n ,
0) = 0 . ,S ( n ,
0) = 0 . , I ( n ,
0) = 0 . , R ( n ,
0) = 0 . , with zero initial fluxes for all the compartments. Transmission conditions at each arc-nodeinterface, satisfying conditions (14) and (15), are given in Table 5. It can be noticed thatit is imposed a symmetric transmission of the infectious on arcs a and a from n .The first simulation is run assuming that only the infectious individuals are movingalong the network, with λ I = 10 and τ I = 0 .
1, while susceptible and recovered subjectsdo not leave the origin node (city). This choice, albeit unrealistic, allows us to observe ina cleaner way the effects of transport and spread of the infection related to the mobility20 a n n n ( a ) ( a ) D en s i t y o f i nd i v i dua l s ( n ) D en s i t y o f i nd i v i dua l s ( n ) D en s i t y o f i nd i v i dua l s ( n ) Figure 7: Numerical results of the network test with 3 nodes connected by 2 arcs in serieswith a symmetric time and spatial dynamics of infectious individuals I in arc a and arc a and stationary susceptible and recovered people. Time evolution of the populationdynamics at node n , node n , from which the infectious disease start spreading, and atnode n .of people. In Fig. 7 we can notice a rapid decay of susceptible individuals at node n (from which the epidemic outbreak starts developing), which get ill because of the smallpercentage of infectious people that do not leave the node. Nevertheless, not the wholesusceptible population of the node is infected: almost 40% of it remains unharmed at theequilibrium state. A different evolution of the epidemic is observed at nodes n and n .Here it is clearly highlighted how the model associates a higher incidence function to thenode with a larger population. In fact, node n has twice as many susceptible subjectsas node n , information that is reproduced by the model as an epidemic that at n hasan incidence double that of n (even though infected subjects travel along arcs a and a in a total symmetry, as shown in Fig. 7). Indeed, the entire population at node n isinfected, while a very little percentage of susceptible at node n remains uninfected, evenif here individuals are the half of n . In addition, we can see that the epidemic peak itselfis reached much earlier at n .In the second simulation, we allow mobility also of susceptible and recovered subjects,to whom the same velocity and relaxation time of infectious individuals is assigned ( λ I =10, τ I = 0 . n in Fig. 7 to the one in Fig. 8. In this second simulation,indeed, hardly anyone remains at node n and the whole population moves symmetricallytowards nodes n and n . In these last nodes, the epidemic spreads in a similar way to21 a n n n ( a - S ) ( a - I ) D en s i t y o f i nd i v i dua l s ( n ) D en s i t y o f i nd i v i dua l s ( n ) D en s i t y o f i nd i v i dua l s ( n ) Figure 8: Numerical results of the network test with 3 nodes connected by 2 arcs in serieswith a symmetric time and spatial dynamics of susceptible S , infectious I and recovered R individuals in arcs a and a . Time evolution of the population dynamics at node n ,node n , from which the infectious disease starts spreading, and at node n .that seen in the previous scenario, with the difference that a percentage of the healthypopulation continues to arrive even when the disease is in regression. Furthermore, it isworth to be mentioned that in this case study, considering the network as a whole, lessindividuals are infected: almost 73% against the 86% of infected population obtained inthe previous simulation. This result is explained by the fact that the mobility added tothe susceptible in this case allows them to “escape” along the arcs, avoiding the contagionat the nodes of the network, especially the one at n . Clearly, this does not mean thatthe mobility of people helps the control of contagion, because in the first scenario onlysusceptible and recovered individuals are (unrealistically) not allowed to move, but all theinfectious are. In the third test case, we consider a more complex network, composed by 4 nodes and3 bidirectional arcs, as shown in Fig. 9. Arcs have length L = 2 and L = L = 5 andare discretized with a grid size ∆ x = 0 .
1. In this test, the epidemic outbreak starts from n (small city) and propagates along the whole network, reaching first n (metropolitancity) and then n and n (big cities). Thus, initial conditions, again null in all the arcs,22 a n n n n a ( a - S ) ( a - I )( a - S ) ( a - I ) D en s i t y o f i nd i v i dua l s ( n ) D en s i t y o f i nd i v i dua l s ( n ) D en s i t y o f i nd i v i dua l s ( n ) D en s i t y o f i nd i v i dua l s ( n ) Figure 9: Numerical results of the network test with 4 nodes connected by 3 arcs withbifurcation with time and spatial dynamics of susceptible S and infectious individuals inarcs a and a . Results concerning arc a are omitted because very similar to those in a .Time evolution of the population dynamics at node n , from which the infectious diseasestarts spreading, node n , node n and node n . Plot ( n ) is zoomed to exhibit also theevolution of infectious at that node.at nodes are imposed as follows: S ( n ,
0) = 0 . , I ( n ,
0) = 0 . , R ( n ,
0) = 0 . ,S ( n ,
0) = 0 . , I ( n ,
0) = 0 . , R ( n ,
0) = 0 . ,S ( n ,
0) = 0 . , I ( n ,
0) = 0 . , R ( n ,
0) = 0 . ,S ( n ,
0) = 0 . , I ( n ,
0) = 0 . , R ( n ,
0) = 0 . , with zero initial fluxes for all the classes. Transmission conditions at each interface aregiven in Table 6 and it can be seen that these coefficients define a more complex dynamicswith respect to the one of the previous tests. Here, we want to verify the effects of a23ontrol measure of social distancing on the population, imposing k = 1 at n , leaving k = 0 in the rest of the network. The same transport parameters are given to all thethree compartments, fixing again λ = 10 and τ = 0 .
1. In Fig. 9 numerical results ofthe evolution of the infectious disease at each node are presented. From the outbreaklocation, n , the population starts traveling through arc a , to reach node n . The formerreceives only later in time a return of part of the susceptible individuals who, after reaching n , come back to their origin point. The resulting mobility of susceptible and infectiouspeople along arc a is shown in the two initial plot of Fig. 9. At node n there is not areal epidemiological dynamics as it represents a node of transit, in which individuals donot stop. Waves of motion along arc a are represented in the third and fourth plot ofFig. 9. These waves are very similar to those observed in arc a , as a consequence of thetransmission coefficients chosen (reason why graphical results regarding a are omitted).There are very few more people moving in the direction of n than there are travelingtowards n , with both the nodes accommodating an initial susceptible population of equalsize. Nevertheless, the spread of the infectious disease is not comparable in these twonodes (cities), due to the control measure imposed at n ( k = 1) and not at n ( k = 0).In fact, it is here confirmed the impact of the enforcement of social distancing measures,which permits to significantly lower the epidemic peak at n . In this work, a novel SIR-type kinetic transport model for the spread of infectiousdiseases on networks is presented. The hyperbolic system describes at a macroscopic levelthe propagation of epidemics at finite speeds, recovering the classical one-dimensionalreaction-diffusion model as relaxation times and characteristic speeds of each compartmentof the population (susceptible, infectious and the recovered individuals) tend to zero andinfinity, respectively. The extension of the model to the treatment of networks that identifyat each node particular limited locations (like a specific city or region), connected bypaths (network arcs) along which individuals can move to reach different destinations, isalso presented. In this context, at each node the standard zero-dimensional SIR modeldescribes the evolution of the epidemic in the areas of major interest that are affected bythe mobility of the population in the network, which is computed in the arcs through thehere proposed SIR-type kinetic transport model. To ensure a correct coupling betweennodes and arcs, proper transmission conditions are defined at each arc-node interface,which guarantee the conservation of the global mass of the system.To solve the system of equations on each arc, a second-order IMEX finite volumemethod that is robust enough to correctly capture the asymptotic behavior of hyperbolicsystem under different kinds of scaling is proposed. In particular, the numerical methodsatisfies the AP property in the stiff regime, hence in the parabolic diffusive limit. Theexpected accuracy of the numerical scheme is confirmed for all the variables of the problemby means of accuracy analysis, even when dealing with stiff regimes, characterized byrelaxation times close to zero. Furthermore, the behavior of the model when consideringspatially heterogeneous environments is investigated concerning different configurations oftransport parameters (characteristic velocities and relaxation times). Results confirm that24he effects of a spatial variability of the contact rate vanish when the reproduction numberof the infectious disease, R , is less than 1, while a temporary persistence of the infectiousis noticed when R >
1. Finally, in order to demonstrate the suitability of the proposedmodel to simulate the spread of epidemic diseases on networks, three tests are carried outconcerning different simple networks. The impact of human mobility in these networksis assessed evaluating different mobility patterns, transmission coefficients at interfacesand social distancing interventions. Numerical results underline that these characteristicshighly influence the course of an epidemic, therefore confirming that restrictions on themobility of people and social distancing measures are very effective in reducing the spreadof an infectious disease.Since data of the spread epidemics are generally highly heterogeneous and affectedby a great deal of uncertainty, future perspectives include the application of uncertaintyquantification methods to assess the impact of stochastic inputs in the proposed SIR-typekinetic transport model. Furthermore, an extension of the model for the inclusion of theage structure of the population is foreseen, being an essential characteristic to correctlydescribe the impact of specific kinds of infectious diseases, like the COVID-19 pneumonia.
Acknowledgements
This work was partially supported by MIUR (Ministero dell’Istruzione, dell’Universit`ae della Ricerca) PRIN 2017, project “Innovative numerical methods for evolutionary partialdifferential equations and applications” , code 2017KKJP4X.
References [1] G. Albi, M. Zanella, and L. Pareschi. Control with uncertain data of socially struc-tured compartmental epidemic models.
ArXiv:2004.13067 , 2020.[2] D. Amadori and G. Guerra. Global weak solutions for systems of balance laws.
Applied Mathematics Letters , 12(6):123 – 127, 1999.[3] U. M. Ascher, S. J. Ruuth, and R. J. Spiteri. Implicit-explicit Runge-Kutta methodsfor time-dependent partial differential equations.
Applied Numerical Mathematics ,25(2-3):151–167, nov 1997.[4] D. Balcan, V. Colizza, B. Gon¸calves, H. Hud, J. J. Ramasco, and A. Vespignani. Mul-tiscale mobility networks and the spatial spreading of infectious diseases.
Proceedingsof the National Academy of Sciences of the United States of America , 106(51):21484–21489, 2009.[5] E. Barbera, G. Consolo, and G. Valenti. Spread of infectious diseases in a hy-perbolic reaction-diffusion susceptible-infected-removed model.
Physical Review E ,88(5):052719, nov 2013.[6] S. Bargmann, P. M. Jordan, and W. Lambert. A second-sound based hyperbolic SIRmodel for high-diffusivity spread.
Physics Letters A , 375(5):898–907, 2011.257] N. Bellomo, R. Bingham, M. A. J. Chaplain, G. Dosi, G. Forni, D. A. Knopoff,J. Lowengrub, R. Twarock, and M. E. Virgillito. A multi-scale model ofvirus pandemic: Heterogeneous interactive entities in a globally connected world. arXiv:2006.03915 , 2020.[8] S. Boscarino, L. Pareschi, and G. Russo. A unified IMEX Runge-Kutta approach forhyperbolic systems with multiscale relaxation.
SIAM Journal on Numerical Analysis ,55(4):2085–2109, 2017.[9] G. Bretti, R. Natalini, and B. Piccoli. Numerical approximations of a traffic flowmodel on networks.
Networks & Heterogeneous Media , 1(1):57–84, 2006.[10] G. Bretti, R. Natalini, and M. Ribot. A hyperbolic model of chemotaxis on a net-work: a numerical study.
ESAIM: Mathematical Modelling and Numerical Analysis ,48(1):231–258, jan 2014.[11] V. Capasso and G. Serio. A generalization of the Kermack-McKendrick deterministicepidemic model.
Mathematical Biosciences , 42(1-2):43–61, nov 1978.[12] R. M. Colombo, M. Garavello, F. Marcellini, and E. Rossi. An age and space struc-tured SIR model describing the Covid-19 pandemic. medRxiv , 2020.[13] L. Fermo and A. Tosin. A fully-discrete-state kinetic theory approach to traffic flowon road networks.
Math. Models Methods Appl. Sci. , 25(3):423–461, 2015.[14] E. Franco. A feedback SIR (fSIR) model highlights advantages and limitations ofinfection-based social distancing. arXiv:2004.13216 , 2020.[15] K. O. Friedrichs and P. D. Lax. Systems of Conservation Equations with a ConvexExtension.
Proceedings of the National Academy of Sciences , 68(8):1686–1688, aug1971.[16] M. Gatto, E. Bertuzzo, L. Mari, S. Miccoli, L. Carraro, R. Casagrandi, and A. Ri-naldo. Spread and dynamics of the COVID-19 epidemic in Italy: Effects of emer-gency containment measures.
Proceedings of the National Academy of Sciences ,117(19):10484–10491, may 2020.[17] G. Giordano, F. Blanchini, R. Bruno, P. Colaneri, A. Di Filippo, A. Di Matteo, andM. Colaneri. Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy.
Nature Medicine , 26:855–860, 2020.[18] L. Gosse and G. Toscani. An asymptotic-preserving well-balanced scheme for thehyperbolic heat equations.
C. R. Math. Acad. Sci. Paris , 334(4):337–342, 2002.[19] H. W. Hethcote. The Mathematics of Infectious Diseases.
SIAM Review , 42(4):599–653, jan 2000.[20] S. Jin, L. Pareschi, and G. Toscani. Diffusive relaxation schemes for multiscalediscrete-velocity kinetic equations.
SIAM J. Numer. Anal. , 35(6):2405–2439, 1998.2621] A. Korobeinikov and P. K. Maini. Non-linear incidence and stability of infectious dis-ease models.
Mathematical Medicine and Biology: A Journal of the IMA , 22(2):113–128, jun 2005.[22] M. U. G. Kraemer, C.-H. Yang, B. Gutierrez, C.-H. Wu, B. Klein, D. M. Pigott,L. du Plessis, N. R. Faria, R. Li, W. P. Hanage, J. S. Brownstein, M. Layan, A. Vespig-nani, H. Tian, C. Dye, O. G. Pybus, and S. V. Scarpino. The effect of human mobilityand control measures on the COVID-19 epidemic in China.
Science , 368(6490):493–497, may 2020.[23] P. L. Lions and G. Toscani. Diffusive limit for finite velocity Boltzmann kineticmodels.
Rev. Mat. Iberoamericana , 13(3):473–513, 1997.[24] S. Merler and M. Ajelli. Human mobility and population heterogeneity in the spreadof an epidemic.
Procedia Computer Science , 1(1):2237–2244, 2010.[25] I. M¨uller and T. Ruggeri.
Rational Extended Thermodynamics . Springer, New York,1998.[26] J. D. Murray.
Mathematical Biology I, II . Springer- Verlag, New York, 2002-2003.[27] G. Naldi and L. Pareschi. Numerical schemes for hyperbolic systems of conservationlaws with stiff diffusive relaxation.
SIAM J. Numer. Anal. , 37(4):1246–1270, 2000.[28] L. Pareschi and G. Russo. Implicit-explicit Runge-Kutta schemes and applicationsto hyperbolic systems with relaxation.
Journal of Scientific Computing , 25(1/2):129–155, 2005.[29] L. Pellis, F. Ball, S. Bansal, K. Eames, T. House, V. Isham, and P. Trapman. Eightchallenges for network epidemic models.
Epidemics , 10:58–62, mar 2015.[30] B. Piccoli and M. Garavello.
Traffic Flow on Networks . American Institute of Math-ematical Sciences, 2006.[31] S. Riley, K. Eames, V. Isham, D. Mollison, and P. Trapman. Five challenges forspatial epidemic models.
Epidemics , 10(2015):68–71, mar 2015.[32] G.-Q. Sun. Pattern formation of an epidemic model with diffusion.
Nonlinear Dy-namics , 69(3):1097–1104, aug 2012.[33] E. F. Toro.
Riemann Solvers and Numerical Methods for Fluid Dynamics . SpringerVerlag, 3rd edition, 2009.[34] J. Wang, F. Xie, and T. Kuniya. Analysis of a reaction-diffusion cholera epidemicmodel in a spatially heterogeneous environment.
Communications in Nonlinear Sci-ence and Numerical Simulation , 80:104951, jan 2020.[35] Y. Wang, J. Wang, and L. Zhang. Cross diffusion-induced pattern in an SI model.